算額（その1408）

算額（その1408）

百四十九 群馬県多野郡新町 稲荷神社 文政3年(1820)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，半円，直角三角形
#Julia, #SymPy, #算額, #和算

半円の中に三角形（直角三角形）を作り，甲円，乙円，丙円，丁円を容れる。大斜（半円の直径），小斜（三角形の一番短い辺）が与えられたとき，丁円の直径はいかほどか。

半円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (x4, r4)
円周上にある直角三角形の頂点の座標を (x0, sqrt(R^2 - x0))
とおき，以下の連立方程式を解く。

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

@syms R, x0, r1, x1, r2, x2, r3, x3, r4, x4,
      大斜, 小斜, 中斜
R = 大斜/2
eq1 = sqrt((R - x0)^2 + (R^2 - x0^2)) - 小斜
eq2 = r2/(x2 + R) - r1/(x1 + R)
eq3 = r3/(R - x3) - r1/(R - x1)
eq4 = x3 - x1 - 2sqrt(r1*r3)
eq5 = x1 - x4 - 2sqrt(r1*r4)
eq6 = x4 - x2 - 2sqrt(r2*r4)
eq7 = x1 - x2 - 2sqrt(r1*r2)
eq8 = dist2(-R, 0, x0, sqrt(R^2 - x0^2), x1, r1, r1)
eq9 = 小斜 + sqrt(大斜^2 - 小斜^2) - 大斜 - 2r1;
# eqxx = (大斜 + sqrt(大斜^2 - 小斜^2) + 小斜)*r1 - 小斜*中斜
# eqxx = dist2(R, 0, x0, sqrt(R^2 - x0^2), x1, r1, r1);
# eqxx = dist2(-R, 0, x0, sqrt(R^2 - x0^2), x2, r2, r2)
# eqxx = dist2(R, 0, x0, sqrt(R^2 - x0^2), x3, r3, r3);
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq10], (x0, r1, x1, r2, x2, r3, x3, r4, x4))

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)
    (x0, r1, x1, r2, x2, r3, x3, r4, x4) = u
    return [
        -小斜 + sqrt(-x0^2 + 大斜^2/4 + (-x0 + 大斜/2)^2),  # eq1
        -r1/(x1 + 大斜/2) + r2/(x2 + 大斜/2),  # eq2
        -r1/(-x1 + 大斜/2) + r3/(-x3 + 大斜/2),  # eq3
        -x1 + x3 - 2*sqrt(r1*r3),  # eq4
        x1 - x4 - 2*sqrt(r1*r4),  # eq5
        -x2 + x4 - 2*sqrt(r2*r4),  # eq6
        x1 - x2 - 2*sqrt(r1*r2),  # eq7
        (r1^2*x0 - r1*x1*sqrt(-4*x0^2 + 大斜^2) - x0*x1^2 + 大斜*(-4*r1^2 - 4*r1*sqrt(-4*x0^2 + 大斜^2) - 8*x0*x1 - 2*x0*大斜 + 4*x1^2 + 4*x1*大斜 + 大斜^2)/8)/大斜,  # eq8
        -2*r1 - 大斜 + 小斜 + sqrt(大斜^2 - 小斜^2),  # eq9
    ]
end;

大斜 = 10
R = 大斜/2
小斜 = 4
iniv = BigFloat[3.45, 1.58, 2.58, 1.05, 0.01, 0.46, 4.29, 0.32, 1.16]
res = nls(H, ini=iniv)

    ([3.4, 1.5825756949558398, 2.582575694955839, 1.0456116707250647, 0.009826491126086901, 0.46246227038447935, 4.293577213300662, 0.31816569849213683, 1.163390997458993], true)

たとえば，大斜，小斜，中斜が 10, 4, 9.16515 のとき，丁円の直径は 0.6363313969842737 である。

全てのパラメータは以下のとおりである。

   x0 = 3.4;  r1 = 1.58258;  x1 = 2.58258;  r2 = 1.04561;  x2 = 0.00982649;  r3 = 0.462462;  x3 = 4.29358;  r4 = 0.318166;  x4 =1.16339

function draw(大斜, 小斜, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 大斜/2
    中斜 = sqrt(大斜^2 - 小斜^2)
    (x0, r1, x1, r2, x2, r3, x3, r4, x4) = res[1]
    @printf("大斜，小斜，中斜が %g, %g, %g のとき，丁円の直径は %g である。\n", 大斜, 小斜, 中斜, 2r4)
    @printf("x0 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x4 =%g\n", x0, r1, x1, r2, x2, r3, x3, r4, x4)
    y0 = sqrt(R^2 - x0^2)
    plot([-R, x0, R], [0, y0, 0], color=:blue, lw=0.5)
    circle(0, 0, R, beginangle=0, endangle=180)
    circle(x1, r1, r1)
    circle(x2, r2, r2, :green)
    circle(x3, r3, r3, :magenta)
    circle(x4, r4, r4, :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(x0, sqrt(R^2 - x0^2), "(x0,sqrt(R^2-x0^2))", :blue, :center, :bottom, delta=delta)
        point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
        point(x2, r2, "乙円:r2\n(x2,r2)", :green, :center, delta=-delta)
        point(x3, r3, "丙円:r3,(x3,r3)", :magenta, :right, delta=-7delta, deltax=13delta)
        point(x4, r4, "丁円:r4,(x4,r4)", :orange, :right, delta=-5delta, deltax=13delta)
        point(R, 0, "大斜/2", :red, :left, :bottom, delta=delta)
        point(0, 2.6r2, "中斜", :blue, :center, mark=false)
        point((R + x0)/2, r1, "小斜", :blue, :center, mark=false, deltax=7delta)
        plot!(xlims=(-R - 3delta, R + 15delta), ylims!(-10delta, R + 3delta))
    end
end;

draw(10, 4, true)