手打ちうどん 飩兵衛

高松市川部町
しっぽくうどんをいただく
温まります

算額（その1445）

算額（その1445）

福島県二本松市 若宮八幡神社 文化5年(1808) 復元奉納平成16年(2004)
http://www.wasan.jp/fukusima/nihonmatuhatiman.html

街角の数学 Street Wasan　～落書き帳「○△□」～ 185．十二歳の少年、算額を奉納す
http://streetwasan.web.fc2.com/math16.9.21.html

キーワード：球9個，外球，3次元
#Julia, #SymPy, #算額, #和算

大球の中に甲球 2 個，小球数個を容れる。小球は互いに外接し，甲球と外接し，大球に内接する。小球の直径が 1 寸のとき，大球の直径はいかほどか。

上から見たのが左図。小球は互いに外接し，大球に内接している。
横から見たのが右図。小球は甲球にも外接している。

大球の半径と中心座標を R, (0, 0)
甲球の半径と中心座標を r1, (0, 0, r1), (0, 0, -r1)
小球の半径と中心座標を r2, (R - r2, 0, 0)
とおき，以下の連立方程式を解く。

小球の個数は 6 個である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive,
      θ::positive, n::positive
n = 6
θ = 180/n
r1 = R/2
x2 = r2/sind(θ)
eq1 = (R - r2)^2 + r1^2 - (r1 + r2)^2
res = solve(eq1, R)[1]
res |> println

    3*r2

大球の半径は，小球の半径の 3 倍である。
小球の直径が 1 寸のとき，大球の直径は 3 寸である。

r2 = 1/2
R = 3r2
r1 = R/2
n = 6
θ = 180/n
x2 = r2/sind(θ)
(R, r1, x2) |> println

    (1.5, 0.75, 1.0)

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 3r2
    r1 = R/2
    n = 6
    θ = 180/n
    x2 = r2/sind(θ)
    @printf("小球の直径が %g のとき，大球の直径は %g，甲球の直径は %g である。\n", 2r2, 2R, 2r1)
    p1 = plot()
    circle(0, 0, R, :orange)
    circle(0, 0, r1, :blue)
    rotate(x2, 0, r2, angle=2θ)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(x2, 0, "小球:r2,(x2,0)", :red, :center, delta=-delta)
        point(0, R, "大球:R", :orange, :center, :bottom, delta=delta)
        point(0, 0, "甲球:r1", :blue, :center, :bottom, delta=delta)
    end
    p2 = plot()
    circle(0, 0, R, :orange)
    circle22(0, 0.75, 0.75, :blue)
    circle2(1, 0, 0.5)
    circle2(0.5, 0, 0.5)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(x2, 0, "小球:r2,(x2,0)", :red, :center, delta=-delta)
        point(0, R, "大球:R", :orange, :center, :bottom, delta=delta)
        point(0, r1, "甲球:r1", :blue, :center, delta=-delta)
    end
    plot!(p1, p2)
end;

draw(1/2, true)

算額（その1444）

算額（その1444）

福島県福島市堂殿 黒沼神社 明治11年(1878)
街角の数学 Street Wasan　～落書き帳「○△□」～ 186．僕は十一歳
http://streetwasan.web.fc2.com/math16.9.22.html
キーワード：円6個，円弧
#Julia, #SymPy, #算額, #和算

円弧（弓形）の中に，甲円 2 個，乙円 2 個，丙円 3 個を容れる。丙円の直径が 1 寸のとき，甲円の直径はいかほどか。

円弧の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, y + r1)
乙円の半径と中心座標を r2, (0, R - r2)
丙円の半径と中心座標を r3, (0, y + r3), (x3, y + r3)
y = R - 2r2 - 2r3
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r1::positive, x1::positive, r2::positive,
      r3::positive, x3::positive
y = R - 2r2 - 2r3
x3 = 2x1
eq1 = x1^2 + (y + r1)^2 - (R - r1)^2
eq2 = x3^2 + (y + r3)^2 - (R - r3)^2
eq3 = x1^2 + (2r3 + r2 - r1)^2 - (r1 + r2)^2
eq4 = x1^2 + (r1 - r3)^2 - (r1 + r3)^2
res = solve([eq1, eq2, eq3, eq4], (r1, r2, R, x1))[1];

