算額（その2012）

算額（その2012）

(20) 兵庫県青垣町遠坂字後岶 熊野神社 明治18年(1885)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円8個

互いに外接し合う，木円 2 個，火円 2 個，土円 2 個，金円 1 個，水円 1 個がある。
木円，火円の直径が与えられたとき，金円の直径を求める術を述べよ。

木円の半径と中心座標を r1, (r1, y1)
火円の半径と中心座標を r2, (r2, y2)
土円の半径と中心座標を r3, (r3, 0)
金円の半径と中心座標を r4, (0, y4)
水円の半径と中心座標を r5, (0, y5)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, y1::negative, r2::positive, y2::positive, r3::positive, r4::positive, y4::negative, r5::positive, y5::positive

eq1 = (r1 - r2)^2 + (y2 - y1)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + y1^2 - (r1 + r3)^2
eq3 = (r2 - r3)^2 + y2^2 - (r2 + r3)^2;

res1 = solve([eq1, eq2, eq3], (y1, y2, r3))[1]

   (-2*sqrt(r1)*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))/sqrt(r2), 2*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2)), (-2*r1^(3/2)*r2^(3/2)*(r1 - r2)^2 + r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2))/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)))

eq4 = r1^2 + (y1 - y4)^2 - (r1 + r4)^2
eq5 = r3^2 + y4^2 - (r3 + r4)^2
eq6 = r2^2 + (y2 - y5)^2 - (r2 + r5)^2
eq7 = r3^2 + y5^2 - (r3 + r5)^2;
#res2 = solve([eq4, eq5, eq6, eq7], (r4, y4, r5, y5))

次いで，r4, y4, r5, y5 を求めるが，前段で求めた y1, y2, r3 は数式が複雑なので代入せず，変数のまま解く。

res2 = solve([eq4, eq5, eq6, eq7], (r4, y4, r5, y5))[4]

   (y1*(y1*(r1 - r3)*sqrt(4*r1*r3 + y1^2)/(-r1^2 + 2*r1*r3 - r3^2 + y1^2) + y1 - y1*(2*r1*r3 - 2*r3^2 + y1^2)/((-r1 + r3 + y1)*(r1 - r3 + y1)))/(2*(r1 - r3)), -y1*(r1 - r3)*sqrt(4*r1*r3 + y1^2)/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2)) + y1*(2*r1*r3 - 2*r3^2 + y1^2)/(2*(-r1 + r3 + y1)*(r1 - r3 + y1)), y2*(y2*(r2 - r3)*sqrt(4*r2*r3 + y2^2)/(-r2^2 + 2*r2*r3 - r3^2 + y2^2) + y2 - y2*(2*r2*r3 - 2*r3^2 + y2^2)/((-r2 + r3 + y2)*(r2 - r3 + y2)))/(2*(r2 - r3)), -y2*(r2 - r3)*sqrt(4*r2*r3 + y2^2)/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2)) + y2*(2*r2*r3 - 2*r3^2 + y2^2)/(2*(-r2 + r3 + y2)*(r2 - r3 + y2)))

res2[1] |> simplify |> println

   y1^2*(-r1 - r3 + sqrt(4*r1*r3 + y1^2))/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))

res2[2] |> simplify |> println

   y1*(2*r1*r3 - r1*sqrt(4*r1*r3 + y1^2) - 2*r3^2 + r3*sqrt(4*r1*r3 + y1^2) + y1^2)/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))

res2[3] |> simplify |> println

   y2^2*(-r2 - r3 + sqrt(4*r2*r3 + y2^2))/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2))

res2[4] |> simplify |> println

   y2*(2*r2*r3 - r2*sqrt(4*r2*r3 + y2^2) - 2*r3^2 + r3*sqrt(4*r2*r3 + y2^2) + y2^2)/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2))

金円の直径の式に y1, y2, r3 を代入すればよいが，SymPy では簡約化できないとてつもなく長い式になる。

