算額（その611）

算額（その611）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/10

等脚台形内に，上円，中円，下円，甲円，乙円と斜線を入れる。中円の直径が 4 寸，子（中円の直径から甲円と乙円の直径の和を引いた長さ）が 1 寸のとき，台形の高さはいかほどか。

台形の下底，上底，高さを a, b, h とする。また，斜線と斜辺の交点座標を (x3, (a - x3)*h/(a - b)), (x4, (a - x4)*h/(a - b)) とする。
上円の半径と中心座標を r1, (0, h - r1)
中円の半径と中心座標を r2, (0, 2r3 + r2)
下円の半径と中心座標を r3, (0, r3)
甲円の半径と中心座標を r4, (0, h - 2r1 - r4)
乙円の半径と中心座標を r5, (0, 2r3 + r5)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, h::positive, x3::positive, x4::positive,
     r1::positive, r2::positive, r3::positive, r4::positive, r5::positive, 子::positive
eq1 = distance(a, 0, b, h, 0, h - r1) - r1^2
eq2 = distance(a, 0, b, h, 0, 2r3 + r2) - r2^2
eq3 = distance(a, 0, b, h, 0, r3) - r3^2
eq4 = distance(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b), 0, h - r1) - r1^2
eq5 = distance(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b), 0, r3) - r3^2
eq6 = distance(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b), 0, h - 2r1 - r4) - r4^2
eq7 = distance(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b), 0, 2r3 + r5) - r5^2
eq8 = 2r2 - 2r4 - 2r5 - 子
eq9 = 2r1 + 2r2 + 2r3 - h;

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 v, r.f_converged
end;

function H(u)
   (a, b, h, x3, x4, r1, r3, r4, r5) = u
   return [
       h^2*(a*r1 + b*h - b*r1)^2/(a^2 - 2*a*b + b^2 + h^2)^2 - r1^2 + (h - h*(a^2 - a*b + h^2 - h*r1)/(a^2 - 2*a*b + b^2 + h^2) - r1)^2,  # eq1
       h^2*(a*h - a*r2 - 2*a*r3 + b*r2 + 2*b*r3)^2/(a^2 - 2*a*b + b^2 + h^2)^2 - r2^2 + (-h*(a^2 - a*b + h*r2 + 2*h*r3)/(a^2 - 2*a*b + b^2 + h^2) + r2 + 2*r3)^2,  # eq2
       h^2*(a*h - a*r3 + b*r3)^2/(a^2 - 2*a*b + b^2 + h^2)^2 - r3^2 + (-h*(a^2 - a*b + h*r3)/(a^2 - 2*a*b + b^2 + h^2) + r3)^2,  # eq3
       h^2*(a*r1*x3^2 - a*r1*x4^2 + b*h*x3^2 - b*h*x4^2 - b*r1*x3^2 + b*r1*x4^2 - 2*h*x3^2*x4 + 2*h*x3*x4^2)^2/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2)^2 - r1^2 + (h - h*(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - a*b*x3^2 - 2*a*b*x3*x4 - a*b*x4^2 - 2*a*x3^2*x4 - 2*a*x3*x4^2 + 2*b*x3^2*x4 + 2*b*x3*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2 - h*r1*x3^2 + 2*h*r1*x3*x4 - h*r1*x4^2)/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2) - r1)^2,  # eq4
       h^2*(a*h*x3^2 - a*h*x4^2 - a*r3*x3^2 + a*r3*x4^2 + b*r3*x3^2 - b*r3*x4^2 - 2*h*x3^2*x4 + 2*h*x3*x4^2)^2/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2)^2 - r3^2 + (-h*(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - a*b*x3^2 - 2*a*b*x3*x4 - a*b*x4^2 - 2*a*x3^2*x4 - 2*a*x3*x4^2 + 2*b*x3^2*x4 + 2*b*x3*x4^2 + h*r3*x3^2 - 2*h*r3*x3*x4 + h*r3*x4^2)/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2) + r3)^2,  # eq5
       h^2*(2*a*r1*x3^2 - 2*a*r1*x4^2 + a*r4*x3^2 - a*r4*x4^2 + b*h*x3^2 - b*h*x4^2 - 2*b*r1*x3^2 + 2*b*r1*x4^2 - b*r4*x3^2 + b*r4*x4^2 - 2*h*x3^2*x4 + 2*h*x3*x4^2)^2/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2)^2 - r4^2 + (h - h*(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - a*b*x3^2 - 2*a*b*x3*x4 - a*b*x4^2 - 2*a*x3^2*x4 - 2*a*x3*x4^2 + 2*b*x3^2*x4 + 2*b*x3*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2 - 2*h*r1*x3^2 + 4*h*r1*x3*x4 - 2*h*r1*x4^2 - h*r4*x3^2 + 2*h*r4*x3*x4 - h*r4*x4^2)/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2) - 2*r1 - r4)^2,  # eq6
       h^2*(a*h*x3^2 - a*h*x4^2 - 2*a*r3*x3^2 + 2*a*r3*x4^2 - a*r5*x3^2 + a*r5*x4^2 + 2*b*r3*x3^2 - 2*b*r3*x4^2 + b*r5*x3^2 - b*r5*x4^2 - 2*h*x3^2*x4 + 2*h*x3*x4^2)^2/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2)^2 - r5^2 + (-h*(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - a*b*x3^2 - 2*a*b*x3*x4 - a*b*x4^2 - 2*a*x3^2*x4 - 2*a*x3*x4^2 + 2*b*x3^2*x4 + 2*b*x3*x4^2 + 2*h*r3*x3^2 - 4*h*r3*x3*x4 + 2*h*r3*x4^2 + h*r5*x3^2 - 2*h*r5*x3*x4 + h*r5*x4^2)/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2) + 2*r3 + r5)^2,  # eq7
       2*r2 - 2*r4 - 2*r5 - 子,  # eq8
       -h + 2*r1 + 2*r2 + 2*r3,  # eq9
   ]
