算額（その1293）

算額（その1293）

十七 群馬県高崎市八幡町 八幡宮 文化7年(1810)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，菱形

菱形の中に，大円 2 個，中円 2 個，小円 1 個を容れる。菱形の対角線の長い方が 12 寸，短い方が 5 寸のとき，小円の直径はいかほどか。

菱形の対角線を 2a, 2b; a > b
大円の半径と中心座標を r1, (r3 + r1, 0)
中円の半径と中心座標を r2, (0, r3 + r2)
小円の半径と中心座標を r3, (0, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive, r3::positive;
sinθ = b/sqrt(a^2 + b^2)
cosθ = a/sqrt(a^2 + b^2)
eq1 = (r3 + r1)^2 + (r3 + r2)^2 - (r1 + r2)^2 |> expand
eq2 = (a - r3 - r1)*sinθ - r1
eq3 = (b - r3 - r2)*cosθ - r2;

SymPy では能力的に一度に解けないので，逐次解いていく。
まず，eq1 を解いて r1 を求める。

ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println

   r3*(r2 + r3)/(r2 - r3)

eq2 の r1 に，上で求めた r1 を代入し，r2 を求める。

eq12 = eq2(r1 => ans_r1)
ans_r2 = solve(eq12, r2)[1]
ans_r2 |> println

   r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))

eq3 に，上で求めた r1, r2 を代入し，r3 を求める。

eq13 = eq3(r1 => ans_r1, r2 => ans_r2) |> simplify |> numerator
ans_r3 = solve(eq13, r3)[2] |> factor
ans_r3 |> println

   -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2

解として得られた r3 を 代入して r2，r1 を確定する。

a = 12/2
b = 5/2
r3 = -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
r2 = r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
r1 = r3*(r2 + r3)/(r2 - r3)
(r1, r2, r3) |> println

   (1.5296184351192643, 0.9631806558860881, 0.49337363357064845)

菱形の対角線の長い方が 12 寸，短い方が 5 寸のとき，小円の直径は 2*0.49337363357064845 = 0.9867472671412969 寸である。「答」では「小円の径一寸三分二厘半微弱」としているが，正しいとは思えない。

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

   a = 6;  b = 2.5;  r1 = 1.52962;  r2 = 0.963181;  r3 = 0.493374

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
   r2 = r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
   r1 = r3*(r2 + r3)/(r2 - r3)
   @printf("菱形の対角線の長い方が %g，短い方が %g のとき，小円の直径は %g である。\n", 2a, 2b, 2r3)
   @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g;  r3 = %g\n", a, b, r1, r2, r3)
   plot([a , 0, -a, 0, a], [0, b, 0, -b, 0], color=:green, lw=0.5)
   circle2(r3 + r1, 0, r1)
   circle22(0, r3 + r2, r2, :blue)
   circle(0, 0, r3, :magenta)
   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(a, 0, " a", :green, :left, :vcenter)
       point(0, b, "b", :green, :center, :bottom, delta=delta)
       point(r3 + r1, 0, "大円:r1\n(r3+r1,0)", :red, :center, delta=-delta)
       point(0, r3 + r2, "中円:r2\n(0,r3+r2)", :blue, :center, :bottom, delta=delta)
       point(0, 0, "小円:r3,(0,0)", :magenta, :right, :vcenter, deltax=-7delta)
   end  
end;

draw(12/2, 5/2, true)