y1 = -2*sqrt(r1)*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))/sqrt(r2)
y2 = 2*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))
r3 = (-2*r1^(3/2)*r2^(3/2)*(r1 - r2)^2 + r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2))/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2))
金径 = y1^2*(-r1 - r3 + sqrt(4*r1*r3 + y1^2))/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))
金径 |> println

   4*r1*(-r1 + sqrt(4*r1*(r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2) - 2*r1^1.5*r2^1.5*(r1 - r2)^2)/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)) + 4*r1*((r1^2*r2^2 + r1*r2^3 - 2*r1^1.5*r2^2.5)/(r1^2 - 2*r1*r2 + r2^2))/r2) - (r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2) - 2*r1^1.5*r2^1.5*(r1 - r2)^2)/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)))*(r1^2*r2^2 + r1*r2^3 - 2*r1^1.5*r2^2.5)/(r2*(r1^2 - 2*r1*r2 + r2^2)*(-2*r1^2 + 4*r1*(r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2) - 2*r1^1.5*r2^1.5*(r1 - r2)^2)/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)) + 8*r1*((r1^2*r2^2 + r1*r2^3 - 2*r1^1.5*r2^2.5)/(r1^2 - 2*r1*r2 + r2^2))/r2 - 2*(r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2) - 2*r1^1.5*r2^1.5*(r1 - r2)^2)^2/((r1 - r2)^4*(r1^2 - 2*r1*r2 + r2^2)^2)))

r1，r2 に具体例を代入すれば，正しい金円の直径が求まるので，式は間違いない。

金径(r1 => 8/2, r2 => 10/2) |> println

   0.804254212447020

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

   r1 = 4;  y1 = 4.22291;  r2 = 5;  y2 = -4.72136;  r3 = 1.11456;  r4 = 0.804254;  y4 = -1.64362;  r5 = 0.871325;  y5 = -1.64362

function draw(r1, r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   y1 = -2*sqrt(r1)*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))/sqrt(r2)
   y2 = 2*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))
   r3 = (-2*r1^(3/2)*r2^(3/2)*(r1 - r2)^2 + r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2))/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2))
   (y1, y2, y3) =  (-2*sqrt(r1)*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))/sqrt(r2), 2*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2)), (-2*r1^(3/2)*r2^(3/2)*(r1 - r2)^2 + r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2))/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)))
   r4 = y1^2*(-r1 - r3 + sqrt(4*r1*r3 + y1^2))/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))
   y4 = y1*(2*r1*r3 - r1*sqrt(4*r1*r3 + y1^2) - 2*r3^2 + r3*sqrt(4*r1*r3 + y1^2) + y1^2)/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))
   r5 = y2^2*(-r2 - r3 + sqrt(4*r2*r3 + y2^2))/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2))
   y5 = y2*(2*r2*r3 - r2*sqrt(4*r2*r3 + y2^2) - 2*r3^2 + r3*sqrt(4*r2*r3 + y2^2) + y2^2)/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2))
   @printf("r1 = %g;  y1 = %g;  r2 = %g;  y2 = %g;  r3 = %g;  r4 = %g;  y4 = %g;  r5 = %g;  y5 = %g\n", r1, y1, r2, y2, r3, r4, y5, r5, y5)
   plot()
   circle2(r1, y1, r1)
   circle2(r2, y2, r2, :blue)
   circle2(r3, 0, r3, :green)
   circle(0, y4, r4, :magenta)
   circle(0, y5, r5, :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(r1, y1, "木円:r1,(r1,y1)", :red, :center, delta=-delta/2)
       point(r2, y2, "火円:r2,(r2,y2)", :blue, :center, delta=-delta/2)
       point(r3, 0, "土円:r3,(r3,0)", :green, :left, delta=-delta/2)
       point(0, y4, " 金円:r4,(0,r4)", :magenta, :left, :vcenter)
       point(0, y5, " 水円:r5,(0,r5)", :orange, :left, :vcenter)
   end
end;

draw(10/2, 8/2, true)