算額（その1128）

算額（その1128）

四十四 岩手県一関市真滝 熊野白山滝神社 文久元年(1861)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円5個，直線

互いに接している大円と中円があり，両者を通る直線と各円の隙間に 3 個の小円を容れる。大円と小円の直径が 10 寸，1 寸のとき，中円の直径はいかほどか。

大円の半径と中心座標を r1, (x1, y1)
中円の半径と中心座標を r2, (0, r2)
小円の半径と中心座標を r3, (0, 2r2 - r3), (x3, 2r2 - 3r3), (x1, y1 + r1 - r3)
とおき，以下の連立方程式を解く。

注：大円と小円の直径が 10 寸，1 寸のとき，y1 = 0 となるが，これは特殊な場合で常に y1 = 0 が成り立つわけではない。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive, x1::positive, y1,
     r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive;
y1 = 2r2 - r1
eq1 = (x1 - x3)^2 + (2r2 - 3r3 - y1)^2 - (r1 + r3)^2
eq2 = x3^2 + (r2 - 3r3)^2 - (r2 + r3)^2
eq3 = x1^2 + (r2 - y1)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r2, x1, x3))[4]  # 4 of 4

   (2*r1*(r1^3*r3 + 2*r1^2*r3^(3/2)*sqrt(r1 - r3) - 4*r1^2*r3^2 - 8*r1*r3^(5/2)*sqrt(r1 - r3) + 4*r1*r3^3 + 8*r3^(7/2)*sqrt(r1 - r3))/(r1^4 - 8*r1^3*r3 + 24*r1^2*r3^2 - 32*r1*r3^3 + 16*r3^4), 2*sqrt(2)*r1*sqrt((r1^4*r3^2 + 4*r1^3*r3^(5/2)*sqrt(r1 - r3) - 16*r1^2*r3^(7/2)*sqrt(r1 - r3) - 16*r1^2*r3^4 + 16*r1*r3^(9/2)*sqrt(r1 - r3) + 32*r1*r3^5 - 16*r3^6)/(r1^4 - 8*r1^3*r3 + 24*r1^2*r3^2 - 32*r1*r3^3 + 16*r3^4))*(r1^4*r3^(5/2)*sqrt(r1 - r3) - r1^4*r3^3 - 8*r1^3*r3^(7/2)*sqrt(r1 - r3) + 8*r1^3*r3^4 + 24*r1^2*r3^(9/2)*sqrt(r1 - r3) - 24*r1^2*r3^5 - 32*r1*r3^(11/2)*sqrt(r1 - r3) + 32*r1*r3^6 + 16*r3^(13/2)*sqrt(r1 - r3) - 16*r3^7)/(r3^3*(r1^5 - 10*r1^4*r3 + 40*r1^3*r3^2 - 80*r1^2*r3^3 + 80*r1*r3^4 - 32*r3^5)), 2*sqrt(2)*sqrt((r1^4*r3^2 + 4*r1^3*r3^(5/2)*sqrt(r1 - r3) - 16*r1^2*r3^(7/2)*sqrt(r1 - r3) - 16*r1^2*r3^4 + 16*r1*r3^(9/2)*sqrt(r1 - r3) + 32*r1*r3^5 - 16*r3^6)/(r1^4 - 8*r1^3*r3 + 24*r1^2*r3^2 - 32*r1*r3^3 + 16*r3^4)))

かなり複雑な式で分母は因数分解できるも場合もあるが(r2, x1)，分子を含めると式は簡単にはできない。術にもあるが，簡単な式にできたと思っても，それは特定の条件下のものに過ぎないと思われる。

大円と小円の直径が 10 寸，1 寸のとき，中円の直径は 5 寸である（r2 = 2.5）。

res[1](r1 =>10/2, r3 => 1/2).evalf() |> println
res[2](r1 =>10/2, r3 => 1/2).evalf() |> println
res[3](r1 =>10/2, r3 => 1/2).evalf() |> println

   2.50000000000000
   7.07106781186548
   2.82842712474619

