算額（その1072）

算額（その1072）

九十六 大船渡市立根 気仙安養寺稲荷堂 文政5年(1822)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

正方形の中に楕円，元円，利円，貞円，享円を容れる。貞円の直径が 1 寸のとき，元円の直径はいかほどか。

注1：利円，享円の条件については何も述べられていない。実際，制限はあるもののどんな値も取りうる。
注2：算額の図（山村の図）では，楕円が楕円らしくないが，それはよくあることかもしれないが，結果に影響を与えるかもしれない。

正方形の一辺の長さを 2a
楕円の長半径，短半径，中心座標を a, b, (0, b); b = a - r1
元円の半径と中心座標を r1, (0, 2a - r1), (r1, a), (r1, a - 2r1)
貞円の半径と中心座標を r2, (a - r3, r3)
楕円と元円の接点座標を (x1, y1)
楕円と貞円の接点座標を (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive,
     r1::positive, x1::positive, y1::positive,
     r2::positive, x2::positive, y2::positive
b = a - r1
eq1 = x2^2/a^2 + (y2 - a + r1)^2/b^2 - 1
eq2 = (x2 - a + r2)^2 + (y2 - r2)^2 - r2^2
eq3 = -b^2*x2/(a^2*(y2 - a + r1)) - (a - r2 - x2)/(y2 - r2)
eq4 = x1^2/a^2 + (y1 - a + r1)^2/b^2 - 1
eq5 = (x1 - r1)^2 + (y1 - a)^2 - r1^2
eq6 = -b^2*x1/(a^2*(y1 - a + r1)) - (r1 - x1)/(y1 - a);

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)
   (a, r1, x1, y1, x2, y2) = u
   return [
       -1 + (-a + r1 + y2)^2/(a - r1)^2 + x2^2/a^2,  # eq1
       -r2^2 + (-r2 + y2)^2 + (-a + r2 + x2)^2,  # eq2
       -(a - r2 - x2)/(-r2 + y2) - x2*(a - r1)^2/(a^2*(-a + r1 + y2)),  # eq3
       -1 + (-a + r1 + y1)^2/(a - r1)^2 + x1^2/a^2,  # eq4
       -r1^2 + (-a + y1)^2 + (-r1 + x1)^2,  # eq5
       -(r1 - x1)/(-a + y1) - x1*(a - r1)^2/(a^2*(-a + r1 + y1)),  # eq6
   ]
end;

r2 = 1/2
iniv = BigFloat[2.1, 0.63, 0.83, 2.6, 1.5, 0.52]
res = nls(H, ini=iniv)

   ([3.5742972321311055, 1.1362300699229417, 1.470601088978576, 4.6602138004791485, 2.755899335723288, 0.8855162507120506], true)

元円の直径は 2*1.1362300699229417 = 2.2724601398458835 である。「答」では 2.243 寸有奇となっており，かなり異なった結論になった。

ちなみに，楕円の長半径，短半径は 3.5743, b = 2.43807 とこれまた，あまりきれいな数ではない。

2*1.1362300699229417

   2.2724601398458835

その他のパラメータは以下のとおりである。

   r2 = 0.5;  a = 3.5743;  b = 2.43807;  r1 = 1.13623;  r2 = 0.5;  x1 = 1.4706;  y1 = 4.66021;  x2 = 2.7559;  y2 = 0.885516

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, r1, x1, y1, x2, y2) = res[1]
   b = a - r1
   @printf("貞円の直径が %g のとき，元円の直径は %g である。\n", 2r2, 2r1)
   @printf("r2 = %g;  a = %g;  b = %g;  r1 = %g;  r2 = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g\n",
       r2, a, b, r1, r2, x1, y1, x2, y2)
   plot([a, a, -a, -a, a], [0, 2a, 2a, 0, 0], color=:green, lw=0.5)
   ellipse(0, a - r1, a, b, color=:blue)
   circle(0, 2a - r1, r1)
   circle2(r1, a, r1)
   circle2(r1, a - 2r1, r1)
   circle2(a - r2, r2, r2, :green)
   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(x2,  y2, "(x2,y2) ", :green, :right, :bottom, delta=2delta, deltax=4delta)
       point(x1,  y1, "(x1,y1)", :green, :left, :bottom, delta=delta/2)
       point(r1, a, "元円:r1,(r1, a)", :red, :center, delta=-delta/2)
       point(r1, a - 2r1, "元円:r1,(r1, a-2r1)", :red, :center, :bottom, delta=delta/2)
       point(0, 2a - r1, "元円:r1,(0,2a-r1)", :red, :center, delta=-delta/2)
       point(a - r2, r2, "貞円:r3,(a-r3,r3)", :green, :center, delta=-delta/2)
       point(0, 2a - 2r1, "2a-2r1", :blue, :center, delta=-delta/2)
       point(0, a - r1, "a-r1", :blue, :center, delta=-delta/2)
   end
end;