算額（その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;
