算額（その1046）

算額（その1046）

八十二 藤沢町保呂羽保呂羽山 保呂羽神社 明治7年(1874)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

正三角形の中に全円，大円，小円と斜線を容れる。小円の直径が 1 寸のとき，大円の直径はいかほどか。

正三角形の一辺の長さを 2a
全円の半径と中心座標を r0, (0, r0); r0 = √3a/3
大円の半径と中心座標を r1, (x1, r1); x1 < 0
小円の半径と中心座標を r2, (x2, r2)
斜線と斜辺の交点座標を (b, √3(b + a)); b < 0
とおき，以下の連立方程式を解く。

一度に解けないので，逐次解いてゆく。

include("julia-source.txt");

using SymPy
@syms a::positive, b::negative, r0::positive,
     r1::positive, x1::negative, r2::positive, x2::positive,
     y2::positive
r0 = √Sym(3)a/3
eq1 = dist2(a, 0, 0, √Sym(3)a, x2, y2, r2)
eq2 = dist2(-a, 0, 0, √Sym(3)a, x1, r1, r1)
eq3 = dist2(a, 0, b, √Sym(3)*(b + a), x2, y2, r2)
eq4 = dist2(a, 0, b, √Sym(3)*(b + a), x1, r1, r1)
eq5 = x2^2 + (r0 - y2)^2 - (r0 + r2)^2
eq6 = x1^2 + (r0 - r1)^2 - (r0 + r1)^2;

eq2, eq6 から a, x1 を求める。

res = solve([eq2, eq6], (a, x1))[2]

   (3*sqrt(3)*r1, -2*sqrt(3)*r1)

eq1, eq3, eq5, eq6 に a, x1 を代入する。

eq1 = 4eq1(a => res[1]);
eq3 = numerator(eq3(a => res[1]));
eq4 = numerator(eq4(a => res[1], x1 => res[2]));
eq5 = eq5(a => res[1]);

eq1 |> println
eq3 |> println
eq4 |> println
eq5 |> println

   81*r1^2 - 18*sqrt(3)*r1*x2 - 18*r1*y2 - 4*r2^2 + 3*x2^2 + 2*sqrt(3)*x2*y2 + y2^2
   81*b^2*r1^2 - 18*sqrt(3)*b^2*r1*x2 + 18*b^2*r1*y2 - 4*b^2*r2^2 + 3*b^2*x2^2 - 2*sqrt(3)*b^2*x2*y2 + b^2*y2^2 + 486*sqrt(3)*b*r1^3 - 324*b*r1^2*x2 - 12*sqrt(3)*b*r1*r2^2 + 18*sqrt(3)*b*r1*x2^2 - 6*sqrt(3)*b*r1*y2^2 + 2187*r1^4 - 486*sqrt(3)*r1^3*x2 - 486*r1^3*y2 - 108*r1^2*r2^2 + 81*r1^2*x2^2 + 54*sqrt(3)*r1^2*x2*y2 + 27*r1^2*y2^2
   252*b^2*r1^2 + 1332*sqrt(3)*b*r1^3 + 5184*r1^4
   x2^2 - (3*r1 + r2)^2 + (3*r1 - y2)^2

eq4 を解いて b を求める。

ans_b = solve(eq4, b)[2]
ans_b |> println

   -16*sqrt(3)*r1/7

eq3, eq5 の b に ans_b を代入する。

eq3 = 49eq3(b => ans_b)/r1^2 |> simplify
eq5 = eq5(b => ans_b) |> expand |> simplify;

eq1 |> println
eq3 |> println
eq5 |> println

   81*r1^2 - 18*sqrt(3)*r1*x2 - 18*r1*y2 - 4*r2^2 + 3*x2^2 + 2*sqrt(3)*x2*y2 + y2^2
   6075*r1^2 - 1350*sqrt(3)*r1*x2 - 9990*r1*y2 - 4332*r2^2 + 225*x2^2 + 1110*sqrt(3)*x2*y2 + 4107*y2^2
   -6*r1*r2 - 6*r1*y2 - r2^2 + x2^2 + y2^2

eq1, eq3, eq5 を解いて r1, x2, y2 を求める。

res2 = solve([eq1, eq3, eq5], (r1, x2, y2))[3];  # 3 of 3

res2[1] |> println
res2[2] |> println
res2[3] |> println

   2*r2*(sqrt(5) + 3)/9
   5*sqrt(3)*r2/6 + 2*sqrt(15)*r2/3
   3*r2/2

res2[1](r2 => 1/2).evalf() |> println
res2[2](r2 => 1/2).evalf() |> println
res2[3](r2 => 1/2).evalf() |> println

   0.581785330833310
   2.01268228522284
   0.750000000000000

大円の半径 r1 は，小円の半径 r2 の 2(√5 + 3)/9 倍である。
小円の直径が 1 寸のとき大円の直径は 2(√5 + 3)/9 = 1.16357066166662 寸である。
「術」の「小円径×(√20 + 6)/9」 と一致した。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   r1 = r2*2(√5 + 3)/9
   x2 = r2*(5√3/6 + 2√15/3)
   y2 = r2*3/2
   b = r1*(-16√3/7)
   a = r1*3√3
   x1 = -2√3*r1
   r0 = a/√3
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   circle(0, √3a/3, √3a/3)
   circle(x1, r1, r1, :green)
   circle(x2, y2, r2, :magenta)
   segment(a, 0, b, √3(b + a))
   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(0, r0, "全円:r0\n(0,r0)", :red, :center, :bottom, delta=delta/2)
       point(x1, r1, "大円:r1\n(x1,r1)", :green, :center, :bottom, delta=delta/2)
       point(x2, y2, "小円:r2\n(x2,y2)", :magenta, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, √3a, " √3a", :blue, :left, :vcenter)
   end
end;