算額（その238）

算額（その238）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第6問　円の中に同じ半径で円弧を 2 個，その中心を通るように作り，その間に緑，赤，黄，浅青，黒円の合計 22 個の円を入れる。外側の円の直径を与えたとき，赤円の直径を求めよ。

この算額は，赤円の直径のみ求めよとなっているので，以下の算額（その55）愛宕神社のものと同じになっている。

宮城県角田市横倉 愛宕神社 明治15年(1882)1月
http://www.wasan.jp/miyagi/yokokuraatago.html

せっかくなので算額にある他の円の直径も求めてみよう。

原点を中心とする大円の半径を r0 とし，以下のように記号を定め方程式を解く。

include("julia-source.txt")

using SymPy
@syms r0, r1::positive, y1::positive, r2::positive, x2::positive, y2::positive, r3::positive;
@syms x3::positive, y3::positive;

#r0 = 1
eq1 = (r0//2)^2 +(r0//2 + r1)^2 - (r0 - r1)^2
eq2 = (r0 - r1)^2 + y1^2 - (r0 + r1)^2
res1 = solve([eq1, eq2], (r1, y1))

   2-element Vector{Tuple{Sym, Sym}}:
    (r0/6, -sqrt(6)*r0/3)
    (r0/6, sqrt(6)*r0/3)

赤円の次に大きい黄円の半径を r2, その中心座標を (x2, y2) とおき，次の連立方程式を解く。

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

r0 = 1
r1 = r0//6
(x1, y1) = (r0//2, r0//2 + r1)
eq3 = (r0//2 - x2)^2 + y2^2 - (r0//2 + r2)^2 |> simplify
eq4 = (r0 - x2)^2 + y2^2 - (r0 - r2)^2 |> simplify
eq5 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2 |> simplify
res = solve([eq3, eq4, eq5], (r2, x2, y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (1/11, 3/11, 6/11)
    (1/3, 1, 2/3)

最初の解が適切である。

次に，現在の円の r2, (x2, y2) がわかっているときに次の円の半径 r3，中心座標 (x3, y3) を求める漸化式を作る。上と同じ連立方程式（記号を変える）を解く。解は r2, (x2, y2) を含む式になる。

using SymPy
@syms r0, r1::positive, x1::positive, y2::positive, r2::positive, x2::positive, y2::positive,
   r3::positive, x3::positive, y3::positive;
r0 = 1
eq3 = (r0//2 - x3)^2 + y3^2 - (r0//2 + r3)^2 |> simplify
eq4 = (r0 - x3)^2 + y3^2 - (r0 - r3)^2 |> simplify
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> simplify
res = solve([eq3, eq4, eq5], (r3, x3, y3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    ((-r2^3 - 3*r2^2*x2 + 2*r2^2 + r2*x2^2 + r2*y2^2 + 3*x2^3 - 2*x2^2 + 3*x2*y2^2 + 2*y2^2 - 2*sqrt(2)*y2*sqrt(-r2^4 - r2^3 + 2*r2^2*x2^2 - 3*r2^2*x2 + 2*r2^2*y2^2 + 2*r2^2 + r2*x2^2 + r2*y2^2 - x2^4 + 3*x2^3 - 2*x2^2*y2^2 - 2*x2^2 + 3*x2*y2^2 - y2^4))/(2*(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4)), 3*(-r2^3 - 3*r2^2*x2 + 2*r2^2 + r2*x2^2 + r2*y2^2 + 3*x2^3 - 2*x2^2 + 3*x2*y2^2 + 2*y2^2 - 2*sqrt(2)*y2*sqrt(-r2^4 - r2^3 + 2*r2^2*x2^2 - 3*r2^2*x2 + 2*r2^2*y2^2 + 2*r2^2 + r2*x2^2 + r2*y2^2 - x2^4 + 3*x2^3 - 2*x2^2*y2^2 - 2*x2^2 + 3*x2*y2^2 - y2^4))/(2*(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4)), 2*y2*(-2*r2^2 - r2 + 2*x2^2 - 3*x2 + 2*y2^2 + 2)/(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4) + sqrt(2)*sqrt(-(-r2^2 - 2*r2 + x2^2 - 2*x2 + y2^2)*(-r2^2 + r2 + x2^2 - x2 + y2^2))*(r2 + 3*x2 - 2)/(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4))
    ((-r2^3 - 3*r2^2*x2 + 2*r2^2 + r2*x2^2 + r2*y2^2 + 3*x2^3 - 2*x2^2 + 3*x2*y2^2 + 2*y2^2 + 2*sqrt(2)*y2*sqrt(-r2^4 - r2^3 + 2*r2^2*x2^2 - 3*r2^2*x2 + 2*r2^2*y2^2 + 2*r2^2 + r2*x2^2 + r2*y2^2 - x2^4 + 3*x2^3 - 2*x2^2*y2^2 - 2*x2^2 + 3*x2*y2^2 - y2^4))/(2*(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4)), 3*(-r2^3 - 3*r2^2*x2 + 2*r2^2 + r2*x2^2 + r2*y2^2 + 3*x2^3 - 2*x2^2 + 3*x2*y2^2 + 2*y2^2 + 2*sqrt(2)*y2*sqrt(-r2^4 - r2^3 + 2*r2^2*x2^2 - 3*r2^2*x2 + 2*r2^2*y2^2 + 2*r2^2 + r2*x2^2 + r2*y2^2 - x2^4 + 3*x2^3 - 2*x2^2*y2^2 - 2*x2^2 + 3*x2*y2^2 - y2^4))/(2*(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4)), 2*y2*(-2*r2^2 - r2 + 2*x2^2 - 3*x2 + 2*y2^2 + 2)/(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4) - sqrt(2)*sqrt(-(-r2^2 - 2*r2 + x2^2 - 2*x2 + y2^2)*(-r2^2 + r2 + x2^2 - x2 + y2^2))*(r2 + 3*x2 - 2)/(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4))

これも，最初の式が適切である。この解に基づき，r2, x2, y2 を与えて r3, x3, y3 を得る関数にする。

nextcircle(r2, x2, y2) = ((-r2^3 - 3*r2^2*x2 + 2*r2^2 + r2*x2^2 + r2*y2^2 + 3*x2^3 - 2*x2^2 + 3*x2*y2^2 + 2*y2^2 - 2*sqrt(2)*y2*sqrt(-r2^4 - r2^3 + 2*r2^2*x2^2 - 3*r2^2*x2 + 2*r2^2*y2^2 + 2*r2^2 + r2*x2^2 + r2*y2^2 - x2^4 + 3*x2^3 - 2*x2^2*y2^2 - 2*x2^2 + 3*x2*y2^2 - y2^4))/(2*(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4)), 3*(-r2^3 - 3*r2^2*x2 + 2*r2^2 + r2*x2^2 + r2*y2^2 + 3*x2^3 - 2*x2^2 + 3*x2*y2^2 + 2*y2^2 - 2*sqrt(2)*y2*sqrt(-r2^4 - r2^3 + 2*r2^2*x2^2 - 3*r2^2*x2 + 2*r2^2*y2^2 + 2*r2^2 + r2*x2^2 + r2*y2^2 - x2^4 + 3*x2^3 - 2*x2^2*y2^2 - 2*x2^2 + 3*x2*y2^2 - y2^4))/(2*(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4)), 2*y2*(-2*r2^2 - r2 + 2*x2^2 - 3*x2 + 2*y2^2 + 2)/(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4) + sqrt(2)*sqrt(-(-r2^2 - 2*r2 + x2^2 - 2*x2 + y2^2)*(-r2^2 + r2 + x2^2 - x2 + y2^2))*(r2 + 3*x2 - 2)/(r2^2 + 6*r2*x2 - 4*r2 + 9*x2^2 - 12*x2 + 8*y2^2 + 4));

この関数を使えば r1 = 1/6 から順次 r2, r3,... を求めることができる。

r1 = 1/6 から始まる類円は，r2 = 1/11, r3 = 1/18, r4 = 1/27... である。

n 番目の類縁の半径の一般項は n^2 + 2n + 3 である。

各円の中心座標を求めて図を描く。

   r0 = 1.00000;  r1 = 0.16667;  x1 = 0.16667;  y1 = 0.81650
    1  (0.16666666666666666,  0.5,                  0.6666666666666666)   1/r1 = 6.0
    2  (0.0909090909090909,   0.27272727272727265,  0.5454545454545454)   1/r2 = 11.000000000000002
    3  (0.05555555555555554,  0.1666666666666666,   0.44444444444444436)  1/r2 = 18.000000000000007
    4  (0.037037037037037014, 0.11111111111111105,  0.37037037037037024)  1/r2 = 27.000000000000018
    5  (0.026315789473684195, 0.07894736842105259,  0.31578947368421045)  1/r2 = 38.00000000000002
    6  (0.01960784313725489,  0.05882352941176468,  0.2745098039215686)   1/r2 = 51.00000000000003
    7  (0.015151515151515148, 0.04545454545454544,  0.24242424242424243)  1/r2 = 66.00000000000001
    8  (0.012048192771084336, 0.03614457831325301,  0.21686746987951808)  1/r2 = 83.00000000000001
    9  (0.00980392156862745,  0.029411764705882353, 0.19607843137254902)  1/r2 = 102.0
   10  (0.008130081300813007, 0.024390243902439025, 0.17886178861788618)  1/r2 = 123.00000000000001
   11  (0.006849315068493151, 0.020547945205479454, 0.1643835616438356)   1/r2 = 146.0

using Plots
using Printf

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 1
   (r1, y1) = (r0/6, sqrt(6)*r0/3)
   @printf("r0 = %.5f;  r1 = %.5f;  x1 = %.5f;  y1 = %.5f\n", r0, r1, r1, y1)
   plot()
   circle(0, 0, r0, :green)
   circle(r0, 0, r0, :green, beginangle=120, endangle=240)
   circle(-r0, 0, r0, :green, beginangle=300, endangle=420)
   circle(r0/2, r0/2 + r1, r1, :magenta)
   circle(-r0/2, r0/2 + r1, r1, :magenta)
   circle(r0/2, -r0/2 - r1, r1, :magenta)
   circle(-r0/2, -r0/2 - r1, r1, :magenta)
   circle(r1, y1, r1, :magenta)
   circle(-r1, y1, r1, :magenta)
   circle(r1, -y1, r1, :magenta)
   circle(-r1, -y1, r1, :magenta)
   circle(r0/2, 0, 1/2, :blue)
   circle(-r0/2, 0, 1/2, :blue)
   if more
       point(r1, sqrt(6)*r0/3, "(r1,y1)", :magenta, :center)
       point(r0/2, r0/2+r1, "(r0/2,r0/2+r1)", :magenta, :center)
       point(r0/2, 0, " r0/2", :blue)
       point(0, 0, " 0", :black)
       point(r0, 0, " r0", :black)
       #point(x2, y2, "(x2,y2)")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   end
   (x1, y1) = (r0/2, r0/2 + r1)
   println("1  ", (r1, x1, y1), "  1/r1 =$(1/r1)")
   for i = 1:10
       (r2, x2, y2) = nextcircle(r1, x1, y1)
       println("$(i+1)  ", (r2, x2, y2), "  1/r2 =$(1/r2)")
       circle(x2, y2, r2, i)
      (r1, x1, y1) = (r2, x2, y2)
   end  
end;