算額（その1119）

算額（その1119）

十 岩手県胆沢町若柳市野々 個人宅 安政2年(1855)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円3個，正方形，斜線2本

正方形の中に斜線 2 本を引き，甲円 2 個，乙円 1 個を容れる。乙円の直径が 1 寸のとき，短い方の斜線の長さはいかほどか。

正方形の一辺の長さを a，斜線と正方形の斜線との交点を (0, b), (c, a)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (a - r2, y2)
とおき，以下の連立方程式を解く。
短い方の斜線の長さは sqrt((a - c)^2  + a^2) である。

include("julia-source.txt")

using SymPy
@syms a::positive, b::positive, c::positive,
     r1::positive, x1::positive,
     r2::positive, y2::positive
eq1 = dist2(0, b, a, a, x1, r1, r1)
eq2 = dist2(0, b, a, a, r1, a - r1, r1)
eq3 = dist2(0, b, a, a, a - r2, y2, r2)
eq4 = dist2(c, a, a, 0, x1, r1, r1)
eq5 = dist2(c, a, a, 0, r1, a - r1, r1)
eq6 = dist2(c, a, a, 0, a - r2, y2, r2);

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, x1, y2) = u
   return [
       a^2*b^2 - 2*a^2*b*r1 + 2*a^2*b*x1 - a^2*r1^2 - 2*a^2*r1*x1 + a^2*x1^2 - 2*a*b^2*x1 + 2*a*b*r1^2 + 2*a*b*r1*x1 - 2*a*b*x1^2 - b^2*r1^2 + b^2*x1^2,  # eq1
       a*(a^3 - 2*a^2*b - 4*a^2*r1 + a*b^2 + 6*a*b*r1 + 2*a*r1^2 - 2*b^2*r1 - 2*b*r1^2),  # eq2
       a*(a^3 - 2*a^2*r2 - 2*a^2*y2 + 2*a*b*r2 - a*r2^2 + 2*a*r2*y2 + a*y2^2 - 2*b*r2*y2),  # eq3
       a*(a^3 - 2*a^2*r1 - 2*a^2*x1 + 2*a*c*r1 - a*r1^2 + 2*a*r1*x1 + a*x1^2 - 2*c*r1*x1),  # eq4
       a*(a*c^2 - 2*a*r1^2 - 2*c^2*r1 + 2*c*r1^2),  # eq5
       -a^2*r2^2 - 2*a^2*r2*y2 + a^2*y2^2 + 2*a*c*r2^2 + 2*a*c*r2*y2 - 2*a*c*y2^2 - c^2*r2^2 + c^2*y2^2,  # eq6
   ]
end;

r2 = 1/2
iniv = BigFloat[2.9, 0.5, 1.2, 0.8, 1.6, 1.9]
res = nls(H, ini=iniv)

   ([2.927373520970951, 0.4800319114286444, 1.2284447668135445, 0.7795448985438439, 1.5736412260872523, 1.8576489200442854], true)

乙円の直径が 1 寸のとき，小斜は 3.3846528098144275 寸である。

「答」では 2寸3分6厘となっているが，そんな長さはどこにもない。

sqrt(2.927373520970951^2 + (2.927373520970951 - 1.2284447668135445)^2)

   3.3846528098144275

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

   r2 = 0.5;  a = 2.92737;  b = 0.480032;  c = 1.22844;  r1 = 0.779545;  x1 = 1.57364;  y2 = 1.85765

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, b, c, r1, x1, y2) = res[1]
   小斜 = sqrt(a^2 + (a - c)^2)
   @printf("乙円の直径が %g のとき，小斜の長さは %g である。\n", 2r2, 小斜)
   @printf("r2 = %g;  a = %g;  b = %g;  c = %g;  r1 = %g;  x1 = %g;  y2 = %g\n", r2, a, b, c, r1, x1, y2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
   circle(x1, r1, r1)
   circle(r1, a - r1, r1)
   circle(a - r2, y2, r2, :blue)
   segment(0, b, a, a, :magenta)
   segment(c, a, a, 0, :magenta) 
   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(0, b, "b ", :magenta, :right, :vcenter)
       point(0, a, "a ", :magenta, :right, :vcenter)
       point(a, 0, " a", :magenta, :center, delta=-delta/2)
       point(c, a, "(c,a)", :magenta, :center, :bottom, delta=delta/2)
       point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(r1, a - r1, "甲円:r1,(r1,a-r1)", :red, :center, delta=-delta/2)
       point(a - r2, y2, "乙円:r2,(a-r2,y2)", :blue, :center, delta=-delta/2)
       plot!(xlims=(-6delta, a+3delta), ylims=(-4delta, a+3delta))
   end
end;