算額（その233）

算額（その233）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（222)
長野県東筑摩郡筑北村青柳 碩水寺豊川稲荷社 慶応4年(1868)

菱形の中に，大円，中円，小円が入っている。中円，小円の径がそれぞれ 8寸，4寸のとき，大円の径を求めよ。

大円の半径 r1, (r3 + r1, 0),
中円の半径 r2, (0, r3 + r2),
小円の半径 r3, (0, 0),
だたし，r2 = 8/2, r3 = 4/2。
以下の方程式を解く。

r2, r3 が与えられているとき，条件式は 1 つである。
r2, r3 が任意の値を取るとしてこの方程式を解くと r1 は r2, r3 を含む式になる。

include("julia-source.txt")

using SymPy

@syms r1::positive, r2::positive, r3::positive;

(r2, r3) = (4, 2)
eq1 = (r3 + r1)^2 + (r3 + r2)^2 - (r1 + r2)^2
solve(eq1, r1)[1] |> println

   r3*(r2 + r3)/(r2 - r3)

すなわち，「中円と小円の半径の和を，中円と小円の半径の差で割り，小円の半径を掛ければ，大円の半径が得られる」
2 * (4 + 2) / (4 - 2) = 6 である。

大円の半径は6，すなわち元の単位での大円の直径は 12寸である。

外形の菱形を書くには，x軸上の頂点のx座標を x, y軸上の頂点のx座標を y として方程式を2本追加し，合わせて3本の連立方程式を解く。

ただし，r2, r3 が未知数のままだと solve() では解を求めることができないので，問にあるように r2 = 4, r3 = 2 の場合について解く。

using SymPy
@syms r1::positive, r2::positive, r3::positive, x::positive, y::positive;

(r2, r3) = (4, 2)
eq1 = (r3 + r1)^2 + (r3 + r2)^2 - (r1 + r2)^2
eq2 = distance(0, y, x, 0, r3 + r1, 0) - r1^2
eq3 = distance(0, y, x, 0, 0, r3 + r2) - r2^2;

res = solve([eq1, eq2, eq3], (r1, x, y))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (6, (13*sqrt(6) + 32)/(2*sqrt(7 - 2*sqrt(6))), 26/5 + 32*sqrt(6)/15)

x はもう少し簡単な式になる。

res[1][2] |> sympy.sqrtdenest |> simplify |> println

   11 + 9*sqrt(6)/2

y は通分ぐらいしか手がない。

res[1][3] |> factor |> println

   2*(39 + 16*sqrt(6))/15

   r1 = 6.0000000
   x = 22.0227038;  y = 10.4255781

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x, y) = (6, 11 + 9*sqrt(6)/2, (78 + 32*sqrt(6))/15)
   @printf("r1 = %.7f\n", r1)
   @printf("x = %.7f;  y = %.7f\n", x, y)
   plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], color=:black, lw=0.5)
   circle(r3 + r1, 0, r1, :red)
   circle(-r3 - r1, 0, r1, :red)
   circle(0, r3 + r2, r2, :blue)
   circle(0, -r3 - r2, r2, :blue)
   circle(0, 0, r3, :green)
   if more == true
       point(r3, 0, " r3", :green)
       point(0, 0, "小円", :green, :center, :top, mark=false)
       point(r3 + r1, 0, " r3+r1", :red)
       point(r3 + r1, r3, "大円", :red, mark=false)
       point(x, 0, " x", :black)
       point(0, r3 + r2, "r3+r2", :blue, :right, :bottom)
       point(0, r3 + r2, " 中円", :blue, :left, :bottom, mark=false)
       point(0, y, "y ", :black, :right, :bottom)
       point(0, r3, "r3", :blue, :right, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;