算額（その1113）

算額（その1113）

百一 岩手県大船渡市猪川町 田茂山町神明杜 田茂山神社（奉納年不明）
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円6個，円弧，正三角形

円弧の中に正三角形 1 個，甲円 2 個，乙円 4 個を容れる。甲円の直径が 1 寸のとき，弦の長さを求めよ。

正三角形の一辺の長さを 2a
円弧の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, y0 + r1)
乙円の半径と中心座標を r2, (x21, y0 + r2), (x22, y22)
弦と y 軸の交点座標を y0
とおき，以下の連立方程式を解く。
弦の長さは 2sqrt(R^2 - y0^2) である。
甲円の中心と正三角形の底辺の右側の頂点を結ぶ直線は円弧の中心を通り，x 軸となす角は 60° なので，弦と y 軸の交点座標 y0 = R/2 である。

include("julia-source.txt")

using SymPy
@syms a::positice, y0::positive, R::positive,
     r1::positive, x1::positive,
     r2::positive, x21::positive, x22::positive, y22::positive
x1 = r1 /√Sym(3)+ a
y0 = √Sym(3)a
R = 2y0
eq1 = x1^2 + (y0 + r1)^2 - (R - r1)^2 |> expand
eq2 = x21^2 + (y0 + r2)^2 - (R - r2)^2 |> expand
eq3 = x22^2 + y22^2 - (R - r2)^2 |> expand
eq4 = (x21 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq5 = dist2(a, y0, 0, R, x22, y22, r2);

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")

   -8*a^2 + 20*sqrt(3)*a*r1/3 + r1^2/3,  # eq1
   -9*a^2 + 6*sqrt(3)*a*r2 + x21^2,  # eq2
   -12*a^2 + 4*sqrt(3)*a*r2 - r2^2 + x22^2 + y22^2,  # eq3
   a^2 + 2*sqrt(3)*a*r1/3 - 2*a*x21 + r1^2/3 - 4*r1*r2 - 2*sqrt(3)*r1*x21/3 + x21^2,  # eq4
   3*a^2 - 3*a*x22 - sqrt(3)*a*y22 - r2^2 + 3*x22^2/4 + sqrt(3)*x22*y22/2 + y22^2/4,  # eq5

eq1 を a について解くことができる。その 2√3 倍が弦の長さである。

ans_a = solve(eq1, a)[2]
ans_a |> println

   r1*(5*sqrt(3) + 9)/12

y0 = √Sym(3)ans_a
R = 2y0
x0 = sqrt(R^2 - y0^2)
弦 = 2x0 |> factor
#= 弦 =# 弦 |> println
弦(r1 => 1/2).evalf() |> println

   r1*(5*sqrt(3) + 9)/2
   4.41506350946110

弦の長さは，甲円の半径の (5√3 + 9)/2，甲円の直径の (5√3 + 9)/2 = (√75 + 9)/4 = 4.415063509461097 倍である。
これは「術」と一致する。

---

以下は，図を描くためにほかのパラメータを求める手順である。

まず，eq2, eq3, eq4, eq5 の a に上で求めた ans_a を代入する。

eq12 = eq2(a => ans_a) |> simplify
eq13 = eq3(a => ans_a) |> simplify
eq14 = eq4(a => ans_a) |> simplify
eq15 = eq5(a => ans_a) |> simplify
eq12 |> println
eq13 |> println
eq14 |> println
eq15 |> println

   -39*r1^2/4 - 45*sqrt(3)*r1^2/8 + 15*r1*r2/2 + 9*sqrt(3)*r1*r2/2 + x21^2
   -13*r1^2 - 15*sqrt(3)*r1^2/2 + 5*r1*r2 + 3*sqrt(3)*r1*r2 - r2^2 + x22^2 + y22^2
   9*sqrt(3)*r1^2/8 + 9*r1^2/4 - 4*r1*r2 - 3*sqrt(3)*r1*x21/2 - 3*r1*x21/2 + x21^2
   r1^2*(5*sqrt(3) + 9)^2/48 - r1*x22*(5*sqrt(3) + 9)/4 - sqrt(3)*r1*y22*(5*sqrt(3) + 9)/12 - r2^2 + 3*x22^2/4 + sqrt(3)*x22*y22/2 + y22^2/4

eq14, eq15 を解いて x21, x22 を求める。

(ans_x21, ans_x22) = solve([eq14, eq15], (x21, x22))[4]
#= x21 =# ans_x21 |> println
#= x22 =# ans_x22 |> println

   2*sqrt(r1)*sqrt(r2) + 3*r1*(1 + sqrt(3))/4
   5*sqrt(3)*r1/6 + 3*r1/2 + 2*sqrt(3)*r2/3 - sqrt(3)*y22/3

eq12 に ans_x21, ans_x22 を代入し，r2 を求める。

eq22 = eq12(x21 => ans_x21, x22 => ans_x22) |> simplify
ans_r2 = solve(eq22, r2)[1] |> simplify |> sympy.sqrtdenest |> simplify
#= r2 =# ans_r2 |> println
ans_r2(r1 => 1/2).evalf() |> println

   3*r1*(4*sqrt(3) + 13)/121
   0.247043841697630

eq13 に ans_x21, ans_x22, ans_r2 を代入し y22 を求める。

eq23 = eq13(x21 => ans_x21, x22 => ans_x22) |> simplify
eq33 = eq23(r2 => ans_r2) |> simplify
ans_y22 = solve(eq33, y22)[1] |> sympy.sqrtdenest |> simplify
#= y22 =# ans_y22 |> println
ans_y22(r1 => 1/2).evalf() |> println

   7*r1*(151 + 93*sqrt(3))/484
   2.25678210302411

以上をまとめると,r1 が与えられたとき各パラメータの値を得ることができる。

r1 = 1/2
a = r1*(5√3 + 9)/12
r2 = 3*r1*(4√3 + 13)/121
y22 = 7*r1*(151 + 93√3)/484
x21 = 2*sqrt(r1)*sqrt(r2) + 3r1*(1 + √3)/4
x22 = 5√3*r1/6 + 3*r1/2 + 2√3*r2/3 - √3y22/3
(a, r2, x21, x22, y22)

その他のパラメータは以下の通りである。

   r1 = 0.5;  R = 2.54904;  x0 = 2.20753;  a = 0.735844;  y0 = 1.27452;  x1 = 1.02452;  r2 = 0.247044;  x21 = 1.72743;  x22 = 0.453996;  y22 = 2.25678

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   a = r1*(5√3 + 9)/12
   r2 = 3*r1*(4√3 + 13)/121
   y22 = 7*r1*(151 + 93√3)/484
   x21 = 2*sqrt(r1)*sqrt(r2) + 3r1*(1 + √3)/4
   x22 = 5√3*r1/6 + 3*r1/2 + 2√3*r2/3 - √3y22/3
   x1 = r1 /√3 + a
   y0 = √3a
   R = 2y0
   x0 = sqrt(R^2 - y0^2)
   θ = atand(y0, x0)
   @printf("甲円の直径が %g のとき，弦の長さは %g である。\n", 2r1, 2x0)
   @printf("r1 = %g;  R = %g;  x0 = %g;  a = %g;  y0 = %g;  x1 = %g;  r2 = %g;  x21 = %g;  x22 = %g;  y22 = %g\n", r1, R, x0, a, y0, x1, r2, x21, x22, y22)
   plot([a, 0, -a, a], y0 .+ [0, √3a, 0, 0], color=:green, lw=0.5)
   circle(0, 0, R, :green, beginangle=θ, endangle=180 - θ)
   circle2(x1, y0 + r1, r1, :blue)
   circle2(x21, y0 + r2, r2, :red)
   circle2(x22, y22, r2, :red)
   segment(-x0, y0, x0, y0)
   segment(0, 0, x0, y0, :gray70)
   segment(0, 0, -x0, y0, :gray70)
   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, y0 + r1, "甲円:r1\n(x1,y0+r1)", :blue, :center, :bottom, delta=delta/2)
       point(x21, y0 + r2, "乙円:r2,(x21,y0+r2) ", :black, :right, :vcenter)
       point(x22, y22, " 乙円:r2,(x22,y22)", :black, :left, :vcenter)
       point(0, R, " R=y0+√3a", :black, :center, :bottom, delta=delta/2)
       point(0, y0, "y0", :black, :center, delta=-delta/2)
       point(a, y0, "(a,y0)", :black, :center, delta=-delta/2)
       point(x0, y0, "(x0,y0)", :black, :right, delta=-delta/2)
   end
end;