算額（その1420）

算額（その1420）

福島県二本松市亀谷 坂上観音堂 寛政5年(1793)令和3年復元奉納
http://www.wasan.jp/fukusima/sakauekanondo.html
キーワード：円5個，円弧
#Julia, #SymPy, #算額, #和算

円弧（弓形）の中に 5 個の円を入れる。そのうちの甲円，乙円，丙円の直径が 16 寸，25 寸，9 寸のとき，外円の直径はいかほどか。

甲円，乙円，丙円が右から順に一つおきになっているが，残りの円も右から丁円，丙円と名付ける。

外円（円弧，弓形を構成する円）の半径と中心座標を R, (0, 0)
円弧を構成する水平な弦と y 軸の交点座標を (0, y)
甲円の半径と中心座標を r1, (x1, y + r1)
乙円の半径と中心座標を r2, (x2, y + r2)
丙円の半径と中心座標を r3, (x3, y + r3)
丁円の半径と中心座標を r4, (x4, y + r4)
戊円の半径と中心座標を r5, (x5, y + r5)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms R, y, r1, x1, r2, x2, r3, x3, r4, x4, r5, x5
eq1 = x1^2 + (y + r1)^2 - (R - r1)^2
eq2 = x2^2 + (y + r2)^2 - (R - r2)^2
eq3 = x3^2 + (y + r3)^2 - (R - r3)^2
eq4 = x4^2 + (y + r4)^2 - (R - r4)^2
eq5 = x5^2 + (y + r5)^2 - (R - r5)^2
eq6 = (x1 - x4)^2 +(r1 - r4)^2 - (r1 + r4)^2
eq7 = (x4 - x2)^2 +(r4 - r2)^2 - (r4 + r2)^2
eq8 = (x2 - x5)^2 +(r2 - r5)^2 - (r2 + r5)^2
eq9 = (x5 - x3)^2 +(r5 - r3)^2 - (r5 + r3)^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (R, y, x1, x2, x3, r4, x4, r5, x5))

using NLsolve

function nls(func, params...; ini = [0.0])
    if typeof(ini) <: Number
        r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
        v = r.zero[1]
    else
        r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
        v = r.zero
    end
    return Float64.(v), r.f_converged
end;

function H(u)
    (R, y, x1, x2, x3, r4, x4, r5, x5) = u
    return [
        x1^2 - (R - r1)^2 + (r1 + y)^2,  # eq1
        x2^2 - (R - r2)^2 + (r2 + y)^2,  # eq2
        x3^2 - (R - r3)^2 + (r3 + y)^2,  # eq3
        x4^2 - (R - r4)^2 + (r4 + y)^2,  # eq4
        x5^2 - (R - r5)^2 + (r5 + y)^2,  # eq5
        (r1 - r4)^2 - (r1 + r4)^2 + (x1 - x4)^2,  # eq6
        (-r2 + r4)^2 - (r2 + r4)^2 + (-x2 + x4)^2,  # eq7
        (r2 - r5)^2 - (r2 + r5)^2 + (x2 - x5)^2,  # eq8
        (-r3 + r5)^2 - (r3 + r5)^2 + (-x3 + x5)^2,  # eq9
    ]
end;
(r1, r2, r3) = (16, 25, 9) ./ 2
iniv = BigFloat[70, 45, 33, -12, -43, 12, 13, 8, -31]
res = nls(H, ini=iniv)

    ([69.73157051282051, 43.72996794871795, 33.686751312703706, -10.660364339463198, -43.92070107858837, 12.139917695473251, 13.976922133962859, 8.642578125, -31.448074801416432], true)

# 直径の分数表示
rationalize(2*69.73157051282051)

    87025//624

# 帯分数で表示するための計算 商と余りの計算
divrem(87025,624)

    (139, 289)

外円の直径は 2*69.73157051282051 = 139.46314102564102 であり，分数で表すと 87025//624 である。
帯分数表示というのは小学校以来使うことはないが，139 + 289/624 であり，「答」の「外円径一百三十九寸六百二十四分之二百八十九寸とピタリ一致する（長精度演算のおかげ）。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2, r3) = (16, 25, 9) ./ 2
    (R, y, x1, x2, x3, r4, x4, r5, x5) = [70, 45, 33, -12, -43, 12, 13, 8, -31]
    (R, y, x1, x2, x3, r4, x4, r5, x5) = [69.73157051282051, 43.72996794871795, 33.686751312703706, -10.660364339463198, -43.92070107858837, 12.139917695473251, 13.976922133962859, 8.642578125, -31.448074801416432]
    x =sqrt( R^2 - y^2)
    θ = atand(y, x)
    plot()
    circle(0, 0, R, beginangle=θ, endangle=180 - θ, n=500)
    circle(x1, y + r1, r1, :blue)
    circle(x2, y + r2, r2, :magenta)
    circle(x3, y + r3, r3, :green)
    circle(x4, y + r4, r4, :orange)
    circle(x5, y + r5, r5, :tomato)
    plot!([-x, 0, x, -x], [y, 0, y, y], color=:gray70, lw=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(x1, y + r1, "甲円:r1,(x1,y+r1)", :blue, :center, delta=-8.5delta)
        point(x2, y + r2, "乙円:r2\n(x2,y+r2)", :magenta, :center, delta=-delta)
        point(x3, y + r3, "丙円:r3,(x3,y+r3)", :green, :left, delta=-5delta, deltax=-5delta)
        point(x4, y + r4, "丁円:r4\n(x4,y+r4)", :orange, :center, delta=-delta)
        point(x5, y + r5, "戊円:r5\n(x5,y+r5)", :tomato, :center, :bottom, delta=delta)
        point(0, y, " y", :black, :left, delta=-delta)
        point(0, R, "R", :red, :center, :bottom, delta=delta)
    end
end;

draw(true)