算額（その1421）

算額（その1421）

算額（その436）では数値解を求めたが，解析解を求めた記事である。

橘田彌曾八元克 天明八年戊申2月
藤田貞資(1789)：神壁算法
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
キーワード：円7個，外円，弦
#Julia, #SymPy, #算額, #和算

外円の中に，甲円 3 個，乙円 2 個，丙円 1 個が入っている。甲円は弦に接している。
外円の直径が 3 寸 6 分 のとき，甲円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (2r1, a + r1)
乙円の半径と中心座標を r2, (r1, y2)
丙円の半径と中心座標を r3, (0, r0 - r3)
として以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, a::positive, r1::negative,
      r3::positive, r2::positive, y2::positive;
@syms r0, a, r1, r3, r2, y2
a = r0 - 2r1 - 2r3
eq1 = (2r1)^2 + (a + r1)^2 - (r0 - r1)^2 |> expand
eq2 = r1^2 + (y2 - a - r1)^2 - (r1 + r2)^2 |> expand
eq3 = r1^2 + y2^2 - (r0 - r2)^2 |> expand
eq4 = r1^2 + (r0 - r3 - y2)^2  - (r3 + r2)^2 |> expand;

一度に解けないので，まず eq1, eq2, eq3 から r2, y2, r3 を求める（求める変数の指定順序により，得られる式の複雑さが異なることがある）。
4 組の解が得られるが，最初のものが適解である。

res = solve([eq1, eq2, eq3], (r2, y2, r3))[1];

# r2
res[1] |> println

    (r0^2 - r1^2 - sqrt(r0*(r0 + r1))*sqrt(r0^2 - 2*r0*r1 - 3*r1^2))/(2*(r0 + r1))

# y2
res[2] |> println

    sqrt(r0*(r0 + r1))/2 + sqrt(r0^2 - 2*r0*r1 - 3*r1^2)/2

# r3
res[3] |> println

    r0/2 - r1/2 - sqrt((r0 - 3*r1)*(r0 + r1))/2

得られた解を eq4 に代入した式を eq14 とする。eq14 は既知の r0，未知の r1 のみを含む式になる。

eq14 = eq4(r2 => res[1], y2 => res[2], r3 => res[3]) |> simplify |> numerator
eq14 |> println

    -r0^3 + 3*r0^2*r1 - r0^2*sqrt(r0*(r0 + r1)) + r0^2*sqrt(r0^2 - 2*r0*r1 - 3*r1^2) + 5*r0*r1^2 + r0*sqrt(r0*(r0 + r1))*sqrt(r0^2 - 2*r0*r1 - 3*r1^2) + r1^3 + r1^2*sqrt(r0*(r0 + r1)) - r1^2*sqrt(r0^2 - 2*r0*r1 - 3*r1^2) - r1*sqrt(r0*(r0 + r1))*sqrt(r0^2 - 2*r0*r1 - 3*r1^2)

eq14 を解いて r1 を得る。

res2 = solve(eq14, r1);

5 個の解が得られるが，最後のものが適解である。

# r1
res2[5] |> println

    r0*(-(8 + 24*sqrt(57))^(1/3)/3 + 1/3 + 32/(3*(8 + 24*sqrt(57))^(1/3)))

res2[5](r0 => 3.6/2).evalf() |> println

    0.500028910096869

甲円の半径 r0 は，外円の半径 r0 の (1/3　- (8 + 24√57)^(1/3)/3 + 32/(3*(8 + 24√57)^(1/3))) 倍である。
外円の直径が 3.6 寸のとき，甲円の直径は 1.000057820193738 寸である。

他のパラメータは以下のようにして芋づる式に計算できる。
下のプログラムの draw 関数内では，一時変数を使ってより短い計算式にしている。
計算順序により有効数字14，5桁以降に若干の誤差が含まれる。

r0 = 3.6/2
r1 = r0*(1/3　- (8 + 24√57)^(1/3)/3 + 32/(3*(8 + 24√57)^(1/3)))
r2 = (r0^2 - r1^2 - sqrt(r0*(r0 + r1))*sqrt(r0^2 - 2*r0*r1 - 3*r1^2))/(2*(r0 + r1))
y2 = sqrt(r0*(r0 + r1))/2 + sqrt(r0^2 - 2*r0*r1 - 3*r1^2)/2
r3 = r0/2 - r1/2 - sqrt((r0 - 3*r1)*(r0 + r1))/2
a = r0 - 2r1 - 2r3
(r0, r1, r2, y2, r3, a)

    (1.8, 0.5000289100968703, 0.2826151986125809, 1.4326296536610128, 0.23471178258110315, 0.3305186146440531)

using Plots

function draw(r0, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    s = (1 + 3sqrt(57))^(1/3)
    r1 = -r0*(2s^2 - s - 16)/3s
    t = sqrt(r0*(r0 + r1))
    u = sqrt(r0^2 - 2r0*r1 - 3r1^2)
    r2 = (r0^2 - r1^2 - t*u)/(2r0 + 2r1)
    y2 = (t + u)/2
    r3 = (r0 - r1 - u)/2
    a = r0 - 2r1 - 2r3
    @printf("a = %g;  r1 = %g;  r3 = %g;  r2 = %g;  y2 = %g\n", a, r1, r3, r2, y2)
    @printf("外円の直径 = %g;  甲円の直径 = %g\n", 2r0, 2r1)
    plot()
    circle(0, 0, r0, :brown)
    circle(0, a + r1, r1)
    circle(2r1, a + r1, r1)
    circle(-2r1, a + r1, r1)
    circle(0, r0 - r3, r3, :blue)
    circle(r1, y2, r2, :green)
    circle(-r1, y2, r2, :green)
    b = sqrt(r0^2 - a^2)
    segment(-b, a, b, a, :magenta)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
        vline!([0], color=:gray80, lw=0.5)
        hline!([0], color=:gray80, lw=0.5)
        point(0, a + r1, " 甲円:r1\n(0,a+r1)", :red, :center, delta=-delta)
        point(2r1, a + r1, " 甲円:r1\n(2r1,a+r1)", :red, :center, delta=-delta)
        point(0, r0 - r3, " 丙円:r3,(0,r0-r3)", :black, :center, :bottom, delta=delta)
        point(r1, y2, " 乙円:r2,(r1,y2)", :black, :left, :vcenter)
        point(0, a, " a", :magenta, :left, :top, delta=-delta/2)
        point(0, r0, " r0", :brown, :left, :bottom, delta=delta/2)
    else
        plot!(showaxis=false)
    end
end;

draw(3.6/2, true)