算額（その794）

算額（その794）

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

福島県郡山市日和田町宮下 天満宮 明治14年(1881)
http://www.wasan.jp/fukusima/miyasitatenma.html

後者の写真が不鮮明で長い間お蔵入りになっていたが，「精要算法」にあった。問題文と答えも同じである。ただし，後者の図では乙円と戊円は接していないが，精要算法にあるように，乙円と戊円も接している。

外円内に 7 個の円が入っている。外円，甲円の直径がそれぞれ 890 寸，267 寸のとき，乙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (0, r4 - R)
戊円の半径と中心座標を r5, (0, r5 + 2r4 - R)
とおき，以下の連立方程式を解く。

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

using SymPy

@syms R::positive, x2::positive, y2::positive, x3::positive, y3::negative, r1::positive, r2::positive, r3::positive, r4::positive, r5::positive;
# (R, r1) = (890,267) .// 2
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq3 = x3^2 + (r5 + 2r4 - R - y3)^2 - (r3 + r5)^2
eq4 = x3^2 + (y3 - r4 + R)^2 - (r3 + r4)^2
eq5 = x2^2 + y2^2 - (R - r2)^2
eq6 = x3^2 + y3^2 - (R - r3)^2
eq7 = x2^2 + (r5 + 2r4 - R - y2)^2 - (r2 + r5)^2
eq8 = r1 + r4 + r5 - R;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8],
    (x2, y2, x3, y3, r2, r3, r4, r5));
res[2] |> println

   (4*sqrt(3)*R*r1*(R - r1)/(3*R^2 - 2*R*r1 + 3*r1^2), -R*(-3*R^2 + 6*R*r1 + r1^2)/(3*R^2 - 2*R*r1 + 3*r1^2), -4*sqrt(3)*R*r1*(-R + r1)/(R^2 + 2*R*r1 + 9*r1^2), -R*(-R^2 + 2*R*r1 + 11*r1^2)/(R^2 + 2*R*r1 + 9*r1^2), -4*R*r1*(-R + r1)/(3*R^2 - 2*R*r1 + 3*r1^2), -4*R*r1*(-R + r1)/(R^2 + 2*R*r1 + 9*r1^2), -R*(-R + r1)/(R + 3*r1), (3*R*r1 - 3*r1^2)/(R + 3*r1))

2 組の解が得られるが，2 番目のものが適解である。

外円，甲円の直径がそれぞれ 890 寸，267 寸のとき，乙円の直径は 8*R*r1*(R - r1)/(3*R^2 - 2*R*r1 + 3*r1^2) = 280 寸である。

2res[2][5] |> println
2res[2][5](R => 890//2, r1 => 267//2) |> println

   -8*R*r1*(-R + r1)/(3*R^2 - 2*R*r1 + 3*r1^2)
   280

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1) = (890,267) ./ 2
   (x2, y2, x3, y3, r2, r3, r4, r5) = (4*sqrt(3)*R*r1*(R - r1)/(3*R^2 - 2*R*r1 + 3*r1^2), -R*(-3*R^2 + 6*R*r1 + r1^2)/(3*R^2 - 2*R*r1 + 3*r1^2), -4*sqrt(3)*R*r1*(-R + r1)/(R^2 + 2*R*r1 + 9*r1^2), -R*(-R^2 + 2*R*r1 + 11*r1^2)/(R^2 + 2*R*r1 + 9*r1^2), -4*R*r1*(-R + r1)/(3*R^2 - 2*R*r1 + 3*r1^2), -4*R*r1*(-R + r1)/(R^2 + 2*R*r1 + 9*r1^2), -R*(-R + r1)/(R + 3*r1), (3*R*r1 - 3*r1^2)/(R + 3*r1))
   plot()
   circle(0, 0, R, :black)
   circle(0, R - r1, r1)
   circle(x2, y2, r2, :green)
   circle(-x2, y2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(-x3, y3, r3, :magenta)
   circle(0, r4 - R, r4, :brown)
   circle(0, r5 + 2r4 - R, r5, :blue)
   if more == true
       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, R, "外円:R,(0,0)", :red, :center, :bottom, delta=delta/2)
       point(0, R-r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta)
       point(x2, y2, "乙円:r2,(x2,y2)", :green, :center, delta=-delta)
       point(x3, y3, "丙円:r3,(x3,y3)", :magenta, :center, delta=-delta)
       point(0, r4 - R, "丁円:r4,(0,r4-R)", :brown, :center, delta=-delta)
       point(0, r5 + 2r4 - R, "戊円:r5,(0,r5+2r4-R)", :blue, :center, delta=-delta)
   end
end;

 

算額（その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