算額（その2002）

算額（その2002）

(7) 滋賀県 錦織良夫氏宅　p.114
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円4個，楕円，算法助術公式96

楕円の中に 4 個の等円を容れる。楕円の長径，短径が与えられたとき，等円の直径を求めよ。

1.　算法助術の公式96

楕円の長半径，短半径を a, b，等円の半径を r とおく。

算法助術では，楕円の長軸 p = 2a，短軸 q = 2b，円の直径 d = 2r を使っているので，変換する。

include("julia-source.txt");

using SymPy

@syms p, q, d, a, b, r
p = 2a
q = 2b
d = 2r
eq1 = p^2*q^2*(p^2 - q^2)^2 - 4d^2*(p^2 + q^2)*(p^2 - q^2)^2 + 12d^4*(p^2 - q^2)^2 + 4d^4*p^2*q^2 - 12d^6*(p^2 + q^2) + 4d^8;

eq2 = eq1 |> simplify |> factor |> x -> x/256;
eq2 |> println

   a^6*b^2 - 4*a^6*r^2 - 2*a^4*b^4 + 4*a^4*b^2*r^2 + 12*a^4*r^4 + a^2*b^6 + 4*a^2*b^4*r^2 - 20*a^2*b^2*r^4 - 12*a^2*r^6 - 4*b^6*r^2 + 12*b^4*r^4 - 12*b^2*r^6 + 4*r^8

r の 8次式なので r について解くと 8 個の解が得られるが，其のうちの　6 番目が適解である。
式は長いので省略する。

res = solve(eq2, r)[6];

長半径，短半径が 10, 4 のとき，等円の半径は 1.95789337209971 である。

res(a => 10, b => 4).evalf() |>  println

   1.95789337209971

2.　連立方程式を解く

楕円の長半径，短半径を a, b
等円の半径と中心座標を r, (r, r)
楕円と等円の接点座標を (x0, y0)
とおき，以下の連立方程式を解き数値解を求める。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r::positive,
     x0::positive, y0::positive
@syms a, b, x0, y0, r
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (x0 - r)/(y0 - r)
eq3 = (x0 - r)^2 + (y0 - r)^2 - r^2;

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 Float64.(v), r.f_converged
end;

function H(u)
   (r, x0, y0) = u
   return [
       -1 + y0^2/b^2 + x0^2/a^2,  # eq1
       (-r + x0)/(-r + y0) - b^2*x0/(a^2*y0),  # eq2
       -r^2 + (-r + x0)^2 + (-r + y0)^2,  # eq3
   ]
end;

(a, b) = (10, 4)
iniv = BigFloat[1.943, 2.13, 3.877]
res = nls(H, ini=iniv)

   ([1.9578933720997147, 2.127795322101097, 3.9084009429380644], true)

楕円の長径，短径が 20，8 のとき，等円の直径は 3.9157867441994294 である。

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x0, y0) = res[1]
   @printf("a = %g;  b = %g;  r = %g;  x0 = %g;  y0 = %g\n", a, b, r, x0, y0)
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle4(r, r, r, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(a, 0, " a", :green, :left, :bottom, delta=delta)
       point(0, b, "b", :green, :center, :bottom, delta=delta)
       point(r, r, "等円:r\n(r,r)", :blue, :center, delta=-delta)
       point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta)
   end
end;

draw(10, 4, true)