算額（その482）

算額（その482）

宮城県丸森町小斎日向 鹿島神社 明治13年

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

外円内に甲円 1 個と，乙円，大円，小円，1/3 円弧が各 2 個入っている。
小円の直径が 111 寸のときに，乙円の直径を求めよ。

1/3 円の半径は外円の半径と同じ長さである。
外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
乙円の半径と中心座標を r2, (x2, y2)
大円の半径と中心座標を r3, (x3,y3)
小円の半径と中心座標を r4, (x4, y4)
1/3円の半径と中心座標を r0, (x0, y0), (-x0, y0); x0 = r0*cosd(30), y0 = -r0*sind(30)
とおき，以下の連立方程式を nlsolve() で解き，数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, x0::positive, y0::positive, r1::positive,
      r2::positive, x2::positive, y2::positive,
      r3::positive, x3::positive, y3::negative,
      r4::positive, x4::positive, y4::negative;

(x0, y0) = r0 .* (cosd(Sym(30)), -sind(Sym(30)))
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 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq5 = x3^2 + (r0 - r1 - y3)^2 - (r1 + r3)^2
eq6 = (x2 - x0)^2 + (y2 - y0)^2 - (r2 + r0)^2
eq7 = (x0 - x3)^2 + (y0 - y3)^2 - (r0 - r3)^2
eq8 = (x3 + x0)^2 + (y3 - y0)^2 - (r0 + r3)^2
eq9 = (x4 + x0)^2 + (y4 - y0)^2 - (r0 + r4)^2
eq10 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10], (r0, r1, r2, x2, y2, r3, x3, y3, x4, y4))

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r0, r1, r2, x2, y2, r3, x3, y3, x4, y4) = u
   return [
       x2^2 + y2^2 - (r0 - r2)^2,  # eq1
       x3^2 + y3^2 - (r0 - r3)^2,  # eq2
       x4^2 + y4^2 - (r0 - r4)^2,  # eq3
       x2^2 - (r1 + r2)^2 + (r0 - r1 - y2)^2,  # eq4
       x3^2 - (r1 + r3)^2 + (r0 - r1 - y3)^2,  # eq5
       (r0/2 + y2)^2 - (r0 + r2)^2 + (-sqrt(3)*r0/2 + x2)^2,  # eq6
       (-r0/2 - y3)^2 - (r0 - r3)^2 + (sqrt(3)*r0/2 - x3)^2,  # eq7
       (r0/2 + y3)^2 - (r0 + r3)^2 + (sqrt(3)*r0/2 + x3)^2,  # eq8
       (r0/2 + y4)^2 - (r0 + r4)^2 + (sqrt(3)*r0/2 + x4)^2,  # eq9
       -(r3 + r4)^2 + (x3 - x4)^2 + (y3 - y4)^2,  # eq10
   ]
end;

r4 = 111/2
iniv = [big"262.1", 104.8, 42.5, 147.2, 162.9, 121.0, 139.7, -20.2, 85.4, -188.1] .* (2r4/111)
res = nls(H, ini=iniv);
println(res);

   (BigFloat[262.0833333333333172109930603571207271282688760651230918115043361933301448197767, 104.8333333333333223290545431295149542914840906048369827515914261371678084910254, 42.49999999999999786805343180884086908114509397785038194649776129778549724683521, 147.2243186433545568780056630764879738313623588458994394289193799498887306062796, 162.9166666666666596227807113723227255431557568042764514701570712819220452794033, 120.9615384615384515366858225567229155263396173464548862680437593159638153814604, 139.6743535847209988034369530154419060452378981480626331409871051079048418882699, -20.16025641025639136090801019725193623031123311832732180465654678920082881995361, 85.44783984006462676161828043480466817900534748463440517739428912144201840260664, -188.0833333333333094401143692287399234355784895913830835735564890504723176520062], true)

乙円の直径は 85 寸である。

   r0 = 262.083;  r1 =　104.833
   r2 =　42.5;  x2 =　147.224;  y2 = 162.917
   r3 =　120.962;  x3 =　139.674;  y3 =　-20.1603
   r4 = 55.5;  x4 = 85.4478;  y4 = -188.083
   乙円の直径 = 85

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r4 = 111/2
   (r0, r1, r2, x2, y2, r3, x3, y3, x4, y4)  = res[1]
   (x0, y0) = r0 .* (cosd(30), -sind(30))
   @printf("r0 = %g;  r1 =　%g\nr2 =　%g;  x2 =　%g;  y2 = %g\nr3 =　%g;  x3 =　%g;  y3 =　%g\nr4 = %g;  x4 = %g;  y4 = %g\n",  r0, r1, r2, x2, y2, r3, x3, y3, r4, x4, y4)
   @printf("乙円の直径 = %g\n", 2r2)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r0 - r1, r1, :green)
   circle(x2, y2, r2, :magenta)
   circle(x3, y3, r3, :blue)
   circle(-x3, y3, r3, :blue)
   circle(x4, y4, r4, :orange)
   circle(-x4, y4, r4, :orange)
   circle(x0, y0, r0, :red, beginangle=90, endangle=210)
   circle(-x0, y0, r0, :red, beginangle=-30, endangle=90)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r0 - r1, "甲円:r1 \n(0,r0-r1) ", :green, :right, :vcenter)
       point(x2, y2, "乙円:r2,(x2,y2) ", :magenta, :right, :vcenter)
       point(x3, y3, "大円:r3,(x3,y3)", :blue, :center, :top, delta=-delta)
       point(x4, y4, " 小円:r4,(x4,y4)", :orange, :left, :vcenter)
       point(x0, y0, "1/3円:r0,(x0,y0)  ", :red, :right, :top, delta=-1.5delta)
   else
       plot!(showaxis=false)
   end
end;