算額（その2003）

算額（その2003）

(補) 滋賀県栗東町字御園 春日神社 文政7年(1824)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円4個，長方形，斜線2本

長方形の中に，対角線 1 本と斜線 1 本を引き，区分された領域に等円を 4 個容れる。長方形の長辺，短辺がそれぞれ 72寸，65 寸のとき等円の直径はいかほどか。

長方形の長辺，短辺を 2a, 2b
斜線と長方形の短辺の交点座標を (a, c), (-a, -c)
等円の半径と中心座標を r, (x0, b - r), (-x0, r - b), (a - r, y0), (-y0, r - a)
とおき，以下の連立方程式を解く。
注：図形は，x 軸，y 軸で対称ではない。（等円の中心は x 軸，y 軸上にはない）

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, c::positive, r::positive, x0, y0
eq1 = dist2(-a, -c, a, c, x0, b - r, r)
eq2 = dist2(-a, -c, a, c, a - r, y0, r)/a
eq3 = dist2(-a, b, a, -b, x0, b - r, r)/b
eq4 = dist2(-a, b, a, -b, a - r, y0, r)/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)
   (r, c, x0, y0) = u
   return [
       a^2*b^2 - 2*a^2*b*r - 2*a*b*c*x0 + 2*a*c*r*x0 - c^2*r^2 + c^2*x0^2,  # eq1
       a*c^2 - 2*a*c*y0 - a*r^2 + a*y0^2 - 2*c^2*r + 2*c*r*y0,  # eq2
       a^2*b - 2*a^2*r + 2*a*b*x0 - 2*a*r*x0 - b*r^2 + b*x0^2,  # eq3
       a*b^2 + 2*a*b*y0 - a*r^2 + a*y0^2 - 2*b^2*r - 2*b*r*y0,  # eq4
   ]
end;

a = 72//2
b = 65//2

iniv = BigFloat[14, 30, 0.1, 0.1]
res = nls(H, ini=iniv)

   ([14.098029747876524, 28.125, 0.6548773444789615, -0.7794330672778218], true)

等円の直径は 2*14.098029747876524 = 28.196059495753047 である。

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, c, x0, y0) = res[1]
   @printf("a = %g;  b = %g;  r = %g;  c = %g;  x0 = %g;  y0 = %g\n", a, b, r, c, x0, y0)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:green, lw=0.5)
   circle(a - r, y0, r, :blue)
   circle(r - a, -y0, r, :blue)
   circle(x0, b - r, r, :blue)
   circle(-x0, r - b, r, :blue)
   segment(-a, -c, a, c, :red)
   segment(-a, b, a, -b, :red)
   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/2)
       point(0, b, "b", :green, :center, :bottom, delta=delta/2)
       point(a - r, y0, "(a-r,y0)", :blue, :center, delta=-delta/2)
       point(r - a, -y0, "(r-a,-y0)", :blue, :center, delta=-delta/2)
       point(x0, b - r, "(x0,b-r)", :blue, :center, delta=-delta/2)
       point(-x0, r - b, "(-x0,r-b)", :blue, :center, delta=-delta/2)
   end
end;

draw(72/2, 65/2, true)