算額(その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
※コメント投稿者のブログIDはブログ作成者のみに通知されます