算額（その1708） - 算額あれこれ

山形県東置賜郡川西町八幡神社明治18年(1885)

道脇義正，浜田敏夫，大山真：算額に表れた類題の研究（3.承前）fig.40 ，数学史研究，通巻60号，p.1-26，1974年1月〜3月.
http://www.wasan.jp/sugakusipdf/sugakusi060.pdf
キーワード：円5個，台形
#Julia #SymPy #算額 #和算 #数学

等脚台形の中に甲円，乙円，丙円，丁円，戊円を容れる。丙円，丁円の直径が与えられたとき，甲円の直径はいかほどか。

横倒しになっているが，等脚台形の下底，上底，高さを \(2a, 2b, h\)
甲円の半径と中心座標を \(r_1, (x_1, r_1)\)
乙円の半径と中心座標を \(r_2, (h - r_2, y_2)\)
丙円の半径と中心座標を \(r_3, (h - r_3, r_3)\)
丁円の半径と中心座標を \(r_4, (r_4, r_4)\)
戊円の半径と中心座標を \(r_5, (r_5, y_5)\)
とおく。
それぞれのパラメータは，連立方程式を解けば求められる。

しかし，甲円の直径を知るだけであれば簡単に求めることができる。
そのために下図を描く。

まず甲円，丙円，丁円の部分を切り取った図を描く。次に上下幅が \(2r_1\) の二本の平行線を描き，丁円の左に接する縦線を描く。丙円のコピー（丙円' と表示した）を左上隅に置く。このようにして描いた丙円'は丁円と接する。図に示す直角三角形の三辺 A, B, C において，B は丙円' と丁円の中心の垂直方向の距離なので B \(= \sqrt{r_3 r_4}\) である。

\(2r_1 = r_4 + \text{B} + r_3\) であることがわかる。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy

@syms r1::positive, r3::positive, r4::positive, x1::positive
eq1 = (x1 - r4)^2 + (r1 - r4)^2 - (r1 + r4)^2
eq2 = (r3 - r4)^2 + (2r1 - r3 - r4)^2 - (r3 + r4)^2;

res = solve([eq1, eq2], (r1, x1))[4]

( (4*sqrt(r3)*r4^(3/2) + 2*r3*r4 + 2*r4^2)/(4*r4), r4 + sqrt(2)*sqrt(2*sqrt(r3)*r4^(3/2) + r3*r4 + r4^2))

# r1 甲円の半径
res[1] |> simplify

\(\displaystyle \sqrt{r_{3}} \sqrt{r_{4}} + \frac{r_{3}}{2} + \frac{r_{4}}{2}\)

res[1](r3 => 0.5/2, r4 => 0.4/2) * 2 |> println

0.897213595499958

術は「甲円の直径 = sqrt(丙円の直径*丁円の直径) + (丙円の直径 + 丁円の直径)/2」で，等価なものである。

sqrt(0.5*0.4) + (0.5 + 0.4)/2

0.8972135954999579

すべてのパラメータを求めるためには以下の連立方程式を解く。
SymPy の性能では解析解を求めることができないので，数値解を求めることにする。

@syms a, b, h, r1, x1, r2, y2, r3, r4, r5, y5
x1 = r4 + 2sqrt(r1*r4)
eq1 = (h - r2 - x1)^2 + (y2 - r1)^2 - (r1 + r2)^2
eq2 = (h - r3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x1 - r4)^2 + (r1 - r4)^2 - (r1 + r4)^2
eq4 = (x1 - r5)^2 + (y5 - r1)^2 - (r1 + r5)^2
eq5 = (r2 - r3)^2 + (y2 - r3)^2 - (r2 + r3)^2
eq6 = (r5 - r4)^2 + (y5 - r4)^2 - (r4 + r5)^2
eq7 = h - (x1 + 2sqrt(r1*r3) + r3)
eq8 = dist2(0, a, h, b, h - r2, y2, r2)/h
eq9 = dist2(0, a, h, b, r5, y5, r5)/h;
eq10 = dist2(0, a, h, b, x1, r1, r1) |> numerator;
# res = solve([eq1, eq2, eq4, eq5, eq6, eq7, eq8, eq9], (r1, h, r2, y2, r5, y5, a, b))

println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
println(eq6, ", # eq6")
println(eq7, ", # eq7")
println(eq8, ", # eq8")
println(eq9, ", # eq9")
println(eq10, ", # eq10")

(-r1 + y2)^2 - (r1 + r2)^2 + (h - r2 - r4 - 2*sqrt(r1*r4))^2, # eq1
(r1 - r3)^2 - (r1 + r3)^2 + (h - r3 - r4 - 2*sqrt(r1*r4))^2, # eq2
4*r1*r4 + (r1 - r4)^2 - (r1 + r4)^2, # eq3
(-r1 + y5)^2 - (r1 + r5)^2 + (r4 - r5 + 2*sqrt(r1*r4))^2, # eq4
(r2 - r3)^2 - (r2 + r3)^2 + (-r3 + y2)^2, # eq5
(-r4 + r5)^2 + (-r4 + y5)^2 - (r4 + r5)^2, # eq6
h - r3 - r4 - 2*sqrt(r1*r3) - 2*sqrt(r1*r4), # eq7
2*a*b*r2 - 2*a*r2*y2 + b^2*h - 2*b^2*r2 - 2*b*h*y2 + 2*b*r2*y2 - h*r2^2 + h*y2^2, # eq8
a^2*h - 2*a^2*r5 + 2*a*b*r5 - 2*a*h*y5 + 2*a*r5*y5 - 2*b*r5*y5 - h*r5^2 + h*y5^2, # eq9
a^2*h^2 - 2*a^2*h*r4 - 4*a^2*h*sqrt(r1*r4) - a^2*r1^2 + 4*a^2*r1*r4 + a^2*r4^2 + 4*a^2*r4*sqrt(r1*r4) + 2*a*b*h*r4 + 4*a*b*h*sqrt(r1*r4) + 2*a*b*r1^2 - 8*a*b*r1*r4 - 2*a*b*r4^2 - 8*a*b*r4*sqrt(r1*r4) - 2*a*h^2*r1 + 2*a*h*r1*r4 + 4*a*h*r1*sqrt(r1*r4) - b^2*r1^2 + 4*b^2*r1*r4 + b^2*r4^2 + 4*b^2*r4*sqrt(r1*r4) - 2*b*h*r1*r4 - 4*b*h*r1*sqrt(r1*r4), # eq10

