算額（その138）

算額（その138）

兵庫県西宮市社家町 西宮えびす神社 天保13年(1842)　
http://www.wasan.jp/hyogo/nisinomiya.html

径が 10 寸の外円の中に，大円 2 個と，小円 3 円が入っている。小円の径を求めよ。

図のように，外円の半径を 1，大円と小円の半径を r1, r2 とし，方程式を解く。

using SymPy

@syms r1::positive, r2::positive;
eq1 = 2r1 - r2 - 1;  # 外円，大円，小円の半径の関係
eq2 = (1 - r2)^2 + (1 - r1)^2 - (r1 + r2)^2;  # 上の大円と右の小円が外接している

res = solve([eq1, eq2], (r1, r2))

   1-element Vector{Tuple{Sym, Sym}}:
    (-1/2 + sqrt(5)/2, -2 + sqrt(5))

(res[1][1].evalf(), res[1][2].evalf())

   (0.618033988749895, 0.236067977499790)

元の問題では外円の直径が 10 寸である。小円は外円の 0.236067977499790 倍。すなわち，小円の径は約 2寸3分6厘 である。

図に描いたとき，真ん中の小円が両脇の小円に比べてやや小さく見えるのは，錯視の効果である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, lw=0.5, linestyle=:solid)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=lw, linestyle=linestyle)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = res[1]
   println("大円の半径 = $r1 = $(r1.evalf())")
   println("小円の半径 = $r2 = $(r2.evalf())")
   plot()
   circle(0, 0, 1)
   circle(0, 1 - r1, r1, :green)
   circle(0, r1 - 1, r1, :green)
   circle(0, 0, r2, :blue)
   circle(1 - r2, 0, r2, :blue)
   circle(r2 - 1, 0, r2, :blue)
   if more == true
       point(0, 1 - r1, "1-r1 ", :green)
       point(1 - r2, 0, "1-r2 ", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

draw(false)

   大円の半径 = -1/2 + sqrt(5)/2 = 0.618033988749895
   小円の半径 = -2 + sqrt(5) = 0.236067977499790

算額（その137）

算額（その137）

千葉県千葉市中央区 延命寺 嘉永5(1852)年　
http://www.wasan.jp/chiba/enmeiji.html

外円の中に，上円，中円，下円が入っている。上円と下円の径が与えられたとき，外円の径はいかほどか。

外円，上円，中円，下円の半径を r0, r1, r2, r3，中円の中心の座標を (x1, y2) として，図のように記号を定め方程式を解く。

using SymPy

@syms r0::positive, r1::positive, r2::positive, y2::negative, r3::positive;
eq1 = r2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2;  # 上円と中円が外接する
eq2 = r2^2 + (y2 - r3 + r0)^2 - (r2 + r3)^2;  # 中円と下円が外接する
eq3 = r2^2 + y2^2 - (r0 - r2)^2;              # 外円に中円が内接する

res = solve([eq2, eq3, eq1], (r0, r2, y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    ((-r1^3 - 3*r1^2*r3 + r1^2*sqrt(r1^2 + 14*r1*r3 + r3^2) - 3*r1*r3^2 - 4*r1*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) - r3^3 + r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2))/(-r1^2 + 10*r1*r3 + r1*sqrt(r1^2 + 14*r1*r3 + r3^2) - r3^2 + r3*sqrt(r1^2 + 14*r1*r3 + r3^2)), 4*r1*r3*(r1^3 - 33*r1^2*r3 - r1^2*sqrt(r1^2 + 14*r1*r3 + r3^2) - 33*r1*r3^2 - 14*r1*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^3 - r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2))/(r1^4 - 12*r1^3*r3 - r1^3*sqrt(r1^2 + 14*r1*r3 + r3^2) + 22*r1^2*r3^2 + r1^2*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) - 12*r1*r3^3 + r1*r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^4 - r3^3*sqrt(r1^2 + 14*r1*r3 + r3^2)), (-r1^2 + 10*r1*r3 + r1*sqrt(r1^2 + 14*r1*r3 + r3^2) - r3^2 + r3*sqrt(r1^2 + 14*r1*r3 + r3^2))/(2*r1 - 2*r3))
    ((r1^3 + 3*r1^2*r3 + r1^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + 3*r1*r3^2 - 4*r1*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^3 + r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2))/(r1^2 - 10*r1*r3 + r1*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^2 + r3*sqrt(r1^2 + 14*r1*r3 + r3^2)), 4*r1*r3*(r1^3 - 33*r1^2*r3 + r1^2*sqrt(r1^2 + 14*r1*r3 + r3^2) - 33*r1*r3^2 + 14*r1*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^3 + r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2))/(r1^4 - 12*r1^3*r3 + r1^3*sqrt(r1^2 + 14*r1*r3 + r3^2) + 22*r1^2*r3^2 - r1^2*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) - 12*r1*r3^3 - r1*r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^4 + r3^3*sqrt(r1^2 + 14*r1*r3 + r3^2)), (-r1^2 + 10*r1*r3 - r1*sqrt(r1^2 + 14*r1*r3 + r3^2) - r3^2 - r3*sqrt(r1^2 + 14*r1*r3 + r3^2))/(2*r1 - 2*r3))

