算額（その352）
群馬県桐生市天神町 天満宮 文化6年(1830)
大山誠: 桐生市天満宮の算額題についてhttps://www.jstage.jst.go.jp/article/jjsme/80/7/80_20/_pdf
図のように大円と直角三角形が交わっており，等円を 3 個配置する。大円の直径が 10 寸のとき，等円の直径を求めよ。
大円の直径を r0 = 10 ，AC の長さを a，B の x 座標を b（b < 0)，等円の直径を r，その中心座標を (r0 - r, 0), (x1, r - a), (x2, y2) とおく。
以下の連立方程式を nlsolve() で解く。
eq6 は，∠EBF = θ とすると，∠OAD = π/4 - θ および tan(π/4 - θ) = (1 - tan(θ)/(1 + tan(θ) であることによる。
include("julia-source.txt");
using SymPy
@syms r0::positive, r::positive, a::positive, b::negative, x1::negative, x2::negative, y2::positive;
r0 = 10
eq1 = (r0 - 2r)^2 + a^2 - r0^2
eq2 = x2^2 + y2^2 - (r0 - r)^2
eq3 = x1^2 + (r - a)^2 - (r0 + r)^2
eq4 = distance(b, -a, r0 - 2r, a, x1, r - a) - r^2
eq5 = distance(b, -a, r0 - 2r, a, x2, y2) - r^2
tanθ = r/(x1 - b)
eq6 = (r0 - 2r)/a - (1 - tanθ)/(1 + tanθ); # (r0 - 2r)/a - tan(PI/4 - θ)
# θ =atan(tanθ)
# eq6 = (r0 - 2r)/a - tan(PI/4 - θ);
# res = solve([eq1, eq2, eq3, eq4, eq5], (r, a, b, x1, x2, y2))
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)
(r, a, b, x1, x2, y2) = u
return [
a^2 + (10 - 2*r)^2 - 100, # eq1
x2^2 + y2^2 - (10 - r)^2, # eq2
x1^2 + (-a + r)^2 - (r + 10)^2, # eq3
-r^2 + (x1 - (4*a^2*b - 2*a*b*r - 4*a*r^2 + 20*a*r + b^2*x1 + 4*b*r*x1 - 20*b*x1 + 4*r^2*x1 - 40*r*x1 + 100*x1)/(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100))^2 + (-a - a*(-4*a^2 + 4*a*r + b^2 - 2*b*x1 - 4*r^2 - 4*r*x1 + 40*r + 20*x1 - 100)/(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100) + r)^2, # eq4
-r^2 + (x2 - (b*(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100) - (b + 2*r - 10)*(2*a^2 + 2*a*y2 + b^2 + 2*b*r - b*x2 - 10*b - 2*r*x2 + 10*x2))/(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100))^2 + (-a*(4*a*y2 + b^2 - 2*b*x2 - 4*r^2 - 4*r*x2 + 40*r + 20*x2 - 100)/(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100) + y2)^2, # eq5
-(-r/(-b + x1) + 1)/(r/(-b + x1) + 1) + (10 - 2*r)/a, # eq6
]
end;
iniv = [big"2.8", 8.9, -19.5, -11.2, -4.3, 5.8]
res = nls(H, ini=iniv);
println(res);
(BigFloat[2.758765525590156335644381495793494653227916752992539435376878299670476979753896, 8.939097947940123521918150817467879027312574890368149878678356808089001454295962, -19.47025838954737912975466155729263249608233106765657839605148838383817497557, -11.16197065424549685348050573396238073769521442013558168263557635706947113807764, -4.33134044315998946461745333303927732461882711841070523331857400641089934088305, 5.803013585959302297582280566300900177092654477073446038766027387779471835660339], true)
r = 2.75877; a = 8.9391; b = -19.4703; x1 = -11.162; x2 = -4.33134; y2 = 5.80301
大円の直径が 10 寸のとき，等円の直径は 2.75877 寸である。
include("julia-source.txt");
using Plots
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r0 = 10
(r, a, b, x1, x2, y2) = res[1]
@printf("r = %g; a = %g; b = %g; x1 = %g; x2 = %g; y2 = %g\n", r, a, b, x1, x2, y2)
plot()
circle(0, 0, r0, :blue)
circle(x1, r - a, r)
circle(x2, y2, r)
circle(r0 - r, 0, r)
if more
plot!([r0 - 2r, r0 - 2r, b], [-a, a, -a], color=:black, lw=0.5)
segment(0, 0, r0 - 2r, a, :orange)
segment(b, -a, x1, r - a, :orange)
segment(x1, r - a, x1, -a, :orange)
point(r0, 0, " r0", :blue, :left, :bottom)
point(r0 - r, 0, "r0-r", :red, :center, :bottom)
point(r0 - 2r, 0, "r0-2r ", :red, :right, :bottom)
point(r0 - 2r, a, " A:(r0-2r,a)", :black, :left, :bottom)
point(b, -a, " B:(b,-a)", :black)
point(r0 - 2r, -a, " C:(r0-2r,-a)", :black, :left)
point(r0 - 2r, 0, "D ", :black, :right)
point(0, 0, " O", :black)
point(x1, r - a, "E:(x1,r-a)", :red, :center, :bottom, mark=false)
point(x1, -a, " F", :black)
point(x2, y2, "(x2,y2)", :red, :center)
hline!([0, -a], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!([b, r0 - 2r, r0 - 2r, b], [-a, -a, a, -a], color=:black, lw=0.5)
plot!(showaxis=false)
end
end;
注: r の解析解は以下の 3次式=0 の解の 10 倍である。
SymPy では虚数解になるが虚部は非常に小さく（実質0），3 個の実数解が得られる。この内適解は t = 0.275876552559016 なので，r = 2.75876552559016 である。
using SymPy
@syms t::positive
eq = 8t^3 - 8sqrt(Sym(2))t^2 +5t - (12 - 8sqrt(Sym(2)))
res2 = solve(eq)
res2[1].evalf() |> println
res2[2].evalf() |> println
res2[3].evalf() |> println
0.45518044897151 + 0.e-21*I
0.275876552559016 + 0.e-24*I
0.683156560842569 - 0.e-22*I
f(t) = 8t^3 - 8sqrt(Sym(2))t^2 +5t - (12 - 8sqrt(Sym(2))) を図にすると以下のようになる。
using Plots
pyplot(size=(500, 500), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(eq, xlims=(0.25, 0.7))
hline!([0])
vline!([0.275876552559016, 0.45518044897151, 0.683156560842569])
