算額（その1317）

算額（その1317）

百四十七 群馬県甘楽郡妙義町下高田 高太神社 大正12年(1923)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，等脚台形
#Julia, #SymPy, #算額, #和算

等脚台形の中に大円 2 個，小円 2 個を容れる。台形の高さが 20 寸，大円の直径が 14.4 寸のとき，小円の直径はいかほどか。

等脚台形の上底，下底の長さをそれぞれ 2b, 2a 
等脚台形の斜辺を延長してできる二等辺三角形の高さを h0
大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (r2, h - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, h::positive, h0::positive,
     r1::positive, r2::positive, y2::positive;
eq1 = r1/(h0 - r1) - r2/(h0 - (h - r2))
eq2 = b/(h0 - h) - a/h0
eq3 = dist2(0, h0, a, 0, r1, r1, r1)
eq4 = (r1 - r2)^2 + (h - r2 - r1)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3, eq4], (r2, a, b, h0))[2]

   ((-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2), 4*(h^(17/2)*r1^(3/2) - 26*h^(15/2)*r1^(5/2) + 249*h^(13/2)*r1^(7/2) - 1192*h^(11/2)*r1^(9/2) + 3160*h^(9/2)*r1^(11/2) - 4768*h^(7/2)*r1^(13/2) + 3984*h^(5/2)*r1^(15/2) - 1664*h^(3/2)*r1^(17/2) + 256*sqrt(h)*r1^(19/2) - h^8*r1^2 + 26*h^7*r1^3 - 249*h^6*r1^4 + 1192*h^5*r1^5 - 3160*h^4*r1^6 + 4768*h^3*r1^7 - 3984*h^2*r1^8 + 1664*h*r1^9 - 256*r1^10)/(h^9 - 28*h^8*r1 + 301*h^7*r1^2 - 1690*h^6*r1^3 + 5544*h^5*r1^4 - 11088*h^4*r1^5 + 13520*h^3*r1^6 - 9632*h^2*r1^7 + 3584*h*r1^8 - 512*r1^9), 4*(-h^(9/2)*sqrt(r1) + 13*h^(7/2)*r1^(3/2) - 52*h^(5/2)*r1^(5/2) + 68*h^(3/2)*r1^(7/2) - 16*sqrt(h)*r1^(9/2) - h^4*r1 + h^3*r1^2 + 24*h^2*r1^3 - 52*h*r1^4 + 16*r1^5)/(h^4 - 18*h^3*r1 + 84*h^2*r1^2 - 120*h*r1^3 + 32*r1^4), (sqrt(h)*r1*(h - 6*r1) + 2*r1^(3/2)*(-h + 2*r1))/(sqrt(h)*(h - 4*r1)))

res[1] |> println

   (-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2)

小円の半径 r2 は，h, r1 を用いてい以下のように計算できる。

f(h, r1) = (-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2)
f(20, 14.4/2) * 2

   6.399999999999852

台形の高さが 20 寸，大円の直径が 14.4 寸のとき，小円の直径は 6.4 寸である。

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

   h = 20;  r1 = 7.2;  r2 = 3.2;  a = 24.6857;  b = 4.51765;  h0 = 24.48

function draw(h, r1, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, a, b, h0) = ((-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2), 4*(h^(17/2)*r1^(3/2) - 26*h^(15/2)*r1^(5/2) + 249*h^(13/2)*r1^(7/2) - 1192*h^(11/2)*r1^(9/2) + 3160*h^(9/2)*r1^(11/2) - 4768*h^(7/2)*r1^(13/2) + 3984*h^(5/2)*r1^(15/2) - 1664*h^(3/2)*r1^(17/2) + 256*sqrt(h)*r1^(19/2) - h^8*r1^2 + 26*h^7*r1^3 - 249*h^6*r1^4 + 1192*h^5*r1^5 - 3160*h^4*r1^6 + 4768*h^3*r1^7 - 3984*h^2*r1^8 + 1664*h*r1^9 - 256*r1^10)/(h^9 - 28*h^8*r1 + 301*h^7*r1^2 - 1690*h^6*r1^3 + 5544*h^5*r1^4 - 11088*h^4*r1^5 + 13520*h^3*r1^6 - 9632*h^2*r1^7 + 3584*h*r1^8 - 512*r1^9), 4*(-h^(9/2)*sqrt(r1) + 13*h^(7/2)*r1^(3/2) - 52*h^(5/2)*r1^(5/2) + 68*h^(3/2)*r1^(7/2) - 16*sqrt(h)*r1^(9/2) - h^4*r1 + h^3*r1^2 + 24*h^2*r1^3 - 52*h*r1^4 + 16*r1^5)/(h^4 - 18*h^3*r1 + 84*h^2*r1^2 - 120*h*r1^3 + 32*r1^4), (sqrt(h)*r1*(h - 6*r1) + 2*r1^(3/2)*(-h + 2*r1))/(sqrt(h)*(h - 4*r1)))
   @printf("台形の高さが %g，大円の直径が %g のとき，小円の直径は %g である。\n", h, 2r1, 2r2)
   @printf("h = %g;  r1 = %g;  r2 = %g;  a = %g;  b = %g;  h0 = %g\n", h, r1, r2, a, b, h0)
   plot([a, 0, -a, a], [0, h0, 0, 0], color=:green, lw=0.5)
   segment(-b, h, b, h, :green)
   circle2(r1, r1, r1)
   circle2(r2, h - r2, r2, :blue)
   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(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(r2, h - r2, " 小円:r2,(r2,h-r2)", :blue, :left, :vcenter)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(0, h, " h", :green, :left, :bottom, delta=delta/2)
       point(0, h0, " h0", :green, :left, :bottom, delta=delta/2)
       point(b, h, " (b,h)", :green, :left, :vcenter)
   end
end;

draw(20, 14.4/2, true)