算額（その36）

算額（その36）

新潟県三島郡 与板八幡宮 寛政12年
新潟県柏崎市与板 与板三柱神社 であろう
http://www.wasan.jp/niigata/yoitahatiman2.html

直線と互いに接する大円と小円があり，更に小円と大円と直線に接する連続する小円がある。大円の径が225寸，一番小さい小円の径が4寸のとき，赤で示した面積が最大になるときの一番大きい小円の径はいくつになるか。

答えを先に述べておく。

算額によれば，径 4 寸から径 100 寸の最大の小円まで全部で 5 個とあるが，本当は 7 個ではないか?

この問題は，算額（その35）の後半に示したプログラムで解くことができる。
大円の半径を R，小さい方から 2 つの小円の半径と中心の x 座標を r1, x1, r2, x2 と置き，以下の方程式を立てる。

using SymPy
@syms R::positive, r1::positive, r2::positive, x1::positive, x2::positive;

R = 225
eq1 = x1^2 + (R - r1)^2 - (R + r1)^2;
eq2 = x2^2 + (R - r2)^2 - (R + r2)^2;
eq3 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2;

solve(eq1, x1) |> println

   Sym[2*sqrt(R)*sqrt(r1)]

res = solve([eq1, eq2, eq3], (r1, x1));
res[2][1] |> simplify |> println
res[2][2] |> simplify |> println

   x2^2*(2*sqrt(R)*sqrt(r2) + R + r2)/(4*(R - r2)^2)
   x2*(sqrt(R)*sqrt(r2) + R)/(R - r2)

SymPy はここまで。JupyterLab などを使っている場合は，ここでカーネルをリスタートする。

次に大きい円を求める関数を定義する。

function nextcircle(r2, x2; R = 225)
   r = x2^2*(2*sqrt(R)*sqrt(r2) + R + r2)/(4*(R - r2)^2)
   x = x2*(sqrt(R)*sqrt(r2) + R)/(R - r2)
   return r, x
end;

この漸化式によれば，半径が 4 寸からスタートして 7 番目の小円の半径がほぼ 100 になることが確認できる。

算額によれば，径 4 寸から径 100 寸の最大の小円まで全部で 5 個とあるが，本当は 7 個ではないか?

R = 225
r = 4 # r1
x = 2*sqrt(R)*sqrt(r) # x1

for i = 2:7
   r, x = nextcircle(r, x, R = 225)
   println("$i, $r, $x")
end

   2, 5.325443786982248, 69.23076923076923
   3, 7.438016528925618, 81.81818181818181
   4, 11.111111111111112, 99.99999999999999
   5, 18.3673469387755, 128.57142857142856
   6, 35.99999999999999, 179.99999999999997
   7, 99.99999999999996, 299.99999999999994

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 225
   plot()
   circle(0, R, R, :green)
   r = 4 # r1
   x = 2*sqrt(R)*sqrt(r) # x1
   circle(x, r, r)
   println("1, r = $r, x = $x")

   for i = 2:7
       r, x = nextcircle(r, x, R=R)
       println("$i, r = $r, x = $x")
       circle(x, r, r)
   end
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, r3, "r3 ", :black, :right)
       point(0, r3+r1, "r3+r1 ", :blue, :right)
       point(0, 1-r2, "1-r2 ", :red, :right)
       point(x4, r4, "(x4,r4)", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

draw(false)

   1, r = 4, x = 60.0
   2, r = 5.325443786982248, x = 69.23076923076923
   3, r = 7.438016528925618, x = 81.81818181818181
   4, r = 11.111111111111112, x = 99.99999999999999
   5, r = 18.3673469387755, x = 128.57142857142856
   6, r = 35.99999999999999, x = 179.99999999999997
   7, r = 99.99999999999996, x = 299.99999999999994