算額（その2049）

算額（その2049）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円6個

天円 3 個，地円 2 個，人円 1 個が互いに接し合っている。天円，人円の直径を与えたときに地円の直径を求めるすべを述べよ。

天円の半径と中心座標を r1, (x1, y1)
地円の半径と中心座標を r2, (r2, 0)
人円の半径と中心座標を r3, (0, y3)
とおき，以下の連立方程式を解く。
SymPy では連立方程式を一度に解くことはできないので，順に解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, x1::positive, y1::positive,
     r2::positive, r3::positive, y3::positive
eq1 = x1^2 + (y3 + r3 + r1 - y1)^2 - 4r1^2 |> expand
eq2 = (x1 - r2)^2 + y1^2 - (r1 + r2)^2 |> expand
eq3 = x1^2 + (y1 - y3)^2 - (r1 + r3)^2 |> expand
eq4 = r2^2 + y3^2 - (r2 + r3)^2 |> expand;

まず，eq1, eq3 から，x1, x2 を求める。

res = solve([eq1, eq3], (x1, y1))[1]  # 1 of 1

   (2*r1*sqrt(r3)*sqrt(2*r1 + r3)/(r1 + r3), (-r1^2 + 2*r1*r3 + r1*y3 + r3^2 + r3*y3)/(r1 + r3))

eq4 に x1, y1 を代入し y3 を求める。

eq14 = eq4(x1 => res[1], y1 => res[2])
ans_y3 = solve(eq14, y3)[1]  # 1 of 1
ans_y3 |> println

   sqrt(r3)*sqrt(2*r2 + r3)

eq2 に x1, y1, y3 を代入し r2 を求める。

eq12 = eq2(x1 => res[1], y1 => res[2], y3 => ans_y3) |> simplify |> numerator;
ans_r2 = solve(eq12, r2)[2]  # 2 of 2
ans_r2 |> println

   2*r1*(r1^5*r3 + r1^4*r3^(3/2)*sqrt(2*r1 + r3) + r1^4*r3^2 + 4*r1^3*r3^(5/2)*sqrt(2*r1 + r3) + 6*r1^2*r3^(7/2)*sqrt(2*r1 + r3) + 2*r1^2*r3^4 + 4*r1*r3^(9/2)*sqrt(2*r1 + r3) + 3*r1*r3^5 + r3^(11/2)*sqrt(2*r1 + r3) + r3^6 + sqrt(r1^10*r3^2 + 2*r1^9*r3^(5/2)*sqrt(2*r1 + r3) + 2*r1^9*r3^3 + 2*r1^8*r3^(7/2)*sqrt(2*r1 + r3) - 7*r1^8*r3^4 - 16*r1^7*r3^(9/2)*sqrt(2*r1 + r3) - 24*r1^7*r3^5 - 32*r1^6*r3^(11/2)*sqrt(2*r1 + r3) - 10*r1^6*r3^6 + 12*r1^5*r3^(13/2)*sqrt(2*r1 + r3) + 52*r1^5*r3^7 + 92*r1^4*r3^(15/2)*sqrt(2*r1 + r3) + 102*r1^4*r3^8 + 112*r1^3*r3^(17/2)*sqrt(2*r1 + r3) + 88*r1^3*r3^9 + 64*r1^2*r3^(19/2)*sqrt(2*r1 + r3) + 41*r1^2*r3^10 + 18*r1*r3^(21/2)*sqrt(2*r1 + r3) + 10*r1*r3^11 + 2*r3^(23/2)*sqrt(2*r1 + r3) + r3^12))/(r1^6 + 4*r1^5*sqrt(r3)*sqrt(2*r1 + r3) + 10*r1^5*r3 + 8*r1^4*r3^(3/2)*sqrt(2*r1 + r3) + 19*r1^4*r3^2 + 12*r1^3*r3^3 - 8*r1^2*r3^(7/2)*sqrt(2*r1 + r3) + 3*r1^2*r3^4 - 4*r1*r3^(9/2)*sqrt(2*r1 + r3) + 2*r1*r3^5 + r3^6)

式が長く，SymPy では簡約化できないが，天円，人円の半径を与えれば，地円の半径が求まる。
天円，人円の直径が 15 寸，7 寸のとき，地円の直径は 5.72928861346286 寸である。

2ans_r2(r1 => 15/2, r3 => 7/2) |> println

   5.72928861346286

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

  r1 = 7.5;  r3 = 3.5;  r2 = 2.86464;  y3 = 5.68353;  x1 = 10.9728;  y1 = 6.45626

