算額（その1126）

算額（その1126）

四十一 岩手県一関市牧沢 牧沢八幡神社 明治8年(1875)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円4個，二等辺三角形

二等辺三角形の中に，甲円，乙円，丙円を容れる。甲円と丙円の直径がわかっているときに，乙円の直径を求めよ。

与えられるものと求めるものが違うが，図形としては 算額（その164）と同じである。

二等辺三角形の底辺の長さを 2x，高さを y
甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (0, 2r1 + r2)
丙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。

一度に解くと有限の時間内に解が求まらないので，逐次解いていく。

include("julia-source.txt")

using SymPy

@syms r1::positive, r2::positive, r3::positive, x3::positive, x::positive, y::positive;

eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2 |> simplify  # 黃円と赤円が外接する
eq2 = dist2(0, y, x, 0, 0, 2r1 + r2, r2)  # 円の中心から斜辺までの距離
eq3 = dist2(0, y, x, 0, 0, r1, r1)  # 円の中心から
eq4 = r3*x - r1*(x - x3);  # 三角形の相似

まず，eq1 を解いて x3 を求める。

ans_x3 = solve(eq1, x3)[1]
ans_x3 |> println

   2*sqrt(r1)*sqrt(r3)

eq4 の x3 に ans_x3 を代入して，x を求める。

eq4 = eq4(x3 => ans_x3)
ans_x = solve(eq4, x)[1]
ans_x |> println

   2*r1^(3/2)*sqrt(r3)/(r1 - r3)

同様にして，eq3 をとき y を求める。

eq3 = eq3(x3 => ans_x3, x=> ans_x);
ans_y = solve(eq3, y)[1]
ans_y |> println

   8*r1^2*r3/(4*r1*r3 - (r1 - r3)^2)

x3, x, y を eq2 に代入して，r1 を求める。

eq2 = eq2(x3 => ans_x3, x=> ans_x, y => ans_y) |> simplify |> numerator |> (x -> x/(16r1^4*r3))
ans_r2 = solve(eq2, r2)[1] |> factor
ans_r2 |> println

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

「甲円と丙円の半径の差を二乗し，丙円の半径の4倍で割る」という，「術」と同じ式が得られる。

以上をまとめると，次のようになる。
r2 = r1^2/(4*r3) - r1/2 + r3/4
y  = 8r1^2*r3/(4r1*r3 - (r1 - r3)^2)
x  = 2r1^(3/2)*sqrt(r3)/(r1 - r3)
x3 = 2*sqrt(r1*r3)

冒頭の図は，甲円と丙円の半径がそれぞれ 2，1/2 の場合のものである。

甲円の直径が 4，丙円の直径が 1 のとき，乙円の直径は 2.25 である。
r1 = 2;  r2 = 1.125;  r3 = 0.5;  x3 = 2;  x = 2.6666666666666674;  y = 9.142857142857142

function draw(r1, r3, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = (r1 - r3)^2/4r3
   y  = 8r1^2*r3/(4r1*r3 - (r1 - r3)^2)
   x  = 2r1^(3/2)*sqrt(r3)/(r1 - r3)
   x3 = 2sqrt(r1*r3)
   @printf("甲円の直径が %g，丙円の直径が %g のとき，乙円の直径は %g である。\n", 2r1, 2r3, 2r2)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  x = %g;  y = %g\n", r1, r2, r3, x3, x, y)
   plot([x, 0, -x, x], [0, y, 0, 0], color=1, lw=0.5)
   circle(0, r1, r1, :orange)
   circle(0, 2r1 + r2, r2, :black)
   circle2(x3, r3, r3)
   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, r1, " r1", :orange, :left, :vcenter)
       point(0, 2r1 + r2, "2r1+r2", :black, :center, delta=-delta/2)
       point(x3, r3, "(x3,r3)", :red, :center, delta=-delta/2)
       point(x, 0, " x", :blue, :left, :bottom, delta=delta/2)
       point(0, y, " y", :blue, :left, :bottom, delta=delta/2)
   end
end;