算額（その438）

算額（その438）

社丈右衛門正常 天明2年壬寅4月(1782)

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

交差する 2 個の大円の中に，等円 3 個，甲円，乙円，丙円がそれぞれ 4 こずつ入っている。大円の直径が与えられたとき，大円の直径が与えられたとき，各円の直径を求めよ。

大円の半径と中心座標を 2r, (0, r), (0, -r)
等円の半径と中心座標を r, (0, 2r), (0, 0), (0, -2r)
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
として以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, r::positive,
     r1::positive, x1::positive, y1::positive,
     r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive, y4::positive,
     r5::positive, x5::positive, y5::positive,
     r6::positive, x6::positive, y6::positive;

r = r0//2
eq1 = x1^2 + (y1 + r)^2 - (2r + r1)^2
eq2 = x1^2 + (y1 - r)^2 - (2r - r1)^2
eq3 = x1^2 + (2r - y1)^2 - (r1 + r)^2

eq4 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq5 = x2^2 + (r - y2)^2 - (2r - r2)^2
eq6 = x2^2 + (y2 + r)^2 - (2r + r2)^2

eq7 = x3^2 + (r - y3)^2 - (2r - r3)^2
eq8 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq9 = x3^2 + (y3 + r)^2 - (2r + r3)^2;

solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (r1, x1, y1, r2, x2, y2, r3, x3, y3))

   2-element Vector{NTuple{9, Sym}}:
    (3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5, 9*r0/82, 20*sqrt(3)*r0/41, 9*r0/41, 27*r0/730, 182*sqrt(3)*r0/365, 27*r0/365)
    (3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5, 9*r0/82, 20*sqrt(3)*r0/41, 9*r0/41, 3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5)

最初の組が適解である。
(3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5, 9*r0/82, 20*sqrt(3)*r0/41, 9*r0/41, 27*r0/730, 182*sqrt(3)*r0/365, 27*r0/365)

r0 = 36; r = 18;  r1 = 10.8; x1 = 24.9415; y1 = 21.6; r2 = 3.95122; x2 = 30.4165; y2 = 7.90244; r3 = 1.33151; x3 = 31.0915; y3 = 2.66301

円の直径だけをいえば，甲円，乙円，丙円の直径は大円の直径のそれぞれ 3/10, 9/82, 27/730 倍である。

大円の直径 = 72 のとき，等円の直径 = 36;  甲円の直径 = 21.6;  乙円の直径 = 7.90244;  丙円の直径 = 2.66301

72 .* [3/10, 9/82, 27/730] |> println

   [21.599999999999998, 7.902439024390244, 2.663013698630137]

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 72/2
   r = r0/2
   (r1, x1, y1, r2, x2, y2, r3, x3, y3) = (3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5, 9*r0/82, 20*sqrt(3)*r0/41, 9*r0/41, 27*r0/730, 182*sqrt(3)*r0/365, 27*r0/365)
   @printf("r0 = %g; r = %g;  r1 = %g; x1 = %g; y1 = %g; r2 = %g; x2 = %g; y2 = %g; r3 = %g; x3 = %g; y3 = %g\n",
           r0, r, r1, x1, y1, r2, x2, y2, r3, x3, y3)
   @printf("大円の直径 = %g;  等円の直径 = %g;  甲円の直径 = %g;  乙円の直径 = %g;  丙円の直径 = %g\n", 2r0, 2r, 2r1, 2r2, 2r3)
   plot()
   circle(0, r, r0, :gray)
   circle(0, -r, r0, :gray)
   circle(0, 2r, r, :blue)
   circle(0, -2r, r, :blue)
   circle(0, 0, r, :blue)
   circle4(x1, y1, r1)
   circle4(x2, y2, r2, :blue)
   circle4(x3, y3, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, y1, "  甲円:r1,(x1,y1)", :black, :left, :vcenter)
       point(x2, y2, "   乙円:r2,(x2,y2)", :black, :left, :vcenter)
       point(x3, y3, "  丙円:r3,(x3,y3)", :black, :left, :vcenter)
       point(0, r, " r")
       point(0, -r, " -r")
       point(0, 2r, " 2r")
       point(0, 3r, " 3r")
       point(0.08r, 2.4r, "等円", :blue, mark=false)
       point(0.8x1, 2.1r, "大円", :gray, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その437）

算額（その437）

橘田彌曾八元克 天明八年戊申2月(1788)

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

大円の中に，甲，乙，丙，丁，戊，己の 6 個の円が入っている。丁円，戊円，己円の直径がそれぞれ 872.3 寸，671 寸，572 寸のとき，大円の直径はいかほどか。

（必要ならば）己円の中心がy 軸上に来るように回転する。中心の x 座標が 0 になる。
大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4)　r4 = 872.3/2
戊円の半径と中心座標を r5, (x5, y5)  r5 = 671/2
己円の半径と中心座標を r6, (x6, y6)  r6 = 572/2, x6 = 0
として以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r0,
      r1, x1, y1, r2, x2, y2, r3, x3, y3,
      r4, x4, y4, r5, x5, y5, r6, x6, y6

eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = x2^2 + y2^2 - (r0 - r2)^2
eq3 = x3^2 + y3^2 - (r0 - r3)^2
eq4 = x4^2 + y4^2 - (r0 - r4)^2
eq5 = x5^2 + y5^2 - (r0 - r5)^2
eq6 = x6^2 + y6^2 - (r0 - r6)^2
eq7 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq8 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq9 = (x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2
eq10 = (x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq11 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq12 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq13 = (x2 - x6)^2 + (y2 - y6)^2 - (r2 + r6)^2
eq14 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq15 = (x3 - x6)^2 + (y3 - y6)^2 - (r3 + r6)^2;

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)
   (r0, r1, x1, y1, r2, x2, y2, r3, x3, y3, x4, y4, x5, y5, y6) = u
   return [
       x1^2 + y1^2 - (r0 - r1)^2,  # eq1
       x2^2 + y2^2 - (r0 - r2)^2,  # eq2
       x3^2 + y3^2 - (r0 - r3)^2,  # eq3
       x4^2 + y4^2 - (r0 - r4)^2,  # eq4
       x5^2 + y5^2 - (r0 - r5)^2,  # eq5
       x6^2 + y6^2 - (r0 - r6)^2,  # eq6
       -(r1 + r2)^2 + (x1 - x2)^2 + (y1 - y2)^2,  # eq7
       -(r1 + r3)^2 + (x1 - x3)^2 + (y1 - y3)^2,  # eq8
       -(r1 + r4)^2 + (x1 - x4)^2 + (y1 - y4)^2,  # eq9
       -(r1 + r5)^2 + (x1 - x5)^2 + (y1 - y5)^2,  # eq10
       -(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2,  # eq11
       -(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2,  # eq12
       -(r2 + r6)^2 + (x2 - x6)^2 + (y2 - y6)^2,  # eq13
       -(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq14
       -(r3 + r6)^2 + (x3 - x6)^2 + (y3 - y6)^2,  # eq15
   ]
end;

(r4, r5, r6) = (872.3, 671, 572) ./2
x6 = 0
iniv = [big"95.0", 57, 12, 40, 41, -40, -43, 34, 36, -47, -72, 18, 75, -8, -79] .* (620.3/28)
res = nls(H, ini=iniv);
println([round(Float64(x), digits=6) for x in res[1]], " 収束:", res[2]);

   [1586.0, 830.761905, 265.84381, 706.902857, 726.916667, -668.763333, -539.24, 623.071429, 672.917143, -688.777143, -994.422, 577.304, 1248.06, 78.08, -1300.0] 収束:true

r0 = 1586
r1 = 830.762; x1 = 265.844;  y1 = 706.903
r2 = 726.917; x2 = -668.763; y2 = -539.24
r3 = 623.071; x3 = 672.917;  y3 = -688.777
r4 = 436.15;  x4 = -994.422; y4 = 577.304
r5 = 335.5;   x5 = 1248.06;  y5 = 78.08
r6 = 286;     x6 = 0;        y6 = -1300

大円の半径は 1586 寸（直径は 3172 寸）である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, r5, r6) = (872.3, 671, 572) ./2
   x6 = 0
   (r0, r1, x1, y1, r2, x2, y2, r3, x3, y3, x4, y4, x5, y5, y6) = res[1]
   @printf("r0 = %g; r1 = %g; x1 = %g; y1 = %g; r2 = %g; x2 = %g; y2 = %g; r3 = %g; x3 = %g; y3 = %g; r4 = %g;  x4 = %g; y4 = %g; r5 = %g;  x5 = %g; y5 = %g; r6 = %g;  x6 = %g;  y6 = %g\n", r0, r1, x1, y1, r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5, r6, x6, y6)
   @printf("大円の半径 = %g;  直径 = %g\n", r0, 2r0)
   plot()
   circle(0, 0, r0, :gray)
   circle(x1, y1, r1)
   circle(x2, y2, r2, :blue)
   circle(x3, y3, r3, :green)
   circle(x4, y4, r4, :magenta)
   circle(x5, y5, r5, :orange)
   circle(x6, y6, r6, :brown)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, y1, " 甲:r1,(x1,y1)", :red, :left, :vcenter)
       point(x2, y2, "乙:r2,(x2,y2)", :blue, :center, :top, delta=-delta)
       point(x3, y3, "丙:r3,(x3,y3)", :green, :center, :top, delta=-delta)
       point(x4, y4, "丁:r4,(x4,y4)", :magenta, :center, :top, delta=-delta)
       point(x5, y5, "戊:r5,(x5,y5)", :orange, :center, :bottom, delta=delta/2)
       point(x6, y6, " 己:r6,(x6,y6)", :black, :left, delta=-delta/2)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;