算額（その1266）

算額（その1266）

百三十 現存するが奉納場所不明 明治14年(1881)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，菱形

菱形の中に大円 1 個，中円 2 個，小円 2 個を容れる。菱長（菱形の対角線の長い方）と大円の直径が与えられたとき，小円の直径を求めるすべを述べよ。

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

include("julia-source.txt");

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

   (a*r1*sqrt(1/(a - r1))/sqrt(a + r1), a*r1/(a + r1), a*r1*(a - r1)/(a^2 + 2*a*r1 + r1^2))

小円の半径 r3 はたいして簡約化できない。

@syms d
r3 = apart(res[3], d) |> factor
r3 |> println

   -a*r1*(-a + r1)/(a + r1)^2

小円の半径 r3 は，菱長 a，大円の半径 r1 の関数である。
これは，「術」の「大円径/(大円径 + 菱長)^2 * 菱長*(菱長 - 大円径)」に一致する。

たとえば，a = 5, r1 = 3 のとき，r3 = -a*r1*(-a + r1)/(a + r1)^2 = 15/32 = 0.46875 である。
菱長，大円の直径が 10, 6 のとき，小円の直径は 0.9375 である。

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

   a = 5;  r1 =　3;  b =　3.75;  r2 =　1.875;  r3 =　0.46875

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

draw(10/2, 6/2, true)