算額（その270）

算額（その270）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額（260）
長野県下高井郡木島平村 天満宮 推定明治21年(1888)

直線の上に丁円，乙円，甲円，乙円の順に 4 円が外接している。この 4 円に外接する丙円がある。甲円，乙円の径がそれぞれ 2 寸，3 寸のとき，丁円の径を求めよ。
丁円の半径，中心座標を r0, (0, r0)
甲円の半径，中心座標を r1, (x1, r1)
乙円の半径，中心座標を r2, (x2, r2), (x3, r2), x3 = x1 + (x1 - x3) = 2x1 - x2
丙円の半径，中心座標を r4, (x1, 2r1 + r4)

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive, r4::positive, x1::positive, x2::positive, x3::positive, y4::positive;

x2 = 2x1 - x3
x4 = x1
y4 = 2r1 + r4
eq1 = x3^2 + (r0 - r2)^2 - (r0 + r2)^2
eq2 = x1^2 + (r0 - y4)^2 - (r0 + r4)^2
eq3 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = (x2 - x4)^2 + (r2 - y4)^2 - (r2 + r4)^2;

r1, r2 を変数として，r0, r4, x1, x3 を求める。

res = solve([eq1, eq2, eq3, eq4], (r0, r4, x1, x3))

   2-element Vector{NTuple{4, Sym}}:
    (r1*r2^2/(2*r1 - r2)^2, -r1^2/(r1 - r2), -4*r1^(3/2)*sqrt(r2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)), -2*sqrt(r1)*r2^(3/2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)))
    (r1*r2^2/(2*r1 - r2)^2, -r1^2/(r1 - r2), 4*r1^(3/2)*sqrt(r2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)), 2*sqrt(r1)*r2^(3/2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)))

二組の解が求まるが，実際の値を与えて検証すると，二番目の解が適切である。

   丁円: r0 = 18.000000;  丙円: r4 = 4.000000;  x1 = 19.595918;  x3 = 14.696938

(r1, r2) = (2, 3)
(r1*r2^2/(2*r1 - r2)^2, -r1^2/(r1 - r2), -4*r1^(3/2)*sqrt(r2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)), -2*sqrt(r1)*r2^(3/2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)))

   (18.0, 4.0, -19.595917942265427, -14.69693845669907)

(r1, r2) = (2, 3)
(r1*r2^2/(2*r1 - r2)^2, -r1^2/(r1 - r2), 4*r1^(3/2)*sqrt(r2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)), 2*sqrt(r1)*r2^(3/2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)))

   (18.0, 4.0, 19.595917942265427, 14.69693845669907)

二番目の解の数式表現を求める。

甲円，乙円の半径を引数として，丁円の半径を求めるプログラムを書く。

直径を与えれば直径を求める。

draw(2, 3, more=false

   丁円: r0 = 18.000000;  丙円: r4 = 4.000000;  x1 = 19.595918;  x3 = 14.696938

using Plots

function draw(r1, r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r4, x1, x3) = (
       r1*r2^2/(2*r1 - r2)^2,
       -r1^2/(r1 - r2),
       4*r1^(3/2)*sqrt(r2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)),
       2*sqrt(r1)*r2^(3/2)*sqrt(1/(4*r1^2 - 4*r1*r2 + r2^2)))
   x2 = 2x1 - x3
   x4 = x1
   y4 = 2r1 + r4
   @printf("丁円: r0 = %.6f;  丙円: r4 = %.6f;  x1 = %.6f;  x3 = %.6f\n", r0, r4, x1, x3)
   plot()
   circle(0, r0, r0)
   circle(x1, r1, r1, :blue)
   circle(x3, r2, r2, :green)
   circle(x2, r2, r2, :green)
   circle(x4, y4, r4, :magenta)
   if more
       point(x1, r1, " 甲円:r1(x1,r1)", :blue)
       point(x2, r2, " 乙円:r2(x2,r2)", :green, :left, :bottom)
       point(x3, r2, "乙円:r2(x3,r2) ", :green, :right, :bottom)
       point(x4, y4, " 丙円:r4(x4,y4)", :magenta, :left, :bottom)
       point(0, r0, " 丁円:r0", :red)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5, xlims=(-r0, x2 + 4r2))
   else
       plot!(showaxis=false)
   end
end;

算額（その269）

算額（その269）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額（260）
長野県下高井郡木島平村 天満宮 推定明治21年(1888)

追記：この算額の図は多分誤っている。算額（その798）の図が正しい。https://blog.goo.ne.jp/r-de-r/e/1fd7e974eb24c0c5118e308a9e790003
すなわち，「小円は 3 個の中円に外接している」のが正しい。そのように仮定すれば，答えも術も矛盾しない。

大円 2 個，中円 5 個，小円 4 個が配置されている。小円は大円と中円に外接している。
大円の直径が 1 尺（10 寸）のとき，小円の直径を求めよ。

術は「近似値」を求めているが，右上と左下の小円の中心間の距離を「大円の半径+小円の直径」とする無茶なものである。

大円，中円の半径を r1, r2 とおく（r2 = r1/2）。小円の直径を r3, 中心座標を (x3, x3) とおき，以下の方程式を解く。条件式が 2 本で，r3, x3 の二変数について解を求める。その解の式には r1 が変数として含まれる。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, x3::positive;

r2 = r1//2
eq1 = x3^2 + (x3 + r2)^2 - (r1 + r3)^2
eq2 = (2r2 - x3)^2 + x3^2 - (r2 + r3)^2;

solve([eq1, eq2], (r3, x3))

   1-element Vector{Tuple{Sym, Sym}}:
    (3*r1/14, 4*r1/7)

直径が 10 寸ならば，小円の直径は 30/14 = 2.142857142857143 である。
術では (sqrt(1/2) - 1/2)*10 = 2.0710678118654755 としているので，そうとう違う。

追記：この術は算額（その498）の図についての正しい答えである。

(sqrt(1/2) - 1/2)*10

   2.0710678118654755

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 10//2
   r2 = r1//2
   (r3, x3) = (3r1/14, 4r1/7)
   @printf("r3 = %.6f;  x3 = %.6f\n", r3, x3)
   @printf("小円の直径 = 2r3 = %.6f\n", 2r3)
   plot()
   circle(0, r2, 2r2)
   circle(0, -r2, 2r2)
   circle(0, 0, r2, :blue)
   circle42(0, 2r2, r2, :blue)
   circle4(x3, x3, r3, :green)
   if more
       point(r2, 0, " r2")
       point(0, r2, " r2")
       point(2r2, 0, " 2r2")
       point(0, 2r2, " 中円:r2\n 2r2", :blue)
       point(3r2, 0, "3r2 ", :green, :right)
       point(x3, x3, " 小円:r3\n (x3,x3)", :green, :center, :bottom)
       point(1.2r2, 2r2, "大円", :red, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;