円に引かれた水平の弦と円弧に囲まれた（狭い方の）隙間に一番大きい円は 1 個ある。二番目以降の円は 2 個ずつある。一番大きい円の直径を 1 としたとき，2番目以降の円の直径を求めよ。
この問題を解く場合に，算法助法の公式29が役に立つ。
外円の中心を原点とする。
外円の半径を R とすれば，一番大きい円の半径 r1 は r1 = R/2 である。
弦の長さを length とすれば，length^2 = 4(2R)*(2r2) = 16R*r2 であるというのが公式29である。
ここでは特別な場合として，弦が円の中心を通る場合について述べる。この場合，length = 2R である。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R, r1, r2, length
formula29 = length^2 - 4(2R)*(2r2)
formula29 = formula29(length => 2R) |> simplify
formula29 |> println
4*R*(R - 4*r2)
R ≠ 0 なので，R - 4r2 = 0 ということである。r2 について解けば， r2 = R/4 すなわち，2 番目に大きい円の半径は，外円の半径の 1/4 である。
r1 の半径は R/2 だったので，更にその 1/2 である。
公式29を直接使わなくても，以下の連立方程式を解けば，同じ結果が得られる。
このやり方だと，2番目に大きい円の中心座標（x 座標）も同時に得ることができる。また，3番目以降の円の半径を求めることもできる。
using SymPy
@syms R::positive, r1::positive,
r2::positive, x2::positive
R = 2r1
eq1 = x2^2 + r2^2 - (R - r2)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
res2 = solve([eq1, eq2], (r2, x2))[1]
(r1/2, sqrt(2)*r1)
3 番目に大きい円の半径も同様に求めることができる。
@syms r3::positive, x3::positive
r2 = res2[1]
x2 = res2[2]
eq3 = x3^2 + r3^2 - (R - r3)^2
eq4 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2
res3 = solve([eq3, eq4], (r3, x3))[1]
(r1/9, 4*sqrt(2)*r1/3)
4 番目以降の円の半径は前式を初期値として，漸化式で表すことができる。
@syms r4::positive, x4::positive
r3 = res3[1]
x3 = res3[2]
eq5 = x4^2 + r4^2 - (R - r4)^2
eq6 = (x4 - x3)^2 + (r3 - r4)^2 - (r3 + r4)^2
res4 = solve([eq5, eq6], (r4, x4))[1]
(r1/50, 7*sqrt(2)*r1/5)
@syms r5::positive, x5::positive
r4 = res4[1]
x4 = res4[2]
eq7 = x5^2 + r5^2 - (R - r5)^2
eq8 = (x5 - x4)^2 + (r4 - r5)^2 - (r4 + r5)^2
res5 = solve([eq7, eq8], (r5, x5))[1]
(r1/289, 24*sqrt(2)*r1/17)
@syms r6::positive, x6::positive
r5 = res5[1]
x5 = res5[2]
eq9 = x6^2 + r6^2 - (R - r6)^2
eq10 = (x6 - x5)^2 + (r5 - r6)^2 - (r5 + r6)^2
res6 = solve([eq9, eq10], (r6, x6))[1]
(r1/1682, 41*sqrt(2)*r1/29)
一番大きい円の半径から順に列挙すると，r1 = R/2 で，
r1, r1/2, r1/9, r1/50, r1/289, r1/1682, ...
となる。一見何らかの規則性がありそうに見えるが，一筋縄ではいかないようだ。いずれにしろ，指数曲線よりは早く減衰するようだ。そのせいかどうか，この場合の累円を求めよという算額はまだ見たことがない。
描画関数プログラムのソースを見る
function draw(r1, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, x2) = (r1/2, sqrt(2)*r1)
(r3, x3) = (r1/9, 4*sqrt(2)*r1/3)
(r4, x4) = (r1/50, 7*sqrt(2)*r1/5)
(r5, x5) = (r1/289, 24*sqrt(2)*r1/17)
(r6, x6) = (r1/1682, 41*sqrt(2)*r1/29)
R = 2r1
@printf("""r1 = %.15g;
r2 = %.15g; x2 = %.15g;
r3 = %.15g; x3 = %.15g;
r4 = %.15g; x4 = %.15g;
r5 = %.15g; x5 = %.15g;
r6 = %.15g; x6 = %.15g;\n""",
r1, r2, x2, r3, x3, r4, x4, r5, x5, r6, x6)
plot()
circle(0, 0, R)
circle22f(0, r1, r1, 1)
circle4f(x2, r2, r2, 2)
circle4f(x3, r3, r3, 3)
circle4f(x4, r4, r4, 4)
circle4f(x5, r5, r5, 5)
segment(-R, 0, R, 0, :green)
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;
以下のアイコンをクリックして応援してください
