算額（その1289）

算額（その1289）

百三十五 群馬県安中市鷺宮 咲前神社 明治20年(1887)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，正方形，斜線

正方形の中に右上の頂点を通る斜線を引き，区分された領域に大円 1 個，等円 3 個を容れる。正方形の一辺の長さが 7.2 寸，大円の直径が 3.6 寸，斜線の長さが 9 寸のとき，等円の直径はいかほどか。

正方形の一辺の長さを a
斜線と正方形の下辺との交点を (b, 0)
大円の半径と中心座標を r1, (a - r1, r1)
等円の半径と中心座標を r2, (r2, y21), (x22, y22), (x23, a - r2)
とおき，以下の連立方程式を解く。

条件が 3 つあげられているが，斜線の長さと正方形の一辺の長さと大円の直径の間には以下に示される関係があり，それぞれを別々に，任意に設定できるものではない。
斜 = sqrt(a^2 + (a - b)^2)
a + (a - b) - 斜 = 2r1
a = 7.2，r1 = 3.6/2 のときには，斜 = 9 になるというだけのことである。

ans_b = solve(eq6, b)[1] 
ans_b |> println

   (a^2 - 4*a*r1 + 2*r1^2)/(a - 2*r1)

a = 7.2，r1 = 3.6/2 のとき，b = 1.8, 斜 = 9 になる。

ans_b(a => 7.2, r1 => 3.6/2) |> println  # 斜 => 9)  

   1.80000000000000

sqrt(a^2 + (a - ans_b)^2)(a => 7.2, r1 => 3.6/2) |> println

   9.00000000000000

斜は a, r1 から計算される。

eq6 を b について解くと，b = (a^2 - 4*a*r1 + 2*r1^2)/(a - 2*r1) となる。

それ以外の値を与えると破綻する。

よって，与えられるべき条件はこの 3 つの条件の内の 2 個だけである。a と r1 を与えるのが自然であろう。以下では，a, r1 のみが与えられる一般解を求めることにする。

SymPy の solve() では解析解を求めることができないので，nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive,
     y21::positive, x22::positivr, y22::positive, x23::positive;
eq1 = dist2(b, 0, a, a, r2, y21, r2)
eq2 = dist2(b, 0, a, a, x22, y22, r2)
eq3 = dist2(b, 0, a, a, x23, a - r2, r2)
eq4 = (x22 - r2)^2 + (y22 - y21)^2 -  4r2^2
eq5 = (x23 - x22)^2 + (a - r2 - y22)^2 - 4r2^2
eq6 = a + (a - b) - sqrt(a^2 + (a - b)^2) - 2r1;

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)
   (r2, y21, x22, y22, x23, b) = u
   return [
       a^2*b^2 - 2*a^2*b*r2 + 2*a^2*b*y21 - a^2*r2^2 - 2*a^2*r2*y21 + a^2*y21^2 - 2*a*b^2*y21 + 2*a*b*r2^2 + 2*a*b*r2*y21 - 2*a*b*y21^2 - b^2*r2^2 + b^2*y21^2,  # eq1
       a^2*b^2 - 2*a^2*b*x22 + 2*a^2*b*y22 - 2*a^2*r2^2 + a^2*x22^2 - 2*a^2*x22*y22 + a^2*y22^2 - 2*a*b^2*y22 + 2*a*b*r2^2 + 2*a*b*x22*y22 - 2*a*b*y22^2 - b^2*r2^2 + b^2*y22^2,  # eq2
       a*(a^3 - 2*a^2*r2 - 2*a^2*x23 + 2*a*b*r2 - a*r2^2 + 2*a*r2*x23 + a*x23^2 - 2*b*r2*x23),  # eq3
       -4*r2^2 + (-r2 + x22)^2 + (-y21 + y22)^2,  # eq4
       -4*r2^2 + (-x22 + x23)^2 + (a - r2 - y22)^2,  # eq5
       2*a - b - 2*r1 - sqrt(a^2 + (a - b)^2),  # eq6
   ]
end;

(a, r1) = (7.2, 3.6/2, 9)
iniv = BigFloat[1.3, 1.5, 2.9, 3.7, 4.5, 1.8]

res = nls(H, ini=iniv)

   ([1.3333333333333335, 1.6000000000000003, 2.9333333333333336, 3.7333333333333334, 4.533333333333334, 1.7999999999999998], true)

正方形の一辺の長さ，大円の直径，斜がそれぞれ 7.2, 3.6, 9 のとき，等円の直径は 8/3 ≒ 2.66667 である。

すべてのパラメータは以下のとおりである。

   a = 7.2;  r1 = 1.8;  斜 = 9;  r2 = 1.33333;  y21 = 1.6;  x22 = 2.93333;  y22 = 3.73333;  x23 = 4.53333;  b = 1.8

function draw(a, r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, y21, x22, y22, x23, b) = res[1]
   斜 = sqrt(a^2 + (a - b)^2)
   @printf("正方形の一辺の長さ，大円の直径，斜がそれぞれ %g, %g, %g のとき，等円の直径は %g である。\n", a, 2r1, 斜, 2r2)
   @printf("a = %g;  r1 = %g;  斜 = %g;  r2 = %g;  y21 = %g;  x22 = %g;  y22 = %g;  x23 = %g;  b = %g\n",
       a, r1, 斜, r2, y21, x22, y22, x23, b)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, y21, r2, :blue)
   circle(x22, y22, r2, :blue)
   circle(x23, a - r2, r2, :blue)
   segment(b, 0, a, a, :orange)
   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)
   end  
end;

draw(7.2, 3.6/2, true)

以下の図は，正方形の一辺の長さ，大円の直径，斜がそれぞれ 7, 4, 9.66667 のときのもので，等円の直径は 2.29811 である。

    a = 7;  r1 = 2;  斜 = 9.66667;  r2 = 1.14906;  y21 = 2.52264;  x22 = 2.73396;  y22 = 4.18679;  x23 = 4.31887;  b = 0.333333

draw(7, 4/2, true)