# r1 甲円の半径
res[1] |> println
res[1](r3 => 1/2).evalf() |> println

    r3*(sqrt(3) + 2)/2
    0.933012701892219

# r2 乙円の半径
res[2] |> println
res[2](r3 => 1/2).evalf() |> println

    2*sqrt(3)*r3/3
    0.577350269189626

# R 円弧の半径
res[3] |> println
res[3](r3 => 1/2).evalf() |> println

    4*r3*(3 + 2*sqrt(3))/3
    4.30940107675850

# x1 甲円の中心の x 座標
res[4] |> println
res[4](r3 => 1/2).evalf() |> println

    r3*(1 + sqrt(3))
    1.36602540378444

丙円の直径が 1 寸のとき，甲円の直径は 1.86602540378444 寸である。

「答曰、甲円径一寸八分有奇」と，切り捨てたであるが，「術曰、置七分五厘開平方加一個乗丙径、得甲径合問」で「甲円直径 = (sqrt(0.75) + 1)*丙円直径」と正確である。

ちなみに，乙円の直径は 1.15470053837925 寸である。

---

和算家の取り組みと違い，連立方程式を解くという方法は，連立方程式はそのままで求めるものを変えるだけで良いということがある。

・ 甲円の半径を与えて，乙円，丙円，円弧の半径を求めるときは以下のようになる。

res2 = solve([eq1, eq2, eq3, eq4], (r2, r3, R, x1))[1];
res2[1] |> println
res2[2] |> println
res2[3] |> println

    4*r1*(-3 + 2*sqrt(3))/3
    2*r1*(2 - sqrt(3))
    8*sqrt(3)*r1/3

・ 乙円の半径を与えて，甲円，丙円，円弧の半径を求めるときは以下のようになる。

res3 = solve([eq1, eq2, eq3, eq4], (r1, r3, R, x1))[1];
res3[1] |> factor |> println
res3[2] |> println
res3[3] |> factor |> println

    r2*(3 + 2*sqrt(3))/4
    sqrt(3)*r2/2
    2*r2*(sqrt(3) + 2)

・ 円弧の半径を与えて，甲円，乙円，丙円の半径を求めるときは以下のようになる。

res4 = solve([eq1, eq2, eq3, eq4], (r1, r2, r3, x1))[1];
res4[1] |> println
res4[2] |> factor |> println
res4[3] |> factor |> println

    sqrt(3)*R/8
    -R*(-2 + sqrt(3))/2
    R*(-3 + 2*sqrt(3))/4

・ なんなら，甲円の中心の x 座標を与えて，甲円，乙円，丙円，円弧の半径を求めるときは以下のようになる。

res5 = solve([eq1, eq2, eq3, eq4], (r1, r2, r3, R))[1];
res5[1] |> println
res5[2] |> factor |> println
res5[3] |> factor |> println
res5[4] |> factor |> println

    x1*(1 + sqrt(3))/4
    -x1*(-3 + sqrt(3))/3
    x1*(-1 + sqrt(3))/2
    2*x1*(sqrt(3) + 3)/3

function draw(r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, r1, x1, r2) = (4r3*(3 + 2√3)/3, r3*(√3 + 2)/2, r3*sqrt(2√3 + 4), 2√3r3/3)
    @printf("丙円の直径が %g のとき，甲円の直径は %g である。\n", 2r3, 2r1)
    @printf("r3 = %g;  R = %g;  r1 = %g;  x1 = %g;  r2 = %g\n", r3, R, r1, x1, r2)
    y = R - 2r2 - 2r3
    x = sqrt(R^2 - y^2)
    θ = atand(y, x)
    x3 = 2x1
    plot()
    circle(0, 0, R, beginangle=θ, endangle=180-θ)
    circle2(x1, y + r1, r1, :blue)
    circle(0, R - r2, r2, :green)
    circle(0, y + r3, r3, :magenta)
    circle2(x3, y + r3, r3, :magenta)
    segment(-x, y, x, y, :orange)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(x1, y + r1, "甲円:r1\n(x1,y+r1)", :blue, :center, delta=-3delta)
        point(0, R - r2, "乙円:r2,(0,R-r2)", :green, :left, delta=-3delta, deltax=-15delta)
        point(0, y + r3, "丙円:r3,(0,y+r3)", :magenta, :right, delta=-3delta, deltax=12delta)
        point(x3, y + r3, "丙円:r3,(x3,y+r3)", :magenta, :right, delta=-3delta, deltax=12delta)
        point(0, R, "R", :red, :center, :bottom, delta=3delta)
        point(0, y, "y", :red, :center, delta=-3delta)
    end
