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

56 岩手県花泉町金沢字大門沢大門神社慶応2年(1866)

安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円1個，半円,正方形
#Julia #SymPy #算額 #和算 #数学

半円の中に甲，乙の半円と丙円，直角三角形，正方形を容れる。正方形の一辺の長さが 15 寸，半円の直径が 64 寸のとき，丙円の直径はいかほどか。

半円の半径と中心座標を \(R,\ (0,\ 0)\)
甲半円の半径と中心座標を \(r_1,\ (x_1,\ y_1)\)
乙半円の半径と中心座標を \(r_2,\ (x_2,\ y_2)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
\(x\) 軸上の半円の端点と正方形の頂点の座標を \( (x_{01},\ 0),\ (x_{02},\ 0),\ (x_{03},\ 0)\)
甲半円と乙半円の直径の交点座標を \( (x_0,\ y_0)\)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt"); # julia-source.txt ソース
function driver(S, R)
function H(u)
function parameters()
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = x3^2 + y3^2 - (R - r3)^2
eq4 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq6 = (x01 - x03)^2 - (4r1^2 + 4r2^2)
eq7 = (x01 - x0)^2 + y0^2 - 4r1^2
eq8 = (x0 - x03)^2 + y0^2 - 4r2^2
eq9 = S/(x02 - x03) - 2r1/(x01 - x03)
eq10 = S/(x01 - x02) - 2r2/(x01 - x03)
eq11 = y1/(x01 - x1) - y0/(x01 - x0)
eq12 = y2/(x2 - x03) - y0/(x0 - x03)
eq13 = (x1 - x0)^2 + (y0 - y1)^2 - r1^2
eq14 = (x0 - x2)^2 + (y0 - y2)^2 - r2^2
return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8,
eq9, eq10, eq11, eq12, eq13, eq14]
end;
(r1, x1, y1, r2, x2, y2, r3, x3, y3, x01, x02, x03, x0, y0) = u
return parameters()
end;
iniv = BigFloat[
21.68, 4.27, 9.34,
10.34, -19.75, 9.34,
3.55, -16.96, 22.89,
23.84, -8.95, -24.21,
-15.30, 18.67]
res = nls(H, ini=iniv)
res[2] || println("収束していない")
return res[1]
end;
(r1, x1, y1, r2, x2, y2, r3, x3, y3, x01, x02, x03, x0, y0) = driver(15, 64/2)
(r1, x1, y1, r2, x2, y2, r3, x3, y3, x01, x02, x03, x0, y0)

(20.0, 6.173949065130318, 10.289915108550531, 12.0, -17.149858514250884, 10.289915108550531, 3.8269453137801985, -12.288540862582682, 25.351780486217674, 23.323807579381203, -5.830951894845301, -23.323807579381203, -10.975909449120566, 20.579830217101062)

正方形の一辺の長さが 15 寸，半円の直径が 64 寸のとき，丙円の直径は 7.653890627560397 寸である。

# 丙円の直径
2r3

7.653890627560397

術は，置八個開平方以減三個余乗方面以半円径除之以減二個以除方面
計算式は，
\(f(方面, 半円径) = 方面/(2 - 方面(3 - √8)/半円径)\)
となり，
\(f(15, 64) = 7.653890627560397\)
で，答と一致した。

f(方面, 半円径) = 方面/(2 - 方面*(3 - √8)/半円径);
f(15, 64)

7.653890627560397

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

function draw(S, R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, showaxis=false, label="", fontfamily="IPAexMincho")
(r1, x1, y1, r2, x2, y2, r3, x3, y3, x01, x02, x03, x0, y0) = driver(S, R)
@printf("丙円の直径 = %.15g\n", 2r3)
plot([x01, x0, x03], [0, y0, 0], color=:green, lw=0.5)
circle(0, 0, R, :magenta, beginangle=0, endangle=180)
θ1 = atand(y0, x01 - x0)
circle(x1, y1, r1, beginangle=-θ1, endangle=180 - θ1)
θ2 = atand(y0, x0 - x03)
circle(x2, y2, r2, :blue, beginangle=θ2, endangle=180 + θ2)
circle(x3, y3, r3, :orange)
segment(x02, 0, x02 - S*cosd(θ1), S*sind(θ1))
segment(x02, 0, x02 + S*cosd(θ2), S*sind(θ2))
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(x01, 0, "(x01,0)", :green, :center, delta=-delta)
point(x02, 0, "(x02,0)", :green, :center, delta=-delta)
point(x03, 0, "(x03,0)", :green, :center, delta=-delta)
point(0, 0, "(0,0)", :green, :left, delta=-delta, deltax=-4delta)
point(R, 0, "(R,0)", :green, :right, delta=-delta, deltax=4delta)
point(x0, y0, "(x0,y0)", :green, :left, :vcenter, deltax=2delta)
point(x1, y1, "甲半円:r1,(x1,y1)", :red, :left, :bottom, delta=delta)
point(x2, y2, "乙半円:r2\n(x2,y2)", :blue, :right, :bottom, delta=delta)
point(x3, y3, "丙円:r3,(x3,y3)", :orange, :right, :bottom, delta=7delta, deltax=-4delta)
end
end;
draw(15, 64/2, true)