算額（その396）

算額（その396）

東京都渋谷区 金王神社（金王八幡宮） 元治元年(1864)
https://gunmawasan.web.fc2.com/files/konnou-HP.pdf

金王八幡宮③「宝物館」（渋谷散歩③）
https://wheatbaku.exblog.jp/30049403/

直線上に大円，中円，小円が載っており，互いに接している。中円，小円の直径がそれぞれ 9 寸，4 寸のとき，大円の直径はいくつか。

大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (0, r2)
小円の半径と中心座標を r3, (x3, r3)
として，以下の連立方程式を立て解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, x1::positive, r2::positive,
     r3::positive, x3::positive;

eq1 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (r1, x1, x3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    ((-2*r2^(11/2)*r3^(3/2) + 8*r2^(9/2)*r3^(5/2) - 12*r2^(7/2)*r3^(7/2) + 8*r2^(5/2)*r3^(9/2) - 2*r2^(3/2)*r3^(11/2) + r2^6*r3 - 3*r2^5*r3^2 + 2*r2^4*r3^3 + 2*r2^3*r3^4 - 3*r2^2*r3^5 + r2*r3^6)/(r2^3 - 3*r2^2*r3 + 3*r2*r3^2 - r3^3)^2, 2*r2^(3/2)*sqrt(r3)/(r2 - r3) - 2*r2*r3/(r2 - r3), 2*sqrt(r2)*sqrt(r3))
    ((2*r2^(11/2)*r3^(3/2) - 8*r2^(9/2)*r3^(5/2) + 12*r2^(7/2)*r3^(7/2) - 8*r2^(5/2)*r3^(9/2) + 2*r2^(3/2)*r3^(11/2) + r2^6*r3 - 3*r2^5*r3^2 + 2*r2^4*r3^3 + 2*r2^3*r3^4 - 3*r2^2*r3^5 + r2*r3^6)/(r2^3 - 3*r2^2*r3 + 3*r2*r3^2 - r3^3)^2, 2*r2^(3/2)*sqrt(r3)/(r2 - r3) + 2*r2*r3/(r2 - r3), 2*sqrt(r2)*sqrt(r3))

2 番目の解が適解である。

   r1 = 18;  x1 = 18;  x3 = 6
   大円の直径 = 36

r1 を求めるために得られた式は長く複雑である。

(r2, r3) = (9, 4) .// 2
(2*r2^(11/2)*r3^(3/2) - 8*r2^(9/2)*r3^(5/2) + 12*r2^(7/2)*r3^(7/2) - 8*r2^(5/2)*r3^(9/2) + 2*r2^(3/2)*r3^(11/2) + r2^6*r3 - 3*r2^5*r3^2 + 2*r2^4*r3^3 + 2*r2^3*r3^4 - 3*r2^2*r3^5 + r2*r3^6)/(r2^3 - 3*r2^2*r3 + 3*r2*r3^2 - r3^3)^2

   18.000000000000014

術では，中円径 ÷ (sqrt(中円径/小円径) - 1)^2 と簡単である。

(中円径, 小円径) = (9, 4)
中円径 / (sqrt(中円径/小円径) - 1)^2

   36.0

この式を得るには，直線上に載っていて外接する 2 円の中心の水平距離に関する公式 d12 = 2sqrt(r1*r2) を用いて以下のように展開する。

r1, r3 から r2 を求める場合も，r1, r2 から r3 を求める場合も同じである。

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (9, 4) .// 2
   (r1, x1, x3) = ((2*r2^(11/2)*r3^(3/2) - 8*r2^(9/2)*r3^(5/2) + 12*r2^(7/2)*r3^(7/2) - 8*r2^(5/2)*r3^(9/2) + 2*r2^(3/2)*r3^(11/2) + r2^6*r3 - 3*r2^5*r3^2 + 2*r2^4*r3^3 + 2*r2^3*r3^4 - 3*r2^2*r3^5 + r2*r3^6)/(r2^3 - 3*r2^2*r3 + 3*r2*r3^2 - r3^3)^2, 2*r2^(3/2)*sqrt(r3)/(r2 - r3) + 2*r2*r3/(r2 - r3), 2*sqrt(r2)*sqrt(r3))
   @printf("r1 = %g;  x1 = %g;  x3 = %g\n", r1, x1, x3)
   @printf("大円の直径 = %g\n", 2r1)
   plot()
   circle(x1, r1, r1)
   circle(0, r2, r2, :blue)
   circle(x3, r3, r3, :orange)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, r1, " 大:r1,(x1,r1)", :red, :center, :bottom, delta=delta)
       point(0, r2, " 中:r2,(0,r2)", :blue, :center, :bottom, delta=delta)
       point(x3, r3, "小:r3,(x3,r3) ", :orange, :right, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;