算額(その1212)
(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)