function driver(r3, r4)
function H(u)
(r1, h, r2, y2, r5, y5, a, b) = u
return [
(-r1 + y2)^2 - (r1 + r2)^2 + (h - r2 - r4 - 2*sqrt(r1*r4))^2, # eq1
#(r1 - r3)^2 - (r1 + r3)^2 + (h - r3 - r4 - 2*sqrt(r1*r4))^2, # eq2
#4*r1*r4 + (r1 - r4)^2 - (r1 + r4)^2, # eq3
(-r1 + y5)^2 - (r1 + r5)^2 + (r4 - r5 + 2*sqrt(r1*r4))^2, # eq4
(r2 - r3)^2 - (r2 + r3)^2 + (-r3 + y2)^2, # eq5
(-r4 + r5)^2 + (-r4 + y5)^2 - (r4 + r5)^2, # eq6
h - r3 - r4 - 2*sqrt(r1*r3) - 2*sqrt(r1*r4), # eq7
h*(2*a*b*r2 - 2*a*r2*y2 + b^2*h - 2*b^2*r2 - 2*b*h*y2 + 2*b*r2*y2 - h*r2^2 + h*y2^2), # eq8
h*(a^2*h - 2*a^2*r5 + 2*a*b*r5 - 2*a*h*y5 + 2*a*r5*y5 - 2*b*r5*y5 - h*r5^2 + h*y5^2), # eq9
(a^2*h^2 - 2*a^2*h*r4 - 4*a^2*h*sqrt(r1*r4) - a^2*r1^2 + 4*a^2*r1*r4 + a^2*r4^2 + 4*a^2*r4*sqrt(r1*r4) + 2*a*b*h*r4 + 4*a*b*h*sqrt(r1*r4) + 2*a*b*r1^2 - 8*a*b*r1*r4 - 2*a*b*r4^2 - 8*a*b*r4*sqrt(r1*r4) - 2*a*h^2*r1 + 2*a*h*r1*r4 + 4*a*h*r1*sqrt(r1*r4) - b^2*r1^2 + 4*b^2*r1*r4 + b^2*r4^2 + 4*b^2*r4*sqrt(r1*r4) - 2*b*h*r1*r4 - 4*b*h*r1*sqrt(r1*r4))/(a^2 - 2*a*b + b^2 + h^2), # eq10
]
end;

iniv = BigFloat[0.56509, 2.17787, 0.48455, 1.22943, 0.15855, 0.55803, 0.63916, 2.14978]
res = nls(H, ini=iniv)
return(res)
end;

driver(0.5/2, 0.4/2)

([0.44860679774997897, 1.7188516351015686, 0.2741870795471826, 0.7736287611917269, 0.1823572433922654, 0.5819499898073206, 0.7276313970915846, 1.1168434577825994], true)

丙円の直径が 0.5，丁円の直径が 0.4 のとき，甲円の直径は 2*0.44860679774997897 = 0.8972135954999579 である。

2*0.44860679774997897

0.8972135954999579

術は「甲円の直径 = sqrt(丙円の直径*丁円の直径) + (丙円の直径 + 丁円の直径)/2」となっており，これに従うと sqrt(0.5*0.4) + (0.5 + 0.4)/2 = 0.8972135954999579 で，全桁一致する。

描画関数プログラムのソースを見る

function draw(r3, r4, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
answer = driver(r3, r4)
if answer[2]
println(answer)
(r1, h, r2, y2, r5, y5, a, b) = answer[1]
else
println("収束しませんでした")
return
end
x1 = r4 + 2sqrt(r1*r4)
r10 = sqrt(r2*r3) + (r2 + r3)/2
println("r10 = ", r10)
plot([0, h, h, 0, 0], [0, 0, b, a, 0], color=:green, lw=0.5)
circle(x1, r1, r1)
circle(h - r2, y2, r2, :blue)
circle(h - r3, r3, r3, :magenta)
circle(r4, r4, r4, :brown)
circle(r5, 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(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
point(h - r2, y2, "乙円:r2,(h-r2,y2)", :blue, :center, delta=-delta)
point(h - r3, r3, "丙円:r3,(h-r3,r3)", :magenta, :center, delta=-delta)
point(r4, r4, "丁円:r4\n(r4,r4)", :brown, :center, delta=-delta)
point(r5, y5, "戊円:r5\n(r5,y5)", :orange, :center, delta=-delta)
point(h, b, "(h,b)", :green, :right, :bottom, delta=delta/2)
point(0, a, "a", :green, :left, :bottom, delta=delta, deltax=delta/2)
end
end;