算額(その173)
岐阜県大垣市外野釜笛 釜笛八幡神社 慶応元年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF
第6問: 外円内にこの円と同じ半径の円弧を上下左右に作り,その隙間に 5 個の円を入れる。外円の径を知って個々の累円の径を求めよ。
外円の径を R とし,累円の径を r1, r2, r3, r4, r5 とおき,連立方程式を解く。
using SymPy
@syms R::positive, r1::positive, r2::positive, r3::positive, r4::positive, r5::positive;
eq1 = R^2 + (R - r1)^2 - (R + r1)^2
eq2 = R^2 + (R - 2r1 - r2)^2 - (R + r2)^2
eq3 = R^2 + (R - 2r1 - 2r2 - r3)^2 - (R + r3)^2
eq4 = R^2 + (R - 2r1 - 2r2 - 2r3 - r4)^2 - (R + r4)^2
eq5 = R^2 + (R - 2r1 - 2r2 - 2r3 - 2r4 - r5)^2 - (R + r5)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, r3, r4, r5))
1-element Vector{NTuple{5, Sym}}:
(R/4, R/12, R/24, R/40, R/60)
[4, 12, 24, 40, 60, ...] は 4[1, 3, 6, 10, 15, ...] で,これは一般項 4(n*(n+1)/2) = 2n*(n+1) であることはよく知られていると思うが,これを SymPy でやるとどうなるかをいかに示しておく。
一般項が何次式になるかわからない場合には次数を多めに取り,一般項を f として定義する。今回の場合は 3次式と仮定する。ついで,使用した次数分の項を terms ベクトルに用意する。f(n => i) - terms[i] を i = 1,2,3,4 とし,連立方程式を解く。
まずはまだるっこしい方法で,記述が面倒くさい。
@syms a, b, c, d, n
f = a*n^3 + b*n^2 + c*n + d
terms = [4, 12, 24, 40]
coeffs = solve([f(n => 1) - terms[1], f(n => 2) - terms[2], f(n => 3) - terms[3], f(n => 4) - terms[4]], (a, b, c, d))
Dict{Any, Any} with 4 entries:
c => 2
d => 0
a => 0
b => 2
f(coeffs) |> simplify |> println
2*n*(n + 1)
内包表記を使えば重複する記述を避けることができる。
x = [f(n => i).evalf() - terms[i] for i in 1:4];
x |> println
Sym[a + b + c + d - 4, 8.0*a + 4.0*b + 2.0*c + d - 12, 27.0*a + 9.0*b + 3.0*c + d - 24, 64.0*a + 16.0*b + 4.0*c + d - 40]
f(solve(x)) |> simplify |> println
2.0*n*(n + 1)
外円の直径を R とすると n 番目の累円の直径は 1/(2n*(n+1))R/4 = R / (2n*(n+1))
[1/(2n*(n+1)) for n in 1:5]
5-element Vector{Float64}:
0.25
0.08333333333333333
0.041666666666666664
0.025
0.016666666666666666
一方,隣り合う累円について漸化式を立てれば以下のようになる。
@syms R, ri, rip1
a = sqrt(2R*ri + ri^2) - sqrt(2R*rip1 + rip1^2) - ri - rip1
a |> println
-ri - rip1 + sqrt(2*R*ri + ri^2) - sqrt(2*R*rip1 + rip1^2)
これを rip1 についてとくと,今の累円の半径 ri の次の累円の半径 rip1 を求める関数を定義できる。
solve(a, rip1)[1] |> println
ri*(R + ri - sqrt(ri*(2*R + ri)))/(R - ri + sqrt(ri*(2*R + ri)))
@syms R, r
f = r*(R + r - sqrt(r*(2R + r)))/(R - r + sqrt(r*(2R + r)))
f |> println
r*(R + r - sqrt(r*(2*R + r)))/(R - r + sqrt(r*(2*R + r)))
これを Julia の関数定義で書くと以下のようになる。
g(r, R) = r*(R + r - sqrt(r*(2R + r)))/(R - r + sqrt(r*(2R + r)));
外円の半径を R とすると,青円の半径は 1/4 なので,次の赤円の半径は以下のように 0.08333333333333333 となる。円の中心座標も計算できる。
g(1//4, 1)
0.08333333333333333
関数 g(r, R) を使って,順次累円を描くことができる。
using Plots
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
if fill
plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
else
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
end
end;
function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, 10, position, color, vertical))
end;
function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;
function draw2(i, r, x, R)
circle(R*cosd(90+60i), R*sind(90+60i), R, :black, beginangle=270+60i, endangle=330+60i)
circle(R*cosd(-30+60i), R*sind(-30+60i), R, :black, beginangle=90+60i, endangle=90+60+60i)
for j = 1:5
#println("x = $(x[j]*cosd(60i)), y = $(x[j]*sind(60i)), r = $(r[j])")
circle(x[j]*cosd(60i), x[j]*sind(60i), r[j], j, fill=true)
end
end
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
n = 5
r = zeros(n)
x = zeros(n)
R = 1
r[1] = R / 4
x[1] = R - r[1]
plot()
circle(0, 0, R, :black)
for i = 2:5
r[i] = g(r[i-1], R)
x[i] = x[i-1] - r[i-1] - r[i]
end
for i in 0:5
draw2(i, r, x, R)
end
if more == true
point(R, 0, " R", :black)
point(R - r[1], 0, " R-r1", :black, :left, :bottom)
point(R - 2r[1] - r[2], 0, " R-2r1-r2", :black)
segment(0, R, R - r[1], 0, :red, linestyle=:dash, linewidth=0.5)
segment(0, R, R - 2r[1] - r[2], 0, :red, linestyle=:dash, linewidth=0.5)
point(0, R, " R", :black)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
※コメント投稿者のブログIDはブログ作成者のみに通知されます