算額（その1323）

算額（その1323）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円3個，正方形，等脚台形
#Julia, #SymPy, #算額, #和算

正方形の中に等脚台形と甲円，乙円，丙円を容れる。甲円の直径が 3 寸，丙円の直径が 1 寸のとき，乙円の直径はいかほどか。

正方形の一辺の長さを a
台形の頂点座標を (0, 0), (0, b), (c, a), (a, d)
甲円の半径と中心座標を r1, (a - r1, a - r1)
乙円の半径と中心座標を r2, (a - r2, r2)
丙円の半径と中心座標を r3, (r3, a - r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     d::positive, r1::positive, r2::positive, r3::positive;
eq1 = (a - c)^2 + (a - d)^2 - b^2
eq2 = ((a - b) + c - 2r3)^2 - ((a - b)^2 + c^2)
eq3 = ((a - d) + (a - c) - 2r1)^2 - ((a - d)^2 + (a - c)^2)
eq4 = (a + d - 2r2)^2 - (a^2 + d^2)
eq5 = (a - b)/c - d/a;

まず，eq1, eq2, eq5 から a, b, d を求める。

res = solve([eq1, eq2, eq5], (a, b, d))[2]  # 2 of 2

   (c*(c^2 - 2*r3^2)/(c^2 - 4*c*r3 + 2*r3^2), (c^2 - 2*c*r3 + 2*r3^2)^2/((c - 2*r3)*(c^2 - 4*c*r3 + 2*r3^2)), 2*r3*(c^3 - c^2*r3 - 2*c*r3^2 + 2*r3^3)/(c^3 - 6*c^2*r3 + 10*c*r3^2 - 4*r3^3))

eq3 の a, d に代入して c を求める。

eq13 = eq3(a => res[1], d => res[3]) |> simplify |> numerator |> factor
ans_c = solve(eq13, c)[1]
ans_c |> println
ans_c(r1 => 3/2, r3 => 1/2) |> println

   2*r3*(r1 - r3)/(r1 - 2*r3)
   2.00000000000000

eq4 に，a, b, d, c の順に代入し，r2 を求める。

eq14 = eq4(a => res[1], b => res[2], d => res[3])(c => ans_c) |> simplify |> numerator |> factor
ans_r2 = solve(eq14, r2)[2]  # 2 of 2
ans_r2 |> println
ans_r2(r1 => 3/2, r3 => 1/2) |> println

   r3*(-r1^2 + 2*r3^2)/(r1^2 - 4*r1*r3 + 2*r3^2)
   3.50000000000000

r2 は r1, r3 の関数である。甲円の直径が 3 寸，丙円の直径が 1 寸のとき，乙円の直径は 7 寸である。

function draw(r1, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   c = 2*r3*(r1 - r3)/(r1 - 2*r3)
   r2 = r3*(-r1^2 + 2*r3^2)/(r1^2 - 4*r1*r3 + 2*r3^2)
   (a, b, d) = (c*(c^2 - 2*r3^2)/(c^2 - 4*c*r3 + 2*r3^2), (c^2 - 2*c*r3 + 2*r3^2)^2/((c - 2*r3)*(c^2 - 4*c*r3 + 2*r3^2)), 2*r3*(c^3 - c^2*r3 - 2*c*r3^2 + 2*r3^3)/(c^3 - 6*c^2*r3 + 10*c*r3^2 - 4*r3^3))
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
   plot!([0, a, c, 0, 0], [0, d, a, b, 0], color=:orange, lw=0.5)
   circle(a - r1, a - r1, r1)
   circle(a - r2, r2, r2, :blue)
   circle(r3, a - r3, 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, a, "(a,a) ", :green, :right, :bottom, delta=delta/2)
       point(a, d, "(a,d)  ", :green, :right, :vcenter)
       point(c, a, "(c,a)", :green, :center, :bottom, delta=delta/2)
       point(0, b, "b ", :green, :right, :vcenter)
       point(a - r1, a - r1, "甲円:r1\n(a-r1,a-r1)", :red, :center, :bottom, delta=delta/2)
       point(a - r2, r2, "乙円:r2,(a-r2,r2)", :blue, :center, :bottom, delta=delta/2)
       point(r3, a - r3, " 丙円:r3,(r3,a-r3)", :magenta, :left, delta=-3delta)
       xlims!(-5delta, a + 3delta)
   end  
end;

draw(3/2, 1/2, true)