算額（その3）
一点で接する半径の等しい弧がある。その弧にはさまれて，それぞれの弧に接し，また隣の円に接するような円を順々に作る。最初の円の直径を四寸，次は五寸とする。黒い部分（最小の円と両弧の間）の面積を最大にした場合，円をいくつはさむことができるか。
「黒い部分（最小の円と両弧の間）の面積を最大にした場合」というのはいらないのではないか？
- i 番目の円の半径 ri，中心座標 (xi, 0)
- i+1 番目の円の半径 rip1，中心座標 (xip1, 0)
- 円の半径を r，中心座標を (0, r) とする。その円の一部が「弧」
using SymPy
@syms r, ri, rip1, xi, xip1, r1::positive, r2::positive, x1::positive, x2::positive
(r, ri, rip1, xi, xip1, r1, r2, x1, x2)
中心座標と半径の関係
xpi1 - xi = rip1 + ri
xpi1 - xi - rip1 - ri = 0 …… （1）
弧が i 番目の円と接することから
(r + ri)^2 = r^2 + xi^2
xi = sqrt(2r * ri + ri^2)
xip1 = sqrt(2r * rip1 + rip1^2) …… （2）
（1) に (2) を代入
a = sqrt(2r * rip1 + rip1^2) - sqrt(2r * ri + ri^2) - rip1 - ri
a |> println
-ri - rip1 - sqrt(2*r*ri + ri^2) + sqrt(2*r*rip1 + rip1^2)
rip1 について解く。解が円の半径の漸化式になる。
b = solve(a, rip1)[1];
b |> println
ri*(r + ri + sqrt(ri*(2*r + ri)))/(r - ri - sqrt(ri*(2*r + ri)))
一番小さい円の半径を r1，2 番目に大きい円の半径 r2 とする。
eq1 = (r + r1)^2 - r^2 - x1^2;
eq2 = (r + r2)^2 - r^2 - x2^2;
eq3 = x2 - x1 - r2 - r1;
results = solve([eq1, eq2, eq3], (x1, x2, r));
r1 = 4, r2 = 5 が与えられているので，弧を持つ円の半径は 720 である。
results[1][3](r1 => 4, r2 => 5) |> println
720
---
using Plots
# 次に大きい円の半径
f(ri, r) = ri*(r + ri + sqrt(ri*(2*r + ri)))/(r - ri - sqrt(ri*(2*r + ri)))
function circle2(ox, oy, r; color=:red)
theta = 0:0.01:2pi
x = r.*cos.(theta)
y = r.*sin.(theta)
plot!(ox .+ x, oy .+ y, color=color)
end;
function draw(r1, r2)
r = 4*r1*r2*(r1 + r2)/(r1 - r2)^2
println("r = $r")
gr(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
plot()
circle2(0, r, r, color=:black)
circle2(0, -r, r, color=:black)
ri = r1
maxr, maxx = r1, 0
no = 0
while r - ri - sqrt(ri*(2*r + ri)) >= 0
x = sqrt((r + ri)^2 - r^2) # 半径 ri の円の中心の x 座標
maxx = max(x, maxx)
maxr = max(ri, maxr)
no += 1
println("$no: x = $x, r = $ri, $(ri + sqrt(ri*(2*r + ri)))")
circle2(x, 0, ri)
plot!([0, 0], [ri, 0], color=:blue)
ri = f(ri, r) # その次に大きい円の半径
end
println(maxx/maxr)
plot!(xlims=(0, maxx + maxr), ylims=(-maxr, maxr),
size=(500, 500 * 2maxr / (maxx + maxr)))
savefig("$r1-$r2.png")
end
答え：9個の円が存在する
draw(4, 5)
r = 720.0
1: x = 76.0, r = 4, 80.0
2: x = 85.0, r = 5.0, 90.0
3: x = 96.42857142857166, r = 6.428571428571429, 102.85714285714286
4: x = 111.42857142857119, r = 8.571428571428573, 120.00000000000001
5: x = 132.0, r = 12.000000000000002, 144.0
6: x = 162.0, r = 18.000000000000004, 180.0
7: x = 210.0, r = 30.000000000000007, 240.00000000000003
8: x = 300.0, r = 60.000000000000014, 360.0
9: x = 540.0, r = 180.00000000000003, 720.0
2.9999999999999996
その他の例
draw(5, 6)
r = 1320.0
1: x = 115.0, r = 5, 120.0
2: x = 126.0, r = 6.0, 132.0
3: x = 139.3333333333324, r = 7.333333333333333, 146.66666666666669
4: x = 155.833333333334, r = 9.166666666666664, 164.99999999999997
5: x = 176.78571428571388, r = 11.785714285714283, 188.57142857142853
6: x = 204.2857142857146, r = 15.714285714285708, 219.99999999999994
7: x = 242.0, r = 21.99999999999999, 263.99999999999994
8: x = 297.0, r = 32.999999999999986, 329.99999999999994
9: x = 385.0, r = 54.99999999999998, 439.99999999999994
10: x = 550.0, r = 109.99999999999996, 659.9999999999999
11: x = 989.9999999999995, r = 329.99999999999983, 1319.9999999999995
3.0
draw(8,9)
r = 4896.0
1: x = 280.0, r = 8, 288.0
2: x = 297.0, r = 9.0, 306.0
3: x = 316.19999999999624, r = 10.2, 326.4
4: x = 338.0571428571484, r = 11.657142857142855, 349.71428571428567
5: x = 363.1648351648321, r = 13.45054945054945, 376.6153846153846
6: x = 392.3076923076906, r = 15.692307692307693, 408.0
7: x = 426.545454545461, r = 18.545454545454547, 445.0909090909091
8: x = 467.3454545454534, r = 22.254545454545465, 489.6000000000001
9: x = 516.7999999999984, r = 27.200000000000006, 544.0000000000001
10: x = 578.0, r = 34.00000000000001, 612.0
11: x = 655.7142857142836, r = 43.71428571428572, 699.4285714285716
12: x = 757.7142857142878, r = 58.28571428571428, 816.0
13: x = 897.600000000003, r = 81.6, 979.2
14: x = 1101.5999999999976, r = 122.4, 1224.0000000000002
15: x = 1428.0, r = 204.0, 1632.0
16: x = 2040.0, r = 408.0, 2448.0
17: x = 3672.0, r = 1224.0, 4896.0
3.0
draw(99,100)
r = 7.8804e6
1: x = 39501.0, r = 99, 39600.0
2: x = 39700.0, r = 100.0, 39800.0
3: x = 39901.01522841007, r = 101.01522842639594, 40002.03045685279
4: x = 40104.07645290719, r = 102.04599606339997, 40206.12244897959
5: x = 40309.21507052156, r = 103.09262166405024, 40412.307692307695
6: x = 40516.463124478025, r = 104.1554321966693, 40620.618556701025
7: x = 40725.853319693466, r = 105.2347631002617, 40831.08808290155
8: x = 40937.41904142758, r = 106.33095854922276, 41043.74999999999
9: x = 41151.194371621066, r = 107.44437172774865, 41258.638743455485
中略
198: x = 3.283499999999912e6, r = 656699.9999999664, 3.940199999999879e6
199: x = 5.910299999999696e6, r = 1.9700999999998182e6, 7.880399999999516e6
3.000000000000123
