算額（その1363）

算額（その1363）

三十三 群馬県佐波郡境町伊与久 雷電神社 文化14年(1817)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円14個，外円，正三角形
#Julia, #SymPy, #算額, #和算

外円の中に正三角形と円 13 個を容れる。甲円の直径が 39.9 寸のとき，丙円の直径はいかほどか。

図形を反時計回りに 90° 回転させたものを考える。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, 0)
乙円の半径と中心座標を r2, (r1 + r2, r2)
丙円の半径と中心座標を r3, (r1 + r3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     r3::positive, y3::positive;
R = 2r1
eq1 = (r1 + r2)^2 + r2^2 - (R - r2)^2
eq2 = (r1 + r3)^2 + y3^2 - (R - r3)^2
eq3 = y3 - r2 - 2sqrt(r2*r3)
res = solve([eq1, eq2, eq3], (r2, r3, y3))[1];

丙円の半径は r1*(-267 - 56*sqrt(21 - 12*sqrt(3)) + 36*sqrt(7 - 4*sqrt(3)) + 190*sqrt(3))/169 となるが，これは二重根号を解消することにより，実にスッキリと簡約化される。

res[2] |> println

   r1*(-267 - 56*sqrt(21 - 12*sqrt(3)) + 36*sqrt(7 - 4*sqrt(3)) + 190*sqrt(3))/169

丙円の半径は，甲円の半径の 3(14√3 - 9)/169 倍である。

res[2] |> sympy.sqrtdenest |> simplify |> println

   3*r1*(-9 + 14*sqrt(3))/169

甲円の直径が 39.9 寸のとき，丙円の直径は 39.9 * 3(14√3 - 9)/169 = 10.800418599549849 寸である。

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

  R = 39.9;  r1 = 19.95;  r2 = 9.25883;  r3 = 5.40021;  y3 = 23.4009

function draw(r1, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, y3) = (r1*(-3 + 2*sqrt(3)), r1*(-267 - 56*sqrt(21 - 12*sqrt(3)) + 36*sqrt(7 - 4*sqrt(3)) + 190*sqrt(3))/169, r1*(-6*sqrt(3) + 12*sqrt(7 - 4*sqrt(3)) + 16*sqrt(21 - 12*sqrt(3)) + 15)/13)
   R = 2r1
   @printf("甲円の直径が %g のとき，丙円の直径は %g である。\n", 2r1, 2r3)
   @printf("R = %g;  r1 = %g;  r2 = %g;  r3 = %g;  y3 = %g\n", R, r1, r2, r3, y3)
   plot([r1, -R, r1, r1], [√3r1, 0, -√3r1, √3r1], color=:brown, lw=0.5)
   circle(0, 0, R, :orange)
   circle(0, 0, r1, :magenta)
   rotate(r1 + r2, r2, r2, :blue)
   rotate(r1 + r2, -r2, r2, :blue)
   rotate(r1 + r3, y3, r3, :green)
   rotate(r1 + r3, -y3, r3, :green)
   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 + r2, r2, "乙円:r2\n(r1+r2,r2)", :blue, :center, delta=-delta)
       point(r1 + r3, y3, "丙円:r3,(r1+r3,y3) ", :green, :right, :vcenter)
       point(r1, 0, "r1 ", :magenta, :right, :bottom, delta=delta/2)
       point(R, 0, " R", :orange, :left, :bottom, delta=delta/2)
       point(r1, √3r1, "(r1,√3r1)", :brown, :left, :vcenter)
   end  
end;

draw(39.9/2, true)

算額（その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)