算額（その466）

算額（その466）

愛媛県松山市桜谷町 伊佐爾波神社 嘉永3年(1850)

平田浩一: 伊佐爾波神社の算額にみる江戸末期の和算，愛媛大学教育学部紀要，第60巻，195-206, 2013.
https://www.ed.ehime-u.ac.jp/~kiyou/2013/pdf/20.pdf

円弧内に青，黃，赤，白，黒の 5 個の円がある。青，赤，黒の 3 個の円の直径が与えられたとき，白円の直径を求めよ。

円弧の半径と中心座標を r0, (x0, y0)
青円の半径と中心座標を r1, (x1, y1)
黃円の半径と中心座標を r2, (x2, y2)
赤円の半径と中心座標を r3, (x3, y3); y3 = r3
白円の半径と中心座標を r4, (x4, y4); y4 = r4
黒円の半径と中心座標を r5, (x5, y5); x5 = 0, y5 = r5
とおき，以下の連立方程式を解き，数値解を求める。

include("julia-source.txt")

using SymPy

@syms r1::positive, x1::positive, y1::positive, 
     r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive,
     r4::positive, x4::positive,
     r5::positive,
     r0::positive, x0::positive, y0::positive;

x5 = 0
y5 = r5
y3 = r3
y4 = r4
eq1 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq3 = (x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq4 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
eq6 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq7 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
eq8 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq9 = (x0 - x1)^2 + (y0 - y1)^2 - (r0 - r1)^2
eq10 = (x0 - x3)^2 + (y0 - y3)^2 - (r0 - r3)^2
eq11 = (x0 - x5)^2 + (y0 - y5)^2 - (r0 - r5)^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11], (x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0))

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0) = u
   return [
-(r1 + r2)^2 + (x1 - x2)^2 + (y1 - y2)^2,  # eq1
-(r1 + r3)^2 + (-r3 + y1)^2 + (x1 - x3)^2,  # eq2
x1^2 - (r1 + r5)^2 + (-r5 + y1)^2,  # eq3
-(r2 + r3)^2 + (-r3 + y2)^2 + (x2 - x3)^2,  # eq4
(r3 - r4)^2 - (r3 + r4)^2 + (x3 - x4)^2,  # eq5
-(r2 + r4)^2 + (-r4 + y2)^2 + (x2 - x4)^2,  # eq6
x2^2 - (r2 + r5)^2 + (-r5 + y2)^2,  # eq7
x4^2 + (r4 - r5)^2 - (r4 + r5)^2,  # eq8
-(r0 - r1)^2 + (x0 - x1)^2 + (y0 - y1)^2,  # eq9
-(r0 - r3)^2 + (-r3 + y0)^2 + (x0 - x3)^2,  # eq10
x0^2 - (r0 - r5)^2 + (-r5 + y0)^2,  # eq11
   ]
end;

青，赤，黒の 3 個の円の直径がそれぞれ 15, 30, 18 のときの解を求める。

(r1, r3, r5) = (15, 30, 18)
iniv = (r1/31) .* [big"15.0", 45, 6, 21, 25, 62, 10, 25, 57, 37, -15]
iniv = r1 .* [big"7.258", 21.774, 2.903, 10.161, 12.097, 30.0, 4.839, 12.097, 27.581, 17.903, -7.258]
res = nls(H, ini=iniv)
names = ("x1", "y1", "r2", "x2", "y2", "x3", "r4", "x4", "r0", "x0", "y0")
if res[2]
   for (name, x) in zip(names, res[1])
       @printf("%s = %g;  ", name, round(Float64(x), digits=6))
   end
   println()
else
   println("収束していない")
end;
@printf("白円の直径は %g である。\n", 2res[1][7])

   x1 = 17.5989;  y1 = 45.9156;  r2 = 6.54237;  x2 = 23.4663;  y2 = 25.1876;  x3 = 59.6904;  r4 = 9.42818;  x4 = 26.0544;  r0 = 59.7271;  x0 = 40.362;  y0 = 7.41423;  

白円の直径は 18.8564 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0) = res[1]
   @printf("x1 = %g;  y1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  r4 = %g;  x4 = %g;  r0 = %g;  x0 = %g;  y0 = %g\n", x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0)
   @printf("白円の直径は %g である。\n", 2r4)
   plot()
   circle(x1, y1, r1)
   circle(x2, y2, r2)
   circle(x3, r3, r3)
   circle(x4, r4, r4)
   circle(0, r5, r5)
   a = sqrt(r0^2 - y0^2)
   segment(x0 - a, 0, x0 + a, 0, :magenta, lw=0.5)
   θ = atand(y0/a)
   circle(x0, y0, r0, :magenta, beginangle=-θ, endangle=181 + θ)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
        point(x1, y1, " 青:r1\n (x1,y1)", :blue, :center, :bottom, delta=delta)
        point(x2, y2, " 黃:r2,(x2,y2)", :black, :left, :vcenter)
        point(x3, r3, " 赤:r3,(x3,r3)", :red, :center, :bottom, delta=delta)
        point(x4, r4, "白:r4,(x4,r4)", :black, :center, :top, delta=-delta)
        point(x5, r5, "黒:r5,(x5,r5)", :black, :center, :top, delta=-delta)
        point(x0, y0, " 円弧:r0,(x0,y0)", :magenta, :left, :bottom, delta=delta)
     hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;