算額（その3）

算額（その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