算額（その882）

算額（その882）

六四 加須市不動岡 総願寺 慶応2年(1866)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

長方形の中に対角線と斜線を引き，区画された領域に甲円，乙円 2 個ずつを入れる。乙円の直径が 1 寸のとき，長方形の長辺の長さはいかほどか。

長方形の長辺，短辺をそれぞれ a, b
斜線と長辺の交点座標を (c, 0)
甲円の半径と中心座標を r1, (a - r1, y1), (a/2, r1)
乙円の半径と中心座標を r2, (r2, b/2 + r1)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms d, a::positive, b::positive, c::positive,
     r1::positive, y1::positive, r2::positive
eq1 = numerator(apart(dist(a, 0, 0, b, a - r1, y1) - r1^2, d))
eq2 = numerator(apart(dist(a, 0, 0, b, a/2, r1) - r1^2, d))
eq3 = numerator(apart(dist(a, 0, 0, b, r2, b/2 + r2) - r2^2, d))
eq4 = numerator(apart(dist(a, b, c, 0, a - r1, y1) - r1^2, d))
eq5 = numerator(apart(dist(a, b, c, 0, a/2, r1) - r1^2, d));

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, b, c, r1, y1) = u
   return [
       -a^2*r1^2 + a^2*y1^2 - 2*a*b*r1*y1,  # eq1
       a^2*b^2 - 4*a^2*b*r1 - 4*b^2*r1^2,  # eq2
       a^2*b^2 - 4*a^2*b*r2 - 4*a*b^2*r2 + 8*a*b*r2^2,  # eq3
       a^2*b^2 - 2*a^2*b*y1 - a^2*r1^2 + a^2*y1^2 - 2*a*b^2*c - 2*a*b^2*r1 + 4*a*b*c*y1 + 2*a*b*r1*y1 + 2*a*c*r1^2 - 2*a*c*y1^2 + b^2*c^2 + 2*b^2*c*r1 - 2*b*c^2*y1 - 2*b*c*r1*y1 - c^2*r1^2 + c^2*y1^2,  # eq4
       a^2*b^2 - 4*a^2*b*r1 - 4*a*b^2*c + 12*a*b*c*r1 + 4*b^2*c^2 - 4*b^2*r1^2 - 8*b*c^2*r1,  # eq5
   ]
end;

r2 = 1/2
iniv = BigFloat[53, 28, 30, 6, 11] ./5
res = nls(H, ini=iniv)

   ([5.23606797749979, 2.618033988749895, 2.8944271909999157, 0.6180339887498949, 1.0], true)

長方形の長辺の長さは 5.23606797749979 である。

算額では「術曰置五個開平方加三個得直長合問」は √5 + 3 = 5.23606797749979 なので，数値解がこれに一致していることがわかる。

その他のパラメータは以下のとおりである。
r2 = 0.5;  a = 5.23607;  b = 2.61803;  c = 2.89443;  r1 = 0.618034;  y1 = 1

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, b, c, r1, y1) = res[1]
   @printf("乙円の直径が %g のとき，長方形の長辺は %g である\n", 2r2, a)
   @printf("r2 = %g;  a = %g;  b = %g;  c = %g;  r1 = %g;  y1 = %g\n", r2, a, b, c, r1, y1)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:magenta, lw=0.5)
   circle(a - r1, y1, r1)
   circle(a/2, r1, r1)
   circle(r2, b/2 + r2, r2, :green)
   circle(r2, b/2 - r2, r2, :green)
   segment(0, 0, a, b, :blue)
   segment(0, b, a, 0, :blue)
   segment(c, 0, a, b, :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", :magenta, :left, :bottom, delta=delta/2)
       point(c, 0, "  c", :magenta, :left, :bottom, delta=delta/2)
       point(0, b, " b", :magenta, :left, :bottom, delta=delta/2)
       point(a - r1, y1, "甲円:r1\n(a-r1,y1)", :red, :center, delta=-delta/2)
       point(a/2, r1, "甲円:r1\n(a/2,r1)", :red, :center, delta=-delta/2)
       point(r2, b/2 + r2, "乙円:r2\n(r2,b/2+r2)", :green, :center, :bottom, delta=delta/2)
   end
end;