算額（その1167）

算額（その1167）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

長方形内に交差する 2 個の楕円を容れる。楕円の長径が 4 寸，短径が 2.4 寸，長方形の短辺が 3 寸のとき，甲斜の長さはいかほどか。
注：甲斜とは，原点対称の 2 個の楕円の交点間の距離の長い方である。

図に示す x1, y1, x2, y2, x01, y01, x02, y02, x03, y03 を未知数として以下の連立方程式を立て，それを解く。

甲斜の長さは 2sqrt(x03^2 + y03^2) で求める。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, side_a::positive, side_b::positive, K::positive,
     x1::positive, y1::positive, x2::positive, y2::negative,
     x01::positive, y01::positive, 
     x02::positive, y02::positive,
     x03::positive, y03::positive
side_b = sqrt((x2 - x1)^2 + (y1 - y2)^2)
side_a = sqrt((x2 + x1)^2 + (y2 + y1)^2)
eq1 = x01^2/a^2 + y01^2/b^2 - 1
eq2 = -b^2*x01/(a^2*y01) - (y2 - y1)/(x2 - x1)
eq7 = (y2 - y01)/(x2 - x01)  - (y2 - y1)/(x2 - x1)
eq3 = x02^2/a^2 + y02^2/b^2 - 1
eq4 = -b^2*x02/(a^2*y02) - (y2 + y1)/(x2 + x1)
eq8 = (y2 - y02)/(x2 - x02) - (y2 + y1)/(x2 + x1)
eq5 = (4a^2 + 4b^2) - (side_a^2 + side_b^2)  # formula-89
eq6 = (y2 - y1)/(x2 - x1) * (y2 + y1)/(x2 + x1) + 1
eq9 = side_b - K
eq10 = y03/x03 - (y2 + y1)/(x2 + x1)
eq11 = x03^2/a^2 + y03^2/b^2 - 1;

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)
   (x1, y1, x2, y2, x01, y01, x02, y02, x03, y03) = u
   return [
       -1 + y01^2/b^2 + x01^2/a^2,  # eq1
       -(-y1 + y2)/(-x1 + x2) - b^2*x01/(a^2*y01),  # eq2
       -1 + y02^2/b^2 + x02^2/a^2,  # eq3
       -(y1 + y2)/(x1 + x2) - b^2*x02/(a^2*y02),  # eq4
       4*a^2 + 4*b^2 - (-x1 + x2)^2 - (x1 + x2)^2 - (y1 - y2)^2 - (y1 + y2)^2,  # eq5
       1 + (-y1 + y2)*(y1 + y2)/((-x1 + x2)*(x1 + x2)),  # eq6
       -(-y1 + y2)/(-x1 + x2) + (-y01 + y2)/(-x01 + x2),  # eq7
       -K + sqrt((-x1 + x2)^2 + (y1 - y2)^2),  # eq9
       -(y1 + y2)/(x1 + x2) + y03/x03,  # eq10
       -1 + y03^2/b^2 + x03^2/a^2,  # eq11
   ]
end;

(a, b, K) = (4/2, 2.4/2, 3)
iniv = BigFloat[0.614, 2.254, 2.315, -0.228, 1.789, 0.539, 1.505, -0.788, 1.315, 0.91]
res = nls(H, ini=iniv)

   ([0.6329571688388319, 2.244853051407938, 2.3204571688388325, -0.2355388027151154, 1.8516704311458714, 0.45351293387424213, 1.5000000000000002, -0.793725393319377, 1.322875655532295, 0.9000000000000001], true)

(x1, y1, x2, y2, x01, y01, x02, y02, x03, y03) = res[1]
甲斜 = 2sqrt(x03^2 + y03^2)

   3.1999999999999993

甲斜は 3.2 寸である。

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

  長辺 = 3.57211;  短辺 = 3;  x1 = 0.632957;  y1 = 2.24485;  x2 = 2.32046;  y2 = -0.235539;  x01 = 1.85167;  y01 = 0.453513;  x02 = 1.5;  y02 = -0.793725;  x03 = 1.32288;  y03 = 0.9

function draw(p, q, a, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (4/2, 2.4/2)
   (x1, y1) = (0.614, 2.254)
   (x2, y2) = (2.315, -0.228)
   (x01, y01) = (1.789, 0.539)
   (x02, y02) = (1.505, -0.788)
   (x1, y1, x2, y2, x01, y01, x02, y02, x03, y03) = res[1]
   side_b = sqrt((x2 - x1)^2 + (y1 - y2)^2)
   side_a = sqrt((x2 + x1)^2 + (y2 + y1)^2)
   @printf("甲斜 = %g;  長辺 = %g;  短辺 = %g\n", 2sqrt(x03^2 + y03^2), side_a, side_b)
   @printf("長辺 = %g;  短辺 = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g;  x01 = %g;  y01 = %g;  x02 = %g;  y02 = %g;  x03 = %g;  y03 = %g\n", side_a, side_b, x1, y1, x2, y2, x01, y01, x02, y02, x03, y03)
   plot([x1, -x2, -x1, x2, x1], [y1, -y2, -y1, y2, y1], color=:blue, lw=0.5)
   ellipse(0, 0, a, b, color=:red)
   ellipse(0, 0, a, b, φ= 2atand(y03/x03), color=:red)
   segment(x03, y03, -x03, -y03, :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(x1, y1, " (x1,y1)", :blue, :left, :vcenter)
       point(x2, y2, " (x2,y2)", :blue, :left, :vcenter)
       point(-x2, -y2, " (-x2,-y2)", :blue, :left, :vcenter)
       point(-x1, -y1, " (-x1,-y1)", :blue, :left, :vcenter)
       point(x01, y01, " (x01,y01)", :red, :left, :vcenter)
       point(x02, y02, " (x02,y02)", :red, :left, :vcenter)
       point(x03, y03, " (x03,y03)", :green, :left, :vcenter)
       point(-x03, -y03, " (-x03,-y03)", :green, :left, :vcenter)
       xlims!(-x2 - 5delta, x2+10delta)
   end
end;