end;

draw(1/2, true)

 

算額（その1443）

算額（その1443）

街角の数学 Street Wasan　～落書き帳「○△□」～ 934. 『算法天生指南』巻之二（その９）
http://streetwasan.web.fc2.com/math20.04.29.3.html
キーワード：円4個，直線上
#Julia, #SymPy, #算額, #和算

問題 III. 直線上に 4 個の円が載っている。甲円，乙円，丁円の直径が 72 寸，24 寸，32 寸のとき，丙円の直径はいかほどか。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (0, r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, x1::positive, r2::positive, x2::positive,
      r3::positive, x3::positive, y3::positive, r4::positive
# x2 = 2sqrt(r2*r4)
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (y3 - r2)^2 - (r2 + r3)^2
eq4 = x2^2 + (r4 - r2)^2 - (r2 + r4)^2
eq5 = x3^2 + (y3 - r4)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (x1, x2, r3, x3, y3))[3];  # 3 of 4

 それぞれは，長い式になっているので，以下のようにして簡約化する。

1. 「変数^(n/2)」は，r1^(1//2) => s, r4^(1//2) => t のように，変数変換を行う。（このようにしないと因数分解できない）
2. 最後に s => r1^(1//2), t => r4^(1//2) と変数変換し，もとに戻す。

1.　x1

変数変換，因数分解，簡約化を経て，以下のようになる。

@syms s, t
ans_x1 = res[1](r1^(1//2) => s, r4^(1//2) => t) |> factor

長い式の平方根は，式を二乗して因数分解すると二乗の分数式になるので，平方根をとって (r2*s + r2*t - 2s*t^2)/(r2 - s*t) である。

ans_x1.args[3]^2 |> factor |> println

    (r2*s + r2*t - 2*s*t^2)^2/(r2 - s*t)^2

根号の外と，分母に同じ項を含むので約分して以下のように簡単になる。

ans_x1 = 2sqrt(r2)* (r2*s + r2*t - 2s*t^2)/(r2 - s*t) *(r2 - s*t)*(s + t)/(r2*s + r2*t- 2s*t^2)
ans_x1 |> println

    2*sqrt(r2)*(s + t)

必要ならば，変数変換をもとに戻す。

ans_x1 = ans_x1(s => r1^(1//2), t => r4^(1//2))
ans_x1 |> println

    2*sqrt(r2)*(sqrt(r1) + sqrt(r4))

簡約化の前後で同値な式であることを確認する。

res[1](r1 => 72/2, r2 => 24/2, r4 => 32/2) |> println
ans_x1(r1 => 72/2, r2 => 24/2, r4 => 32/2) |> println

    69.2820323027551
    69.2820323027551

2.　x2

@syms s, t
ans_x2 = res[2](r1^(1//2) => s, r4^(1//2) => t) |> factor

長い式の平方根は x1 と同じ式 (r2*s + r2*t - 2s*t^2)/(r2 - s*t) であるから，以下のように簡約化される。

ans_x2 = 2sqrt(r2)*t* (r2*s + r2*t - 2s*t^2)/(r2 - s*t) *(r2 - s*t)/(r2*s + r2*t- 2s*t^2)
ans_x2 |> println

    2*sqrt(r2)*t

ans_x2 = ans_x2(s => r1^(1//2), t => r4^(1//2))
ans_x2 |> println

    2*sqrt(r2)*sqrt(r4)

res[2](r1 => 36, r2 => 12, r4 => 16).evalf() |> println
ans_x2(r1 => 36, r2 => 12, r4 => 16).evalf() |> println

    27.7128129211020
    27.7128129211020

3.　r3

@syms s, t
ans_r3 = res[3](r1^(1//2) => s, r4^(1//2) => t) |> factor
ans_r3 |> println

    -r2^2*(s + t)^2/(4*s*t*(r2 - s*t))

ans_r3 = ans_r3(s => r1^(1//2), t => r4^(1//2))
ans_r3 |> println

    -r2^2*(sqrt(r1) + sqrt(r4))^2/(4*sqrt(r1)*sqrt(r4)*(-sqrt(r1)*sqrt(r4) + r2))

res[3](r1 => 36, r2 => 12, r4 => 16).evalf() |> println
ans_r3(r1 => 36, r2 => 12, r4 => 16).evalf() |> println

    12.5000000000000
    12.5000000000000

4.　x3

@syms s, t
ans_x3 = res[4](r1^(1//2) => s, r4^(1//2) => t) |> factor

長い式の平方根は x1 と同じ式であるから，以下のように簡約化される。

ans_x3 = sqrt(r2)* (r2*s + r2*t - 2s*t^2)/(r2 - s*t)
ans_x3 |>  println

    sqrt(r2)*(r2*s + r2*t - 2*s*t^2)/(r2 - s*t)

ans_x3 = ans_x3(s => r1^(1//2), t => r4^(1//2))
ans_x3 |> println

    sqrt(r2)*(sqrt(r1)*r2 - 2*sqrt(r1)*r4 + r2*sqrt(r4))/(-sqrt(r1)*sqrt(r4) + r2)

res[4](r1 => 36, r2 => 12, r4 => 16).evalf() |> println
ans_x3(r1 => 36, r2 => 12, r4 => 16).evalf() |> println

    20.7846096908265
    20.7846096908265

5.　y3

@syms s, t
ans_y3 = res[5](r1^(1//2) => s, r4^(1//2) => t) |> factor
ans_y3 |> println

    r2*(r2*s^2 + 2*r2*s*t + r2*t^2 - 8*s^2*t^2)/(4*s*t*(r2 - s*t))

ans_y3 = ans_y3(s => r1^(1//2), t => r4^(1//2))
ans_y3 |> println

    r2*(2*sqrt(r1)*r2*sqrt(r4) + r1*r2 - 8*r1*r4 + r2*r4)/(4*sqrt(r1)*sqrt(r4)*(-sqrt(r1)*sqrt(r4) + r2))

res[5](r1 => 36, r2 => 12, r4 => 16).evalf() |> println
ans_y3(r1 => 36, r2 => 12, r4 => 16).evalf() |> println

    35.5000000000000
    35.5000000000000

一時変数を設定すると若干短く定義できる。
    s = sqrt(r1) + sqrt(r4)
    t = (r2 - sqrt(r1*r4))
    x1 = 2sqrt(r2)*s
    x2 = 2sqrt(r2*r4)
    r3 = -r2^2*s^2/(4sqrt(r1*r4)*t)
    x3 = sqrt(r2)*(sqrt(r1)*r2 - 2sqrt(r1)*r4 + r2*sqrt(r4))/t
    y3 = r2*(2sqrt(r1*r4)*r2 + r1*r2 - 8r1*r4 + r2*r4)/(4sqrt(r1*r4)*t)

function draw(r1, r2, r4, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    s = sqrt(r1) + sqrt(r4)
    t = (r2 - sqrt(r1*r4))
    x1 = 2sqrt(r2)*s
    x2 = 2sqrt(r2*r4)
    r3 = -r2^2*s^2/(4sqrt(r1*r4)*t)
    x3 = sqrt(r2)*(sqrt(r1)*r2 - 2sqrt(r1)*r4 + r2*sqrt(r4))/t
    y3 = r2*(2sqrt(r1*r4)*r2 + r1*r2 - 8r1*r4 + r2*r4)/(4sqrt(r1*r4)*t)
    @printf("r1 = %g;  r2 = %g;  r4 = %g;  r3 = %g\n", r1, r2, r4, r3)
    @printf("x1 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", x1, x2, r3, x3, y3)
    plot()
    circle(x1, r1, r1)
    circle(x2, r2, r2, :blue)
    circle(x3, y3, r3, :green)
    circle(0, r4, r4, :magenta)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
        point(x2, r2, "乙円:r2\n(x2,r2)", :blue, :center, delta=-delta)
        point(x3, y3, "丙円:r3\n(x3,y3)", :green, :center, delta=-delta)
        point(0, r4, "丁円:r4,(0,r4)", :magenta, :center, delta=-delta)
    end
end;

draw(72/2, 24/2, 32/2, true)