算額（その485）

算額（その485）

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

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

楕円の中に大円，中円，小円が入っている。楕円の短径(短軸の半分;図参照)が与えられたとき，小円の直径を求めよ。

大円の半径と中心座標を r1, (x1, 0), (-x1, 0)
中円の半径と中心座標を r2, (0, r2), (x2, r2)
小円の半径と中心座標を r3, (x3, y3)
楕円の長径と短径を a, b
楕円と右上の中円の接点を (x, y)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, 
     x::positive, y::positive,
     r1::positive, x1::positive,
     r2::positive, x2::positive,
     r3::positive, x3::positive, y3::positive;

b = 55
x3 = x1
x2 = 2x1
eq1 = (x3 - x1)^2 + y3^2 - (r1 - r3)^2
eq2 = x1^2 + r2^2 - (r1 - r2)^2
eq3 = (x3 + x1)^2 + y3^2 - (r1 + r3)^2
eq4 = (x1 + x2)^2 + r2^2 - (r1 + r2)^2
eq5 = (x2 - x3)^2 + (y3 - r2)^2 - (r2 + r3)^2
eq6 = (x - x1)^2 + y^2 - r1^2
eq7 = (x - x2)^2 + (y - r2)^2 - r2^2
eq8 = x^2/a^2 + y^2/b^2 - 1;
#res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (a, r1, x1, r2, r3, y3, x, y))

   y3^2 - (r1 - r3)^2,  # eq1
   r2^2 + x1^2 - (r1 - r2)^2,  # eq2
   4*x1^2 + y3^2 - (r1 + r3)^2,  # eq3
   r2^2 + 9*x1^2 - (r1 + r2)^2,  # eq4
   x1^2 + (-r2 + y3)^2 - (r2 + r3)^2,  # eq5
   -r1^2 + y^2 + (x - x1)^2,  # eq6
   -r2^2 + (-r2 + y)^2 + (x - 2*x1)^2,  # eq7
   y^2/3025 - 1 + x^2/a^2,  # eq8

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)
   (a, r1, x1, r2, r3, y3, x, y) = u
   return [
       y3^2 - (r1 - r3)^2,  # eq1
       r2^2 + x1^2 - (r1 - r2)^2,  # eq2
       4*x1^2 + y3^2 - (r1 + r3)^2,  # eq3
       r2^2 + 9*x1^2 - (r1 + r2)^2,  # eq4
       x1^2 + (-r2 + y3)^2 - (r2 + r3)^2,  # eq5
       -r1^2 + y^2 + (x - x1)^2,  # eq6
       -r2^2 + (-r2 + y)^2 + (x - 2*x1)^2,  # eq7
       y^2/3025 - 1 + x^2/a^2,  # eq8
   ]
end;

b = 55
iniv = [big"73.0", 49, 22, 19.5, 10, 39, 56, 35]
res = nls(H, ini=iniv)

   (BigFloat[70.61439286442122986759326634852342606346806573057936088754412457172871369525662, 48.1044778079481411118448683589097310128322176790819406320263186975907510564275, 21.51297648015799054903759651066595576338165507866319055547355868260632047583453, 19.24179112325079279837972287417723421441707964112536578143846733491061211029649, 9.62089556158975492218928774413918955995629190756242766925925034366901321701604, 38.48358224635801843450000033108172430091547897712345222108392054779560714508957, 57.36779066297351601295569753085714297026688222855118678814873700650448108748758, 32.06981579433871183683575416702468448547226311180548288477472348488073672944933], false)

短径が 55 のとき，
a = 70.6144;  r1 = 48.1045;  x1 = 21.513;  r2 = 19.2418;  r3 = 9.6209;  y3 = 38.4836;  x = 57.3678;  y = 32.0698
小円の直径 = 19.2418
となる。
しかし，この結果は初期値によって微妙に変化する。条件式のどれかが悪いか，決定的な条件式が書けているのだと思われる。

術では 短径sqrt(3)/5 となっており，短径 = 55 のときには小円の直径は 19.05255888325765 になる。

55sqrt(3)/5

   19.05255888325765

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 55
   (a, r1, x1, r2, r3, y3, x, y) = res[1]
   x3 = x1
   x2 = 2x1
   @printf("a = %g;  r1 = %g;  x1 = %g;  r2 = %g;  r3 = %g;  y3 = %g;  x = %g;  y = %g\n",
       a, r1, x1, r2, r3, y3, x, y)
   @printf("小円の直径 = %g\n", 2r3)
   plot()
   circle(x1, 0, r1, :black)
   circle(-x1, 0, r1, :black)
   circle4(x2, r2, r2)
   circle(0, r2, r2)
   circle(0, -r2, r2)
   circle4(x3, y3, r3, :blue)
   ellipse(0, 0, a, b, color=:green)
   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(x, y, "(x,y)", :green, :left, :bottom, delta=delta)
       point(x1, 0, "x1", :black, :center, delta=-delta)
       point(-x1, 0, "-x1", :black, :center, delta=-delta)
       point(x2, 0, "x2", :red, :center, delta=-delta)
       point(0, r2, " r2", :red, :left, :vcenter)
       point(x3, y3, " 小円:r3\n (x3,y3)", :blue, :left, :vcenter)
       point(x2, r2, "中円:r2(x2,r2)", :red, :center, delta=-delta)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(x1 + r1, 0, "x1+r1 ", :black, :right, :top, delta=-delta/2)
       point(0, b, " b", :green, :left, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;