算額（その173）

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