算額（その793）

算額（その793）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

直線上に互いに接する 3 個の円（甲円，丙円，丁円）が載っており，その 3 個の円のいずれとも接するように乙円が載っている。甲円，丙円，丁円の直径がわかっているときに，乙円の直径を求めよ。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (0, r3)
丁円の半径と中心座標を r4, (x4, r4)
とおき，以下の連立方程式を解く。
なお，SymPy の能力では一般解は求まらないので，与えられた条件における解を求める。なお，後半で r2 の一般解を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     x2::positive, y2::positive, r3::positive,
     r4::positive, x4::positive
(r1, r3, r4) = (100, 64, 48) .// 2
eq1 = (x1 - x2)^2 + (y2 - r1)^2 - (r1 + r2)^2 |> expand
eq2 = (x1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^2 |> expand
eq3 = x2^2 + (y2 - r3)^2 - (r2 + r3)^2 |> expand
eq4 = (x4 - x2)^2 + (y2 - r4)^2 - (r2 + r4)^2 |> expand
eq5 = x4^2 + (r3 - r4)^2 - (r3 + r4)^2 |> expand;
res = solve([eq1, eq2, eq3, eq4, eq5], (x1, r2, x2, y2, x4))
(x1, r2, x2, y2, x4) = res[1]
@printf("乙円の直径 = %g\n", 2res[1][2])
println("(x1, r2, x2, y2, x4) = ",res[1])
@printf("x1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  x4 = %g\n", x1, r2, x2, y2, x4)

   乙円の直径 = 72.9
   (x1, r2, x2, y2, x4) = (72*sqrt(3), 729/20, 26*sqrt(3), 1671/20, 32*sqrt(3))
   x1 = 124.708;  r2 = 36.45;  x2 = 45.0333;  y2 = 83.55;  x4 = 55.4256

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3, r4) = (100, 64, 48) .// 2
   (x1, r2, x2, y2, x4) = (72*sqrt(3), 729/20, 26*sqrt(3), 1671/20, 32*sqrt(3))
   plot()
   circle(x1, r1, r1)
   circle(x2, y2, r2, :blue)
   circle(0, r3, r3, :green)
   circle(x4, r4, r4, :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(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
       point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, delta=-delta)
       point(0, r3, "丙円:r3,(0,r3)", :green, :center, delta=-delta)
       point(x4, r4, "丁円:r4\n(x4,r4)", :magenta, :center, delta=-delta)
   end
end;

1.　算法助術の公式47

直線上に載って互いに接している円の半径をそれぞれ r1，r2 とする。その 2 つの円に載っている円の半径を r3 とする。
上に載っている円のてっぺんから直線までの距離 h を求める。

公式集は，直径についての式になっているので，これを半径についての式に書き直す。
それぞれの円の直径を d1 = 2r1, d2 = 2r2, d3 = 2r3 とする。

d1*h + d2*h - d1*d2 = 2sqrt(d1*d2*d3*h)
r1*h + r2*h - 2r1*r2 = 2sqrt(2r1*r2*r3*h)

using SymPy
@syms r1, r2, r3, r4, h, d1, d2, d3, d4
eq47d = d1*h + d2*h - d1*d2 - 2sqrt(d1*d2*d3*h)
eq47r = 2r1*h + 2r2*h - 2r1*2r2 - 2sqrt(2r1*2r2*2r3*h)
eq47r = r1*h + r2*h - 2r1*r2 - 2sqrt(2r1*r2*r3*h)

2.　応用

乙円のてっぺんから直線までの距離を h とする。
甲円，乙円，丙円，丁円の半径をそれぞれ r1, r2, r3, r4 とする。
まず，甲円，丁円，乙円に対して公式47を適用する。
ただし，そのまま使って連立方程式を解くには SymPy の能力的に無理があるので，辺々を 2 乗して使う。

using SymPy
@syms r1, r2, r3, r4, h

eq11 = (r1*h + r4*h - 2r1*r4)^2 - 4(2r1*r4*r2*h);

同様に，丙円，丁円，乙円に対して公式47を適用する。

eq12 = (r3*h + r4*h - 2r3*r4)^2 - 4(2r3*r4*r2*h);

筆算で解く場合にはまず h について解き，それをどちらかの式に代入して r2 を求めるのだろうが，SymPy で連立方程式を解くと，自動的にやってくれる。

res = solve([eq11, eq12], (h, r2));
res[2]

   (4*r1*r3*(r4^2*(r1 + r3 + 2*r4)/(2*(r1*r3 - r4^2)) + r4 + r4^2*sqrt(r1*r3)*(2*r1*r3 + r1*r4 + r3*r4)/(2*r1*r3*(r1*r3 - r4^2)))/(2*r1*r3 + r1*r4 + r3*r4), r4^2*(r1 + r3 + 2*r4)/(4*(r1*r3 - r4^2)) + r4^2*sqrt(r1*r3)*(2*r1*r3 + r1*r4 + r3*r4)/(4*r1*r3*(r1*r3 - r4^2)))

2 番目のものが適解である。

甲円の直径が 100 寸, 丙円 の直径が 64 寸, 丁円の直径が 48 寸のとき, h =120，r2 = 36.45 である。直径は 72.9 である。

res[2][1](r1 => 100/2, r3 => 64/2, r4 => 48/2) |> println
2res[2][2](r1 => 100/2, r3 => 64/2, r4 => 48/2) |> println

   120.000000000000
   72.9000000000000

r2 の式は，長いままで，SymPy では簡約化できない。

res[2][2] |> println

   r4^2*(r1 + r3 + 2*r4)/(4*(r1*r3 - r4^2)) + r4^2*sqrt(r1*r3)*(2*r1*r3 + r1*r4 + r3*r4)/(4*r1*r3*(r1*r3 - r4^2))

術文は，「天 = sqrt(甲*丙), 地 = 4(天 - 丁)として，乙 = ((甲 + 丙)/天 + 2)丁^2/地

@syms d1::positive, d3::positive, d4::positive
h = sqrt(d1*d3)d4/(sqrt(d1*d3) - d4);
d2 = ((d1 + d3)/sqrt(d1*d3) + 2)*d4^2/4(sqrt(d1*d3)- d4)
d2(d1 => 100, d3 => 64, d4 => 48).evalf() |> println

   72.9000000000000