算額（その1362）

算額（その1362）

三十一 群馬県藤岡市立石 立石神社 文化13年(1816)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，等脚台形
#Julia, #SymPy, #算額, #和算

等脚台形の中に全円 1 個，等円 4 個を容れる。全円の直径が与えられたとき，等円の直径を求める術をのべよ。

等脚台形の下底，上底，高さを 2a, 2b, h
全円の半径と中心座標を r1, (0, r1)
等円の半径と中心座標を r2, (x2, r2), (r2, h - r2)
とおき，以下の連立方程式の数値解を求める。
SymPy では，性能的に解析解を求めることができない。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, h::positive,
     r1::positive, r2::positive, x2::positive;
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = r2^2 + (h - r2 - r1)^2 - (r1 + r2)^2
eq3 = dist2(b, h, a, 0, x2, r2, r2)
eq4 = dist2(b, h, a, 0, r2, h - r2, r2)
eq5 = dist2(b, h, a, 0, 0, r1, r1);

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, h, r2, x2) = u
   return [
       x2^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq1
       r2^2 - (r1 + r2)^2 + (h - r1 - r2)^2,  # eq2
       h*(a^2*h - 2*a^2*r2 + 2*a*b*r2 - 2*a*h*x2 + 2*a*r2*x2 - 2*b*r2*x2 - h*r2^2 + h*x2^2),  # eq3
       h*(2*a*b*r2 - 2*a*r2^2 + b^2*h - 2*b^2*r2 - 2*b*h*r2 + 2*b*r2^2),  # eq4
       h*(a^2*h - 2*a^2*r1 + 2*a*b*r1 - h*r1^2),  # eq5
   ]
end;

r1 = 1
iniv = BigFloat[1.5, 0.5, 2.5, 0.28, 1.1]
res = nls(H, ini=iniv)

   ([1.4706834198711607, 0.4706834198711606, 2.5293165801288393, 0.28017604199929, 1.0586331602576788], true)

全円の直径が 2 のとき，等円の直径は 0.56035208399858 である。

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

    r1 = 1;  a = 1.47068;  b = 0.470683;  h = 2.52932;  r2 = 0.280176;  x2 = 1.05863

function draw(r1, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, h, r2, x2) = res[1]
   @printf("全円の直径が %g のとき，等円の直径は %g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  a = %g;  b = %g;  h = %g;  r2 = %g;  x2 = %g\n", r1, a, b, h, r2, x2)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:green, lw=0.5)
   circle(0, r1, r1, :blue)
   circle2(x2, r2, r2)
   circle2(r2, h - r2, r2)
   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, r1, "全円:r1,(0,r1)", :blue, :center, delta=-delta/2)
       point(x2, r2, "等円:r2\n(x2,r2)", :red, :center, delta=-delta/2)
       point(r2, h - r2, "等円:r2\n(r2,h-r2)", :red, :center, :bottom, delta=delta/2)
       point(b, h, "(b,h)", :green, :right, :bottom, delta=delta/2)
       point(a, 0, "a", :green, :left, :bottom, delta=delta/2)
   end  
end;

draw(1, true)