算額（その111）

算額（その111）

愛媛県松山市 伊佐爾波神社 嘉永3年(1850)正月
http://www.wasan.jp/ehime/isaniwa7.html

円内に 2 つの甲円，3 つの乙円が入っている。それぞれの円の半径を求めよ。

外円の半径を 1 とする。甲円，乙円の中心座標と半径を (r1, y1, r1)，(r2, y2, r2) とおき，以下の方程式を立てて解く。

using SymPy
@syms r1, r2, y1, y2;
eq3 = r1^2 + (1 - r2 - y1)^2 - (r1 + r2)^2
eq4 = (r1 - r2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq5 = r1^2 + y1^2 - (1 - r1)^2
eq6 = r2^2 + y2^2 - (1 - r2)^2;

res = solve([eq3, eq4, eq5, eq6], (r1, r2, y1, y2))

   5-element Vector{NTuple{4, Sym}}:
    (0, 0, 1, 1)
    ((-142*sqrt(2) - 197 + 65*sqrt(-13 + 16*sqrt(2)) + 46*sqrt(-26 + 32*sqrt(2)))/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -16*sqrt(2)/49 + 2*sqrt(-13 + 16*sqrt(2))/49 + 13/49 + 4*sqrt(-26 + 32*sqrt(2))/49, (-10*sqrt(-26 + 32*sqrt(2)) - 13*sqrt(-13 + 16*sqrt(2)) + 30*sqrt(2) + 43)/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -4*sqrt(2)/7 - 2/7 + sqrt(-13 + 16*sqrt(2))/7)
    (-(197 + 142*sqrt(2) + 65*sqrt(-13 + 16*sqrt(2)) + 46*sqrt(-26 + 32*sqrt(2)))/(sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -16*sqrt(2)/49 - 4*sqrt(-26 + 32*sqrt(2))/49 - 2*sqrt(-13 + 16*sqrt(2))/49 + 13/49, (13*sqrt(-13 + 16*sqrt(2)) + 30*sqrt(2) + 43 + 10*sqrt(-26 + 32*sqrt(2)))/(sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -4*sqrt(2)/7 - sqrt(-13 + 16*sqrt(2))/7 - 2/7)
    ((-65*sqrt(13 + 16*sqrt(2)) + 46*sqrt(26 + 32*sqrt(2)) - 142*sqrt(2)*I + 197*I)/(sqrt(13 + 16*sqrt(2)) - 9*I + 4*sqrt(2)*I), 13/49 + 16*sqrt(2)/49 - 2*sqrt(-16*sqrt(2) - 13)/49 + 4*sqrt(-32*sqrt(2) - 26)/49, (-10*sqrt(26 + 32*sqrt(2)) + 13*sqrt(13 + 16*sqrt(2)) - 43*I + 30*sqrt(2)*I)/(sqrt(13 + 16*sqrt(2)) - 9*I + 4*sqrt(2)*I), -2/7 + 4*sqrt(2)/7 - I*sqrt(13 + 16*sqrt(2))/7)
    ((-65*sqrt(13 + 16*sqrt(2)) + 46*sqrt(26 + 32*sqrt(2)) - 197*I + 142*sqrt(2)*I)/(sqrt(13 + 16*sqrt(2)) - 4*sqrt(2)*I + 9*I), 13/49 + 16*sqrt(2)/49 - 4*sqrt(-32*sqrt(2) - 26)/49 + 2*sqrt(-16*sqrt(2) - 13)/49, (-10*sqrt(26 + 32*sqrt(2)) + 13*sqrt(13 + 16*sqrt(2)) - 30*sqrt(2)*I + 43*I)/(sqrt(13 + 16*sqrt(2)) - 4*sqrt(2)*I + 9*I), -2/7 + 4*sqrt(2)/7 + I*sqrt(13 + 16*sqrt(2))/7)

5 組の解が得られるが，2番目のものだけが適切である。なお，@syms において各変数を ::positive と宣言すると解が求まらないので注意。

for i = 1:5
   println("r1 = $(res[i][1].evalf())")
   println("r2 = $(res[i][2].evalf())")
   println("y1 = $(res[i][3].evalf())")
   println("y2 = $(res[i][4].evalf())")
   println()
end

   r1 = 0
   r2 = 0
   y1 = 1.00000000000000
   y2 = 1.00000000000000

   r1 = 0.494520180083009
   r2 = 0.288374102478171
   y1 = 0.104688298457764
   y2 = -0.650578046850383

   r1 = -45.1219371780525
   r2 = -0.681329898313661
   y1 = 9.55216595103462
   y2 = -1.53709459586173

   r1 = 0.31370849898476 - 0.463999436044365*I
   r2 = 0.727090142815704 + 0.445454898991966*I
   y1 = -0.82842712474619 - 0.560096865715887*I
   y2 = 0.522407749927483 - 0.852695809075959*I

   r1 = 0.31370849898476 + 0.463999436044365*I
   r2 = 0.727090142815704 - 0.445454898991966*I
   y1 = -0.82842712474619 + 0.560096865715887*I
   y2 = 0.522407749927483 + 0.852695809075959*I

using Plots
using Printf

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; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")

   (r1, r2, y1, y2) = (
       (-142*sqrt(2) - 197 + 65*sqrt(-13 + 16*sqrt(2)) + 46*sqrt(-26 + 32*sqrt(2)))/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9),
       -16*sqrt(2)/49 + 2*sqrt(-13 + 16*sqrt(2))/49 + 13/49 + 4*sqrt(-26 + 32*sqrt(2))/49,
       (-10*sqrt(-26 + 32*sqrt(2)) - 13*sqrt(-13 + 16*sqrt(2)) + 30*sqrt(2) + 43)/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9),
       -4*sqrt(2)/7 - 2/7 + sqrt(-13 + 16*sqrt(2))/7)

   plot()
   circle(0, 0, 1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   circle(0, 1 - r2, r2, :blue)
   circle(r1, y1, r1, :green)
   circle(-r1, y1, r1, :green)
   if more == true
       @printf("r1 = %.5f; r2 = %.5f;  y1 = %.5f;  y2 = %.5f", r1, r2, y1, y2)
       point(r1, y1, "(r1,y1)", :green, :center)
       point(0, 1 - r2, "1-r2 ", :blue, :right)
       point(r2, y2, "(r2,y2)", :blue, :center)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.49452; r2 = 0.28837;  y1 = 0.10469;  y2 = -0.65058