愛媛県松山市桜谷町 伊佐爾波神社 文久元年(1861)
愛媛和算研究会：愛媛の和算研究ー現代解法を通してー，明朗社，愛媛県砥部町，平成29年8月1日.
https://ehimewasan.com/wp-content/uploads/2024/12/e5bd925b7d9f489aff074333f26ce7ae.pdf
キーワード：円12個
#Julia #SymPy #算額 #和算 #数学
直径の等しい 3 個の円が交わり，その隙間に天円 1 個，地円 4 個，人円 4 個を容れる。地円の直径が与えられたとき，人円の直径はいかほどか。
外円の半径と中心座標を \(R,\ (0,\ 0),\ (r_1 + R,\ 0)\)
天円の半径と中心座標を \(r_1,\ (0,\ 0)\)
地円の半径と中心座標を \(r_2,\ (r_2,\ y_2)\)
人円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive, r2::positive, y2::positive, r3::positive, x3::positive, y3::positive
R = r1 + 2r2
eq1 = r2^2 + y2^2 - (r1 + r2)^2
eq2 = x3^2 + y3^2 - (r1 + r3)^2
eq3 = (x3 - r2)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq4 = (r1 + R - r2)^2 + y2^2 - (R + r2)^2
eq5 = (r1 + R - x3)^2 + y3^2 - (R + r3)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r3, r1, y2, x3, y3))[1];
# r3: 人円の半径
ans_r3 = res[1]
@show(ans_r3)
ans_r3 = r2*(-2*sqrt(2) - 2 + sqrt(2 + 2*sqrt(2)) + 2*sqrt(1 + sqrt(2)))/2
\(\displaystyle \frac{r_{2} \left(- 2 \sqrt{2} - 2 + \sqrt{2 + 2 \sqrt{2}} + 2 \sqrt{1 + \sqrt{2}}\right)}{2}\)
ans_r3(r2 => 1) |> N
0.23824452512475222
ans_r3（の，\(r_2\) に対する係数）は SymPy では簡約化できない。
算額問題ではよく出てくる式の形式の 1 つである \(\frac{a}{b + \sqrt{c + \sqrt{2}d}}\) を仮定して，以下のような簡単なプログラムで式中の係数が判明する。
長精度計算をしているので，誤差範囲で偶然一致したというようなものではない。
また，この式でうまく行かない場合は，別のタイプの術式で同じように係数を探索すればよい。
s2 = √BigFloat(2)
ans_c = (-2s2 - 2 + sqrt(2 + 2s2) + 2*sqrt(1 + s2))/2
for a = -5:5, b = -5:5, c = -5:5, d = -5:5
c + s2*d < 0 && continue
x = a/(b + sqrt(c + s2*d))
abs( (x - ans_c)/x) < 1e-50 && println( (a, b, c, d))
end
(1, 2, 2, 2)
得られた式は \(\frac{1}{2 + \sqrt{2 + 2\sqrt{2}}}\) であり，術で述べられているものと全く同じものである。
s2 = √BigFloat(2)
1/(2 + sqrt(2 +2*s2))
0.2382445251247522245822714240941272481393047449605175069004029891858178318819107
術 = 1/(2 + sqrt(2 + 2√Sym(2)))
\(\displaystyle \frac{1}{2 + \sqrt{2 + 2 \sqrt{2}}}\)
分母の有理化をすれば，\(ans_{r3}\) と同じ形になる。
@syms d
apart(術, d) |> factor
\(\displaystyle \frac{- 2 \sqrt{2} - 2 + \sqrt{2} \sqrt{1 + \sqrt{2}} + 2 \sqrt{1 + \sqrt{2}}}{2}\)
術(r2 => 1).evalf()
\(0.238244525124752\)
# r1: 天円の半径
ans_r1 = res[1]
@show(ans_r1)
ans_r1 = r2*(-2*sqrt(2) - 2 + sqrt(2 + 2*sqrt(2)) + 2*sqrt(1 + sqrt(2)))/2
\(\displaystyle \frac{r_{2} \left(- 2 \sqrt{2} - 2 + \sqrt{2 + 2 \sqrt{2}} + 2 \sqrt{1 + \sqrt{2}}\right)}{2}\)
# y2
ans_y2 = res[3]
@show(ans_y2)
ans_y2 = r2*sqrt(2 + 2*sqrt(2))
\(r_{2} \sqrt{2 + 2 \sqrt{2}}\)
# x3
ans_x3 = res[4]
@show(ans_x3)
ans_x3 = r2*(-sqrt(1/2 + sqrt(2)/2) + 1 + sqrt(2))
\(\displaystyle r_{2} \left(- \sqrt{\frac{1}{2} + \frac{\sqrt{2}}{2}} + 1 + \sqrt{2}\right)\)
ans_x3 |> factor
@show(ans_x3)
ans_x3 = r2*(-sqrt(1/2 + sqrt(2)/2) + 1 + sqrt(2))
\(\displaystyle r_{2} \left(- \sqrt{\frac{1}{2} + \frac{\sqrt{2}}{2}} + 1 + \sqrt{2}\right)\)
\(y_3 = r_2\) である。
これは，方程式を解いてみないとわからない。
# y3
ans_y3 = res[5]
@show(ans_y3)
ans_y3 = r2
\(r_{2}\)
描画関数プログラムのソースを見る
function draw(r2, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, showaxis=false, label="", fontfamily="IPAMincho")
r3 = r3 = r2*(-2*sqrt(2) - 2 + sqrt(2 + 2*sqrt(2)) + 2*sqrt(1 + sqrt(2)))/2
r1 = sqrt(2)*r2
y2 = r2*sqrt(2 + 2*sqrt(2))
x3 = r2*(-sqrt(1/2 + sqrt(2)/2) + 1 + sqrt(2))
y3 = r2
@printf("天円の直径が %g のとき，人円の直径は %g，地円の直径は %g である。\n", 2r2, 2r3, 2r1)
R = r1 + 2r2
plot()
circle(0, 0, R)
circle2(r1 + R, 0, R)
circle(0, 0, r1, :blue)
circle4(r2, y2, r2, :green)
circle4(x3, y3, r3, :black)
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, R, "R", :red, :center, :bottom, delta=delta)
point(r1 + R, 0, "r1+R", :red, :center, delta=-delta)
point(0, 0, "地円:r1\n(0,0)", :green, :center, delta=-delta)
point(r2, y2, "天円:r2,(r2,y2)", :blue, :left, :bottom, delta=5delta, deltax=-7delta)
point(x3, y3, "人円:r3,(r3,y3)", :black, :left, :vcenter, deltax=3delta)
end
end;
draw(1/2, true)
天円の直径が 1 のとき，人円の直径は 0.238245，地円の直径は 1.41421 である。
以下のアイコンをクリックして応援してください