end;

r2 = 4//2
子 = 1
iniv = BigFloat[8.5, 0.5, 16.0, 4.5, 0.9, 0.8, 5.2, 0.2, 1.3]
res = nls(H, ini=iniv)

   (BigFloat[8.472135954999579392818347337462552470881319131313934445299371062594802345966918, 0.472135954999579392818347337462552470881231767596567280477615683449053318730822, 16.00000000000000000000000000000000000000005685400937294517590029875627505751969, 4.47213595499957939281834733746255247088125870170645007332837022256876548756112, 0.8944271909999158785636694674925104941762430878938323936182682049502660273589708, 0.7639320225002103035908263312687237645593763638467272007538116532706423285987476, 5.236067977499789696409173668731276235440652063157959271834138496107495199876993, 0.1909830056250525758977065828171809411398420533761430006087923002379663556053446, 1.309016994374947424102293417182819058860157946623856999391207699762033644413767], true)

中円の直径が 4 寸，子が 1 寸のとき，等脚台形の高さは 16 寸である。

その他のパラメータは以下の通りである。
r2 = 2;  子 = 1;  a = 8.47214;  b = 0.472136;  h = 16;  x3 = 4.47214;  x4 = 0.894427;  r1 = 0.763932;  r3 = 5.23607;  r4 = 0.190983;  r5 1.30902

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, h, x3, x4, r1, r3, r4, r5) = res[1]
   @printf("高さ = %g;  r2 = %g;  子 = %g;  a = %g;  b = %g;  h = %g;  x3 = %g;  x4 = %g;  r1 = %g;  r3 = %g;  r4 = %g;  r5 %g\n",
       h, r2, 子, a, b, h, x3, x4, r1, r3, r4, r5)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:gray60, lw=0.5)
   circle(0, h - r1, r1, :blue)
   circle(0, 2r3 + r2, r2, :magenta)
   circle(0, r3, r3, :green)
   circle(0, h - 2r1 - r4, r4, :brown)
   circle(0, 2r3 + r5, r5, :red)
   segment(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b))
   segment(-x3, (a - x3)*h/(a - b), x4, (a - x4)*h/(a - b))
   segment(0, 2r3 + 2r5, 0, h - 2r1 - 2r4, lw=3)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, h - r1, "上円:r1,(0,h-r1)    ", :blue, :right, :vcenter)
       point(0, h - 2r1 - r4, "     甲円:r4,(0,h-2r1-r4)", :brown, :left, :vcenter)
       point(0, 2r3 + r2, "          中円:r2,(0,2r3+r2)", :magenta, :left, :vcenter)
       point(0, 2r3 + r5, "          乙円:r5,(0,2r3+r5)", :red, :left, :vcenter)
       point(0, r3, " 下円:r3,(0,r3)", :green, :left, :vcenter)
       point(x3, (a - x3)*h/(a - b), " x3", :black, :left, :vcenter)
       point(x4, (a - x4)*h/(a - b), " x4", :black, :left, :vcenter)
       point(a, 0, "a", :black, :left, :bottom, delta=delta/2)
       point(b, h, "(b,h)", :black, :left, :bottom, delta=delta/2)
       point(0, 2r3 + 1.5r2, "子  ", :black, :right, :vcenter, mark=false)
   end
end;