円に引かれた水平の弦と円弧に囲まれた（狭い方の）隙間に一番大きい円は 1 個ある。二番目以降の円は 2 個ずつある。外円と一番大きい円の直径が与えられたとき，2番目以降の円の直径を求めよ。
「和算の心（その003）」では特別な場合として，弦が円の中心を通る場合について述べた。
今回は，一般の場合について述べる。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive,
r2::positive, x2::positive,
x::positive
y = R - 2r1
eq1 = x2^2 + (y + r2)^2 - (R - r2)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = x/(2R - 2r1) - 2r1/x
eq3 = x^2 + y^2 - R^2
res2 = solve([eq1, eq2, eq3], (r2, x2, x))[1]
(-r1*(-R + r1)/R, 2*r1*sqrt(R - r1)/sqrt(R), 2*sqrt(r1)*sqrt(R - r1))
ans_r2 = res2[1](R=>3, r1=> 1/2)
ans_r2 |> println
ans_x2 = res2[2](R=>3, r1=> 1/2)
ans_x2 |> println
0.416666666666667
0.52704627669473*sqrt(3)
@syms r3::positive, x3::positive
eq4 = x3^2 + (y + r3)^2 - (R - r3)^2
eq5 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2
res3 = solve([eq4, eq5], (r3, x3))[1]
res3 |> println
res3[1](R => 3, r1 => 1/2, r2 => ans_r2, x2 => ans_x2).evalf() |> println
res3[2](R => 3, r1 => 1/2, r2 => ans_r2, x2 => ans_x2).evalf() |> println
(-(4*R*r1 - 4*r1^2 + x2^2 - 2*x2*(-sqrt(r2)*sqrt( (-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2)))/(4*(-R + r1 - r2)), -sqrt(r2)*sqrt( (-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2))
0.255102040816327
1.56492159287190
@syms r4::positive, x4::positive
eq6 = x4^2 + (y + r4)^2 - (R - r4)^2
eq7 = (x4 - x3)^2 + (r3 - r4)^2 - (r3 + r4)^2
res4 = solve([eq6, eq7], (r4, x4))[1]
(-(4*R*r1 - 4*r1^2 + x3^2 - 2*x3*(-sqrt(r3)*sqrt( (-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r3 + x3^2))/(-R + r1 - r3) + x3*(-R + r1)/(-R + r1 - r3)))/(4*(-R + r1 - r3)), -sqrt(r3)*sqrt( (-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r3 + x3^2))/(-R + r1 - r3) + x3*(-R + r1)/(-R + r1 - r3))
\(r_3, r_4\) の式を比較すれば，変数名の置き換えだけで次々と半径，と中心座標(\(x\) 座標)が求まることがわかる。
そこで，次のような関数を定義しておけば，簡単に計算を続けることができる。
\(R = 2, r_1 = 1/2\) のとき，円の半径は以下のようになる（分数は近似値）。
1: 0.5 1//2
2: 0.375 3//8
3: 0.18 9//50
4: 0.0688775510204081 27//392
5: 0.0240928019036288 79//3279
6: 0.00816312819134647 141//17273
7: 0.00273597297804469 2//731
8: 0.000913659014267574 2//2189
9: 0.00030473867949933 1//3281
10: 0.000101600202938299 1//9833
この場合もまた「和算の心（その003）」と同じく，簡単な一般項は求めがたい。
\(r_1 = R/2\) のときも，この式が使える。
nextcircle(r2, x2, R, r1) = (-(4*R*r1 - 4*r1^2 + x2^2 - 2*x2*(-sqrt(r2)*sqrt( (-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2)))/(4*(-R + r1 - r2)), -sqrt(r2)*sqrt( (-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2))
描画関数プログラムのソースを見る
function draw(r1, R, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, x2, x) = (-r1*(-R + r1)/R, 2*r1*sqrt(R - r1)/sqrt(R), 2*sqrt(r1)*sqrt(R - r1))
y = R - 2r1
plot()
circlef(0, 0, R, :aliceblue)
circle(0, 0, R, :black)
circlef(0, R - r1, r1, 1)
@printf("%2d: %.15g ", 1, r1)
println(rationalize(r1, tol=1e-7))
for i = 2:10
@printf("%2d: %.15g ", i, r2)
println(rationalize(r2, tol=1e-7))
circle2f(x2, y + r2, r2, i)
(r2, x2) = nextcircle(r2, x2, R, r1)
end
segment(-x, y, x, y)
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)
end
end;
以下のアイコンをクリックして応援してください