術では 「+ 0.5」 とあるのを山村が 「- 0.1」 と訂正した。確かにそのようにすれば中円の直径として 5 が得られる。
しかし，これは「大円径 = 10，小円径 = 1」のときのみの付け焼き刃に過ぎず，一般解ではない。

大円径 = 10
小円径 = 1
A = (大円径 - 小円径)*小円径
B = √A/大円径 - 0.1  # 
中円径 = 小円径/B

   5.0

たとえば，大円径が 10，小円径が 0.8 のときの計算をすると破綻が露見する。
正しい答えは，「大円の直径が 10， 小円の直径が 0.8 のとき，中円の直径は 3.49792834161003」である。

その他のパラメータは以下の通りである。
r1 = 5;  r3 = 0.4;  r2 = 1.74896;  x1 = 5.91433;  x3 = 2.07766

2res[1](r1 =>10/2, r3 => 0.8/2).evalf() |> println

   3.49792834161003

山村の訂正では，中円の半径は 4.670354708490701 になる。どちらが正しいかは，図を見ればわかる。

大円径 = 10
小円径 = 0.8
A = (大円径 - 小円径)*小円径
B = √A/大円径 - 0.1  # 
中円径 = 小円径/B

   4.670354708490701

function draw(r1, r3, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(r1 - r3)
   r2 = 2*r1*(r1^3*r3 + 2*r1^2*r3^(3/2)*t - 4*r1^2*r3^2 - 8*r1*r3^(5/2)*t + 4*r1*r3^3 + 8*r3^(7/2)*t)/(r1 - 2*r3)^4
   x3 = 2*sqrt(2)*sqrt((r1^4*r3^2 + 4*r1^3*r3^(5/2)*t - 16*r1^2*r3^(7/2)*t - 16*r1^2*r3^4 + 16*r1*r3^(9/2)*t + 32*r1*r3^5 - 16*r3^6)/(r1^4 - 8*r1^3*r3 + 24*r1^2*r3^2 - 32*r1*r3^3 + 16*r3^4))
   x1 = r1*x3*(r1^4*r3^(5/2)*t - r1^4*r3^3 - 8*r1^3*r3^(7/2)*t + 8*r1^3*r3^4 + 24*r1^2*r3^(9/2)*t - 24*r1^2*r3^5 - 32*r1*r3^(11/2)*t + 32*r1*r3^6 + 16*r3^(13/2)*t - 16*r3^7)/(r3^3*(r1 - 2*r3)^5)
   y1 = 2r2 - r1
   @printf("大円の直径が %g， 小円の直径が %g のとき，中円の直径は %g である。\n", 2r1, 2r3, 2r2)
   @printf("r1 = %g;  r3 = %g;  r2 = %g;  x1 = %g;  x3 = %g\n", r1, r3, r2, x1, x3)
   plot()
   circle(x1, y1, r1)
   circle(0, r2, r2, :blue)
   circle(0, 2r2 - r3, r3, :green)
   circle(x3, 2r2 - 3r3, r3, :green)
   circle(x1, y1 + r1 - r3, r3, :green)
   hline!([2r2 - 2r3])
   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(x1, y1, "大円:r1,(x1,y1)", :red, :center, delta=-delta/2)
       point(0, r2, "中円:r2,(0,r2)", :blue, :center, delta=-delta/2)
       point(0, 2r2 - r3, "小円:r3,(0,2r2-r3)", :green, :left, :vcenter, delta=6delta)
       point(x3, 2r2 - 3r3, "小円:r3,(x3,2r2-3r3)", :green, :left, :vcenter, deltax=5delta)
       point(x1, y1 + r1 - r3, "小円:r3,(x1,y1+r1-r3)", :green, :left, :vcenter, delta=6delta)
   end
end;