function draw(r1, r3, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 2*r1*(r1^5*r3 + r1^4*r3^(3/2)*sqrt(2*r1 + r3) + r1^4*r3^2 + 4*r1^3*r3^(5/2)*sqrt(2*r1 + r3) + 6*r1^2*r3^(7/2)*sqrt(2*r1 + r3) + 2*r1^2*r3^4 + 4*r1*r3^(9/2)*sqrt(2*r1 + r3) + 3*r1*r3^5 + r3^(11/2)*sqrt(2*r1 + r3) + r3^6 + sqrt(r1^10*r3^2 + 2*r1^9*r3^(5/2)*sqrt(2*r1 + r3) + 2*r1^9*r3^3 + 2*r1^8*r3^(7/2)*sqrt(2*r1 + r3) - 7*r1^8*r3^4 - 16*r1^7*r3^(9/2)*sqrt(2*r1 + r3) - 24*r1^7*r3^5 - 32*r1^6*r3^(11/2)*sqrt(2*r1 + r3) - 10*r1^6*r3^6 + 12*r1^5*r3^(13/2)*sqrt(2*r1 + r3) + 52*r1^5*r3^7 + 92*r1^4*r3^(15/2)*sqrt(2*r1 + r3) + 102*r1^4*r3^8 + 112*r1^3*r3^(17/2)*sqrt(2*r1 + r3) + 88*r1^3*r3^9 + 64*r1^2*r3^(19/2)*sqrt(2*r1 + r3) + 41*r1^2*r3^10 + 18*r1*r3^(21/2)*sqrt(2*r1 + r3) + 10*r1*r3^11 + 2*r3^(23/2)*sqrt(2*r1 + r3) + r3^12))/(r1^6 + 4*r1^5*sqrt(r3)*sqrt(2*r1 + r3) + 10*r1^5*r3 + 8*r1^4*r3^(3/2)*sqrt(2*r1 + r3) + 19*r1^4*r3^2 + 12*r1^3*r3^3 - 8*r1^2*r3^(7/2)*sqrt(2*r1 + r3) + 3*r1^2*r3^4 - 4*r1*r3^(9/2)*sqrt(2*r1 + r3) + 2*r1*r3^5 + r3^6)
   y3 = sqrt(r3)*sqrt(2*r2 + r3)
   (x1, y1) = (2*r1*sqrt(r3)*sqrt(2*r1 + r3)/(r1 + r3), (-r1^2 + 2*r1*r3 + r1*y3 + r3^2 + r3*y3)/(r1 + r3))
   @printf("天円，人円の直径が %g，%g のとき，地円の直径は %g である。\n", 2r1, 2r3, 2r2)
   @printf("r1 = %g;  r3 = %g;  r2 = %g;  y3 = %g;  x1 = %g;  y1 = %g\n", r1, r3, r2, y3, x1, y1)
   plot()
   circle(0, y3 + r3 + r1, r1)
   circle2(x1, y1, r1)
   circle2(r2, 0, r2, :blue)
   circle(0, 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(0, y3 + r3 + r1, "天円:r1,(0,y3+r3+r1)", :red, :center, delta=-delta/2)
       point(x1, y1, "天円:r1,(x1,y1)", :red, :center, delta=-delta/2)
       point(r2, 0, "地円:r2\n(r2,0)", :blue, :center, delta=-delta/2)
       point(0, y3, "人円:r3\n(0,y3)", :green, :center, delta=-delta/2)
   end
end;

draw(15/2, 7/2, true)

算額（その2048）

算額（その2048）

四十九 群馬県安中市板鼻 鷹巣神社 文政11年(1828)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円2個，楕円2個

等楕円と等円が 2 個，互いに接し合っている。楕円の長径，短径が与えられたとき，等円の直径を求める術を述べよ。

楕円の長半径，短半径，中心座標を a, b, (0, b)
等円の半径と中心座標を r, (r, 0)
共通接点の座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r::positive,
     x0::positive, y0::positive
eq1 = x0^2/a^2 + (y0 - b)^2/b^2 - 1
eq2 = (x0 - r)^2 + y0^2 - r^2
eq3 = -b^2*x0/(a^2*(y0 - b)) + (x0 - r)/y0
res = solve([eq1, eq2, eq3], (r, x0, y0))[2];

等円の半径は，以下の式のようになる。計算の順序を考えると見た目ほどは複雑ではないといえるかもしれない。

res[1]

なお，当然とも言えるが，a, b は勝手にどのような値でも取れるわけではない。
また，算額では想定していない，円と楕円の接点が別の場所になる解もある。

function draw(a, b, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x0, y0) = (sqrt(2)*sqrt((a^2 + 3*b^2 - sqrt(a^4 - 10*a^2*b^2 + 9*b^4))/(a^2 - b^2))*(3*a^2 - 3*b^2 + sqrt(a^4 - 10*a^2*b^2 + 9*b^4))/(8*a), sqrt(2)*a*sqrt((a^2 + 3*b^2 - sqrt(a^4 - 10*a^2*b^2 + 9*b^4))/(a^2 - b^2))/2, b*(3*a^2 - 3*b^2 + sqrt(a^4 - 10*a^2*b^2 + 9*b^4))/(2*(a^2 - b^2)))
   @printf("a = %g;  b = %g;  r = %g;  x0 = %g;  y0 = %g\n", a, b, r, x0, y0)
   plot()
   #ellipse(0, b, a, b, color=:red)
   #ellipse(0, -b, a, b, color=:red)
   circle2(r, 0, r, :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(r, 0, " 等円:r,(r,0)", :blue, :left, :vcenter)
       point(0, b, "   等楕円:a,b,(0,b)", :red, :left, :vcenter)
       point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta/2)
   end
   ellipse(0, b, a, b, color=:red)
   ellipse(0, -b, a, b, color=:red)
end;

draw(1, 0.25, true)