最初の解は不適切解である。

2番めの解の r0 もかなり複雑な式になり，simplify() でも簡単にできない。簡約化は，このページの最後に示す。
r0 = (r1^3 + 3*r1^2*r3 + r1^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + 3*r1*r3^2 - 4*r1*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^3 + r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2))/(r1^2 - 10*r1*r3 + r1*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^2 + r3*sqrt(r1^2 + 14*r1*r3 + r3^2))

sqrt(r1^2 + 14*r1*r3 + r3^2) が何回も出てくるので，term = sqrt(r1^2 + 14*r1*r3 + r3^2) とすると，若干短くはなる。
r0 = (r1^3 + 3*r1^2*r3 + r1^2*term + 3*r1*r3^2 - 4*r1*r3*term + r3^3 + r3^2*term)/(r1^2 - 10*r1*r3 + r1*term + r3^2 + r3*term)

「算額の意義」には，計算過程抜きで，上円と下円の直径を与えて外円の直径を求める式が書かれている。

記号を会わせると以下のようになる。
A = r3 / r1, B = (A + 1) / 2 とすると
r0 = (√(B^2 + 3A) + B) × r1

この r0 を求める式を A, B を用いずに書き下すと
r1/2 + r3/2 + sqrt(r1^2 + 14*r1*r3 + r3^2)/2
となり，sqrt(r1^2 + 14*r1*r3 + r3^2) が現れるので，SymPy による長い式も簡単になることと思われる。このページの最後に示す。

SymPy による r0 を求める式による解と両方を求める関数を書いて，正しいことを確認すると，両者は同じ答えを返すことが分かった。

いずれにしろ，外円の直径を求める場合に，中円の直径は不要であるということは特筆しておくべきことであろう。

function getr0(r1, r3)
   A = r3 / r1
   B = (A + 1) / 2
   result1 = (sqrt(B^2 + 3A) + B) * r1
   result2 = (r1^3 + 3*r1^2*r3 + r1^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + 3*r1*r3^2 - 4*r1*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^3 + r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2))/(r1^2 - 10*r1*r3 + r1*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^2 + r3*sqrt(r1^2 + 14*r1*r3 + r3^2))
   term = sqrt(r1^2 + 14*r1*r3 + r3^2)
   result3 = (r1^3 + 3*r1^2*r3 + r1^2*term + 3*r1*r3^2 - 4*r1*r3*term + r3^3 + r3^2*term)/(r1^2 - 10*r1*r3 + r1*term + r3^2 + r3*term)
   return (result1, result2, result3)
