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

『算法新書』文政13年(1830)

一関市博物館和算に挑戦平成26年度出題問題(3)[上級問題]（高校生以上）
https://www.city.ichinoseki.iwate.jp/museum/wasan/h26/hard.html
キーワード：円5個，外円，弦
#Julia #SymPy #算額 #和算 #数学

外円の中に 2 本の弦を設け，甲円，乙円，丙円，丁円を容れる。甲円，乙円，丙円，の直径が 6 寸，3 寸，4 寸のとき，丁円の直径はいかほどか。

外円の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (x_1,\ r_1)\)
乙円の半径と中心座標を \(r_2,\ (r_2,\ y_2)\)
丙円の半径と中心座標を \(r_3,\ (r_3,\ r_3)\)
丁円の半径と中心座標を \(r_4,\ (r_4,\ r_4)\)
とおき，以下の連立方程式を解く。

算法助術の公式48「外円の弦の上に載っている互いに接し外円に接している 2 円の半径と，円と反対方向にある矢の関係式」を適用する。今の場合，矢は第 3 の円の直径である。
甲円，乙円，丙円のセット，甲円，乙円，丁円のセットに対して，eq1, eq2 を立てる。\(r_1,\ r_2,\ r_3\) が既知で，\(r_4,\ R\) を求める。

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

using SymPy

@syms R::positive, r1::positive, r2::positive, r3::positive, r4::positive
eq1 = (r1 + r2)*2r3 + 2r1*r2 - 2sqrt(2r1*r2*R*2r3)

\( \displaystyle　- 4 \sqrt{R} \sqrt{r_{1}} \sqrt{r_{2}} \sqrt{r_{3}} + 2 r_{1} r_{2} + 2 r_{3} \left(r_{1} + r_{2}\right)\)

eq2 = (r1 + r2)*2r4 + 2r1*r2 - 2sqrt(2r1*r2*R*2r4)

\( \displaystyle　- 4 \sqrt{R} \sqrt{r_{1}} \sqrt{r_{2}} \sqrt{r_{4}} + 2 r_{1} r_{2} + 2 r_{4} \left(r_{1} + r_{2}\right)\)

res = solve([eq1, eq2], (r4, R))[2]; # 2 of 2

