算額（その352）

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