end;
getr0(10, 5) |> println
getr0(40, 23) |> println
getr0(123.456, 78.9) |> println

   (21.86140661634507, 21.861406616345057, 21.861406616345057)
   (92.75561198780075, 92.7556119878007, 92.7556119878007)
   (299.8209532704345, 299.82095327043413, 299.82095327043413)

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, lw=0.5, linestyle=:solid)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=lw, linestyle=linestyle)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(r1, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   term = sqrt(r1^2 + 14*r1*r3 + r3^2)
   r0 = (r1^3 + 3*r1^2*r3 + r1^2*term + 3*r1*r3^2 - 4*r1*r3*term + r3^3 + r3^2*term)/(r1^2 - 10*r1*r3 + r1*term + r3^2 + r3*term)
   r2 = 4*r1*r3*(r1^3 - 33*r1^2*r3 + r1^2*term - 33*r1*r3^2 + 14*r1*r3*term + r3^3 + r3^2*term)/(r1^4 - 12*r1^3*r3 + r1^3*term + 22*r1^2*r3^2 - r1^2*r3*term - 12*r1*r3^3 - r1*r3^2*term + r3^4 + r3^3*term)
   y2 =(-r1^2 + 10*r1*r3 - r1*term - r3^2 - r3*term)/(2*r1 - 2*r3)
   # 簡約化
   # r0 = r1/2 + r3/2 + sqrt(r1^2 + 14*r1*r3 + r3^2)/2
   # r2 = 4*r1*r3*(2*r1 + 2*r3 - sqrt(r1^2 + 14*r1*r3 + r3^2))/(r1 - r3)^2
   println("外円の半径 = $r0")
   println("中円の半径 = $r2")
   println("中円の中心座標 = ($r2, $y2)")
   plot()
   circle(0, 0, r0)
   circle(r2, y2, r2, :orange)
   circle(-r2, y2, r2, :orange)
   circle(0, r0 - r1, r1)
   circle(0, r3 - r0, r3, :blue)
   if more == true
       point(0, r0 - r1, "r0-r1 ", :red, :right)
       point(r2, y2, "(r2,y2;r2)", :orange, :center, :top)
       point(0, r3 - r0, "r3-r0", :blue, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

draw(60, 20, false)

   外円の半径 = 112.1110255092798
   中円の半径 = 47.33384694432096
   中円の中心座標 = (47.33384694432096, -44.222051018559554)

上円，下円の半径を変えてみる。

draw(23.5, 30.9, false)

   外円の半径 = 81.22119954240218
   中円の半径 = 40.184945549364414
   中円の中心座標 = (40.184945549364414, 8.315304744143447)

-----

SymPy により得られた r0 を簡約化する。つまり，分母を有理化する。有理化は自動的には行われず，途中でも単に simplify では式が簡約化されず expand, factor, simplify を適用順も考慮して組み合わせて式変形を行う。

r0 = (r1^3 + 3*r1^2*r3 + r1^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + 3*r1*r3^2 - 4*r1*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^3 + r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2))/(r1^2 - 10*r1*r3 + r1*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^2 + r3*sqrt(r1^2 + 14*r1*r3 + r3^2));

r0 を簡約化する。

まず，分母について整理する。

d = denominator(r0);
d |> println

   r1^2 - 10*r1*r3 + r1*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^2 + r3*sqrt(r1^2 + 14*r1*r3 + r3^2)

sqrt(r1^2 + 14*r1*r3 + r3^2) の項の符号を反転させたものを掛ける。

d1 = (r1^2 - 10*r1*r3 + r1*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^2 + r3*sqrt(r1^2 + 14*r1*r3 + r3^2)) *
    (r1^2 - 10*r1*r3 - r1*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^2 - r3*sqrt(r1^2 + 14*r1*r3 + r3^2));

簡約化と因子分解を行う。

d2 = d1 |> simplify |> factor;
d2 |> println

   -36*r1*r3*(r1 - r3)^2

次に分母にかけたものと同じものを分子にも掛ける。

n1 = numerator(r0) * ((r1^2 - 10*r1*r3 + r3^2) - (r1 + r3) * sqrt(r1^2 + 14*r1*r3 + r3^2));

一度展開してから因子分解する。

n2 = n1 |> expand |> factor;
n2 |> println

   -18*r1*r3*(r1 - r3)^2*(r1 + r3 + sqrt(r1^2 + 14*r1*r3 + r3^2))

r0 は n2/d2 である。

res = n2/d2;
println(res)

   r1/2 + r3/2 + sqrt(r1^2 + 14*r1*r3 + r3^2)/2

これが簡約化された式で，http://www.wasan.jp/chiba/enmeiji.html に示されたものと同じである。

式変形に間違いがないか確認する。

res(r1 => 30, r3 => 10).evalf()

   56.0555127546399

r0(r1 => 30, r3 => 10).evalf()

   56.0555127546399

同様に，r2 も簡約化する。

r2 = 4*r1*r3*(r1^3 - 33*r1^2*r3 + r1^2*sqrt(r1^2 + 14*r1*r3 + r3^2) - 33*r1*r3^2 + 14*r1*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^3 + r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2))/(r1^4 - 12*r1^3*r3 + r1^3*sqrt(r1^2 + 14*r1*r3 + r3^2) + 22*r1^2*r3^2 - r1^2*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) - 12*r1*r3^3 - r1*r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^4 + r3^3*sqrt(r1^2 + 14*r1*r3 + r3^2));