# r4 丁円の半径
ans_r4 = res[1] |> factor
ans_r4 |> display
ans_r4 |> println
ans_r4(r1 => 6//2, r2 => 3//2, r3 => 4//2) |> println

\( \displaystyle　\frac{r_{1}^{2} r_{2}^{2}}{r_{3} \left(r_{1} + r_{2}\right)^{2}}\)

r1^2*r2^2/(r3*(r1 + r2)^2)
1/2

甲円，乙円，丙円，の直径が 6 寸，3 寸，4 寸のとき，丁円の直径は 1 寸である。

# R
res[2] |> display
res[2] |> println
res[2](r1 => 6//2, r2 => 3//2, r3 => 4//2) |> println

\( \displaystyle　\frac{\left(r_{1} r_{2} + r_{1} r_{3} + r_{2} r_{3}\right)^{2}}{4 r_{1} r_{2} r_{3}}\)

(r1*r2 + r1*r3 + r2*r3)^2/(4*r1*r2*r3)
81/16

上では，eq1, eq2 から直接 solve() によって \(r_4\) を求めたが，以下の関係式を導くこともできる。

(1/r3)*(1/r4) = (1/r1 + 1/r2)^2

\( \displaystyle　\frac{1}{r_{3} r_{4}} = \left(\frac{1}{r_{2}} + \frac{1}{r_{1}}\right)^{2}\)

以下では，図を描くためのパラメータを求める。

まず，甲円と乙円の中心の \(x\) 座標を求める。
\(y\) 座標は \(2r_3 - R + r_1，2r_3 - R + r_2\) である。

@syms x1::positive, x2::negative

# x1
eq3 = x1^2 + (2r3 - R + r1)^2 - (R - r1)^2
ans_x1 = solve(eq3, x1)[1]
ans_x1 |> println

2*sqrt(r3)*sqrt(R - r1 - r3)

# x2
eq4 = x2^2 + (2r3 - R + r2)^2 - (R - r2)^2
ans_x2 = solve(eq4, x2)[1]
ans_x2 |> println

-2*sqrt(r3)*sqrt(R - r2 - r3)

次に，丁円の中心の座標 (x4, y4) と，弦と円の接点の座標 (x01, sqrt(R^2 - x01^2)), (x02, sqrt(R^2 - x02^2)) を求める。

SymPy では能力的に解析解を求めることができないので，数値解を求める。

@syms r1, r2, r3, r4, R, x1, x2, x4, y4, x01, x02
eq5 = dist(x01, sqrt(R^2 - x01^2), x02, sqrt(R^2 - x02^2), x1, 2r3 - R + r1) - r1^2
eq6 = dist(x01, sqrt(R^2 - x01^2), x02, sqrt(R^2 - x02^2), x2, 2r3 - R + r2) - r2^2
eq7 = dist(x01, sqrt(R^2 - x01^2), x02, sqrt(R^2 - x02^2), x4, y4) - r4^2
# eq8 = x4^2 + y4^2 - (R - r4)^2
eq8 = y4/x4 + ( (x01- x02)/(sqrt(R^2 - x01^2) - sqrt(R^2 - x02^2)));
# res3 = solve([eq5, eq6, eq7, eq8], (x4, y4, x01, x02))

function H(u)
(x4, y4, x01, x02) = u
return [
dist(x01, sqrt(R^2 - x01^2), x02, sqrt(R^2 - x02^2), x1, 2r3 - R + r1) - r1^2,
dist(x01, sqrt(R^2 - x01^2), x02, sqrt(R^2 - x02^2), x2, 2r3 - R + r2) - r2^2,
dist(x01, sqrt(R^2 - x01^2), x02, sqrt(R^2 - x02^2), x4, y4) - r4^2,
# x4^2 + y4^2 - (R - r4)^2
y4/x4 + ( (x01- x02)/(sqrt(R^2 - x01^2) - sqrt(R^2 - x02^2)))
]
end;
(r1, r2, r3, r4, R, x1, x2) = 3.0, 1.5, 2.0, 0.5, 5.0625, 0.7071067811865476, -3.5355339059327378
iniv = BigFloat[-3, 3.5, -0.2, -4.8]
res = nls(H, ini=iniv)

([-2.8677108348121094, 3.548611111111111, -0.20395999275247695, -4.902922315817033], true)

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

function draw(r1, r2, r3, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r4 = r1^2*r2^2/(r3*(r1 + r2)^2)
R = (r1*r2 + r1*r3 + r2*r3)^2/(4*r1*r2*r3)
x1 = 2*sqrt(r3)*sqrt(R - r1 - r3)
x2 = -2*sqrt(r3)*sqrt(R - r2 - r3)
(x4, y4, x01, x02) = [-2.8677108348121094, 3.548611111111111, -0.20395999275247695, -4.902922315817033]
y01 = sqrt(R^2 - x01^2)
y02 = sqrt(R^2 - x02^2)
y = 2r3 - R
x = sqrt(R^2 - y^2)
plot()
circle(0, 0, R, :green)
circle(x1, y + r1, r1)
circle(x2, y + r2, r2, :blue)
circle(0, r3 - R, r3, :magenta)
segment(-x, y, x, y, :brown)
circle(x4, y4, r4, :tomato)
segment(x01, y01, x02, y02, :orange)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
vline!([0], color=:gray80, lw=0.5)
hline!([0], color=:gray80, lw=0.5)
point(0, 0, "外円:R,(0,0)", :green, :center, delta=-delta)
point(x1, y + r1, "甲円:r1,(x1,y+r1)", :red, :center, delta=-delta)
point(x2, y + r2, "乙円:r2,(x2,y+r2)", :blue, :center, delta=-delta)
point(0, r3 - R, "丙円:r3,(0,r3-R)", :magenta, :center, delta=-delta)
point(x4, y4, "丁円:r4,(x4,y4)", :tomato, :right, :bottom, delta=4delta)
point(x01, y01, "(x01,y01)", :green, :center, :bottom, delta=delta)
point(x02, y02, "(x02,y02)", :green, :left, delta=-delta)
point(x, y, "(x,y)", :green, :right, :bottom, delta=delta)
end
end;