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

12 岩手県胆沢町若柳市野々個人宅安政2年(1855)

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

正方形の中に甲円を 2 個と，大小（長短）の斜線を引き，区画された領域に乙円を容れる。乙円の直径が 1 寸のとき，短い方の斜線の長さはいかほどか。

正方形の一辺の長さを \(a\)
長い方の斜線と正方形の辺との交点を \( (x_{02}, a)\)
短い方の斜線の甲円の円周上の端点を \( (x_{01}, y_{01})\)
甲円の半径と中心座標を \(r_1, (r_1, r_1), (r_1, a - r_1)\)
乙円の半径と中心座標を \(r_2, (a - r_2, y_2)\)
とおき，以下の連立方程式を解く。
短い方の斜線の長さは \(\sqrt{ (a - x_{01})^2 + (a - y_{01})^2}\) である。

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

using SymPy
@syms a::positive, r1::positive, r2::positive,
y2::positive, x01::positive,
y01::positive, x02::positive
r1 = a/4
eq1 = dist2(x01, y01, a, a, r1, a - r1, r1)
eq2 = dist2(x01, y01, a, a, a - r2, y2, r2)
eq3 = dist2(x02, a, a, 0, r1, a - r1, r1)
eq4 = dist2(x02, a, a, 0, a - r2, y2, r2)
eq5 = (x01 - r1)^2 + (y01 - r1)^2 - r1^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, y2, x01, y01, x02))

function H(u)
(a, y2, x01, y01, x02) = u
return [
a^2*(a^2 + 3*a*x01 - 5*a*y01 - 3*x01*y01 + 4*y01^2), # eq1
a^4 - 2*a^3*r2 - 2*a^3*x01 - 2*a^3*y2 - a^2*r2^2 + 2*a^2*r2*x01 + 2*a^2*r2*y01 + 2*a^2*r2*y2 + a^2*x01^2 + 4*a^2*x01*y2 + a^2*y2^2 + 2*a*r2^2*x01 - 2*a*r2*x01*y01 - 2*a*r2*x01*y2 - 2*a*r2*y01*y2 - 2*a*x01^2*y2 - 2*a*x01*y2^2 - r2^2*x01^2 + 2*r2*x01*y01*y2 + x01^2*y2^2, # eq2
a^2*(-a^2 + a*x02 + 4*x02^2), # eq3
-a^2*r2^2 - 2*a^2*r2*y2 + a^2*y2^2 + 2*a*r2^2*x02 + 2*a*r2*x02*y2 - 2*a*x02*y2^2 - r2^2*x02^2 + x02^2*y2^2, # eq4
-a^2/16 + (-a/4 + x01)^2 + (-a/4 + y01)^2, # eq5
]
end;

r2 = 1/2
iniv = BigFloat[10, 6.40388, 3.2, 4.9, 3.90388]
res = nls(H, ini=iniv)

([2.780776406404415, 1.7807764064044151, 0.8898484500494128, 1.3625804391381635, 1.0855823048033113], true)

乙円の直径が 1 のとき，小斜の長さは 2.36365994544375 という結果になった。

術は「乙円の直径を \( (\sqrt{49.13} + 11.9)/8\) 倍すれば小斜」ということで，\( (\sqrt{49.13} + 11.9)/8 = 2.363659945443753 \)となり，一致した。

全てのパラメータは以下のとおりである。
r2 = 0.5; a = 2.78078; x01 = 0.889848; y01 = 1.36258; x02 = 1.08558

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

function draw(r2, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, y2, x01, y01, x02) = res[1]
r1 = a/4
小斜 = sqrt( (a - x01)^2 + (a - y01)^2)
@printf("乙円の直径が %g のとき，小斜の長さは %.15g である。\n", 2r2, 小斜)
@printf("r2 = %g; a = %g; x01 = %g; y01 = %g; x02 = %g\n", r2, a, x01, y01, x02)
plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
circle(r1, r1, r1)
circle(r1, a - r1, r1)
circle(a - r2, y2, r2, :blue)
segment(x01, y01, a, a, :orange, lw=1)
segment(x02, a, a, 0, :magenta)
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(0, a, "a ", :green, :right, :vcenter)
point(a, 0, " a", :green, :center, delta=-delta/2)
point(a, a, "(a,a)", :orange, :center, :bottom, delta=delta/2)
point(x02, a, "(x02,a)", :magenta, :center, :bottom, delta=delta/2)
point(x01, y01, "(x01,y01)", :orange, :right, delta=-delta/2)
point(r1, r1, "甲円:r1,(r1,r1)", :red, :center, delta=-delta/2)
point(r1, a - r1, "甲円:r1,(r1,a-r1)", :red, :center, delta=-delta/2)
point(a - r2, y2, "乙円:r2,(a-r2,y2)", :blue, :center, delta=-delta/2)
plot!(xlims=(-6delta, a+3delta), ylims=(-4delta, a+3delta))
end
end;