d1 = denominator(r2);
d1 |> println

   r1^4 - 12*r1^3*r3 + r1^3*sqrt(r1^2 + 14*r1*r3 + r3^2) + 22*r1^2*r3^2 - r1^2*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) - 12*r1*r3^3 - r1*r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^4 + r3^3*sqrt(r1^2 + 14*r1*r3 + r3^2)

分母の有理化（sqrt(r1^2 + 14*r1*r3 + r3^2) の項の符号を反転させたものを掛ける）

d2 = (r1^4 - 12*r1^3*r3 + r1^3*sqrt(r1^2 + 14*r1*r3 + r3^2) + 22*r1^2*r3^2 - r1^2*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) - 12*r1*r3^3 - r1*r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^4 + r3^3*sqrt(r1^2 + 14*r1*r3 + r3^2)) *
(r1^4 - 12*r1^3*r3 - r1^3*sqrt(r1^2 + 14*r1*r3 + r3^2) + 22*r1^2*r3^2 + r1^2*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) - 12*r1*r3^3 + r1*r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^4 - r3^3*sqrt(r1^2 + 14*r1*r3 + r3^2));

d3 = d2 |>simplify |> factor;
d3 |> println

   -36*r1*r3*(r1 - r3)^6

分子にも sqrt(r1^2 + 14*r1*r3 + r3^2 の項の符号を反転させたものを掛ける

n1 = numerator(r2) * (r1^4 - 12*r1^3*r3 - r1^3*sqrt(r1^2 + 14*r1*r3 + r3^2) + 22*r1^2*r3^2 + r1^2*r3*sqrt(r1^2 + 14*r1*r3 + r3^2) - 12*r1*r3^3 + r1*r3^2*sqrt(r1^2 + 14*r1*r3 + r3^2) + r3^4 - r3^3*sqrt(r1^2 + 14*r1*r3 + r3^2));

n2 = n1 |> expand |> factor
n2 |>  println

   -144*r1^2*r3^2*(r1 - r3)^4*(2*r1 + 2*r3 - sqrt(r1^2 + 14*r1*r3 + r3^2))

res = n2 / d3 |> factor;
res |> println

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

式変換が正しく行われたことを確認する

res(r1 => 30, r3 => 10).evalf()

   23.6669234721606

r2(r1 => 30, r3 => 10).evalf()

   23.6669234721606