算額（その1107）

算額（その1107）

二　岩手県花巻市大田 清水寺 嘉永三年(1850)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円6個，二等辺三角形

二等辺三角形の中に，中円 4 個，小円 1 個を容れる。中円 3 個に外接する大円の直径を求めよ。

二等辺三角形の底辺の長さを 2a，高さを h
大円の半径と中心座標を r0, (x0, r1*(1 + √3))
中円の半径と中心座標を r1, (r1, r1), (0, r1*(1 + √3)), (0, r1*(3 + √3))
小円の半径と中心座標を r2, (0, r1*(4 + √3) + r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive, r0::positive, x0::positive
eq1 = dist2(0, b, a, 0, 0, (4 + √Sym(3))r1 + r2, r2)
eq2 = dist2(0, b, a, 0, 0, (3 + √Sym(3))r1, r1)
eq3 = dist2(0, b, a, 0, r1, r1, r1)
eq4 = x0^2 + r1^2 - (r0 + r1)^2
eq5 = (x0 - r1)^2 + r1^2*(1 + √Sym(3))^2 - (r0 + r1)^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (a, b, r1, r0, x0) = u
   return [
       (a^2*b^2 - 8*a^2*b*r1 - 2*sqrt(3)*a^2*b*r1 - 2*a^2*b*r2 + 8*sqrt(3)*a^2*r1^2 + 19*a^2*r1^2 + 2*sqrt(3)*a^2*r1*r2 + 8*a^2*r1*r2 - b^2*r2^2)/(a^2 + b^2),  # eq1
       (a^2*b^2 - 6*a^2*b*r1 - 2*sqrt(3)*a^2*b*r1 + 6*sqrt(3)*a^2*r1^2 + 11*a^2*r1^2 - b^2*r1^2)/(a^2 + b^2),  # eq2
       a*b*(a*b - 2*a*r1 - 2*b*r1 + 2*r1^2),  # eq3
       r1^2 + x0^2 - (r0 + r1)^2,  # eq4
       r1^2*(1 + sqrt(3))^2 - (r0 + r1)^2 + (-r1 + x0)^2,  # eq5
   ]
end;
r2 = 1/2
iniv = BigFloat[1.96, 7.3, 0.85, 2.4, 3.17]
res = nls(H, ini=iniv)

   ([1.9558948090462631, 7.299498801620881, 0.8491981862085499, 2.431851652578137, 3.1692507766256446], true)

大円の半径は，小円の半径の 2.431851652578137/0.5 = 4.863703305156274 倍である。
小円の直径が 1 寸のとき，大円の直径は 4.863703305156274 寸である。

「術」では，「置四十八開平方加八個開平方加一個乗小円径」とあるが，これを山村は √48 + 1 倍とケアレスミスをしている。
術は，√(√48 + 8) + 1 倍といっており，これは二重根号を外して 1 + √2 + √6 となるが，いずれにしろそれは 4.863703305156273 倍ということで，上の数値解と一致している。

sqrt(√Sym(48) + 8)+1 |>sympy.sqrtdenest |> println
1 + sqrt(2) + sqrt(6) |> println

   1 + sqrt(2) + sqrt(6)
   4.863703305156273

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

   r2 = 0.5;  a = 1.95589;  b = 7.2995;  r1 = 0.849198;  r0 = 2.43185;  x0 = 3.16925

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, b, r1, r0, x0) = res[1]
   @printf("小円の直径が %g のとき，大円の直径は %g である。\n", 2r2, 2r0)
   @printf("r2 = %g;  a = %g;  b = %g;  r1 = %g;  r0 = %g;  x0 = %g\n", r2, a, b, r1, r0, x0)
   plot([a, 0, -a, a], [0, b, 0, 0], color=:magenta, lw=0.5)
   circle2(r1, r1, r1)
   circle(0, r1*(1 + √3), r1, :red)
   circle(0, r1*(3 + √3), r1, :red)
   circle(0, r1*(4 + √3)+ r2, r2, :blue)
   circle(x0, r1*(2 + √3), r0, :green)
   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(r1, r1, "中円:r1\n(r1,r1)", :red, :center, delta=-delta/2)
       point(0, r1*(1 + √3), "中円:r1\n(0,r1*(1+√3))", :red, :center, delta=-delta/2)
       point(0, r1*(3 + √3), "中円:r1\n(0,r1*(3+√3))", :red, :center, delta=-delta/2)
       point(0, r1*(4 + √3) + r2, "小円:r1\n(0,r1*(4+√3)+r2)", :red, :center, delta=-delta/2)
       point(x0, r1*(2 + √3), "大円:r0,(x0,r1*(2+√3))", :green, :center, delta=-delta/2)
       point(a, 0, "a", :magenta, :left, :bottom, delta=delta/2)
       point(0, b, " b", :magenta, :left, :vcenter)
   end
end;

算額（その1106）

算額（その1106）

百四 岩手県大船渡市田茂山 根城八幡宮 天保12年(1841)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円4個，斜線，直線上

直線上に 1 個の小円を挟んで 2 個の大円が載り，その3円に外接して大円 1 個が載っている。小円の直径が与えられたとき，大円の直径を求めよ。

大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (0, r2)
斜線の端点の座標を (0, 0), (x1, y1)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms r1::positive, x1::positive, y1::positive, r2::positive
eq1 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = dist2(0, 0, x1, y1, 0, 2r2 + r1, r1)
eq3 = dist2(0, 0, x1, y1, x1, r1, r1)
res = solve([eq1, eq2, eq3], (r1, x1, y1))[1]

   (r2*(-1 + sqrt(17))/2, r2*sqrt(-2 + 2*sqrt(17)), r2 + sqrt(17)*r2)

大円の半径 r1 は，小円の半径の (√17 - 1)/2 倍である。
「術」は「置四個二分五厘開平方内減八分余乗小円径」と書いているが，「八分」ではなく「五分」の誤記（山村の誤読？）である。「四個二分五厘」は 17/4 なので，√(17/4) - 1/2 = (√17 - 1)/2 である。

小円の直径が 1 寸のとき，大円の直径は (√17 - 1)/2 = 1.5615528128088303 寸である。

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

  r2 = 0.5;  r1 = 0.780776;  x1 = 1.24962;  y1 = 2.56155

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (r1, x1, y1) = (r2*(-1 + sqrt(17))/2, r2*sqrt(-2 + 2*sqrt(17)), r2 + sqrt(17)*r2)
   @printf("小円の直径が %g のとき，大円の直径は %g である。\n", 2r2, 2r1)
   @printf("r2 = %g;  r1 = %g;  x1 = %g;  y1 = %g\n", r2, r1, x1, y1)
   plot([-x1, 0, x1], [y1, 0, y1], color=:green, lw=0.5)
   circle(0, 2r2 + r1, r1)
   circle2(x1, r1, r1)
   circle(0, r2, r2, :blue)
   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, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(0, 2r2 + r1, "大円:r1,(0,2r2+r1)", :red, :center, delta=-delta/2)
       point(0, r2, "小円:r2\n(0,r2)", :blue, :center, :bottom, delta=delta/2)
       point(x1, y1, " (x1,y1)", :green, :left, :vcenter)
   end
end;

算額（その1105）

算額（その1105）

百四 岩手県大船渡市田茂山 根城八幡宮 天保12年(1841)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：楕円6個，正六角形

正六角形の中に等楕円を 6 個容れる。楕円の長径・短径が与えられたとき，正六角形の一辺の長さを求めよ。

計算式を簡単にするため，図形を時計方向に 30° 回転させたもので考える。x 軸上に長径を持つ楕円について以下のように記号を定める。
正六角形の一辺の長さを R
楕円の長半径，短半径，中心座標を a, b, (x0, 0)
隣の楕円との接点，正三角形の一辺との接点の座標を (x1, y1), (x2, y2); y1 = √3(R - x1), y2 = √3x2/3
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positive, b::positive, R::positive, x0::positive,
     x1::positive, y1::positive, x2::positive, y2::positive
y1 = √Sym(3)*(R - x1)
y2 = √Sym(3)*x2/3
eq1 = (x1 - x0)^2/a^2 + y1^2/b^2 - 1
eq2 = -b^2*(x1 - x0)/(a^2*y1) + √Sym(3)
eq3 = (x2 - x0)^2/a^2 + y2^2/b^2 - 1
eq4 = -b^2*(x2 - x0)/(a^2*y2) - 1/√Sym(3);
res = solve([eq1, eq2, eq3, eq4], (R, x0, x1, x2))[2]

   ((a^2 + 3*b^2 + sqrt(9*a^4 + 30*a^2*b^2 + 9*b^4)/3)/sqrt(a^2 + 3*b^2), sqrt(a^2 + 3*b^2), (a^2*sqrt(9*a^4 + 30*a^2*b^2 + 9*b^4) + (a^2 + 3*b^2)*(3*a^2 + b^2))/(sqrt(a^2 + 3*b^2)*(3*a^2 + b^2)), 3*b^2/sqrt(a^2 + 3*b^2))

各パラメータは以下のように簡約化して求めることができる。

t  = a^2 + 3b^2
u  = 3a^2 + b^2
x0 = sqrt(t)
R  = x0 + sqrt(u/3)
x1 = (a^2*sqrt(3t*u) + t*u)/(x0*u)
x2 = 3b^2/x0;

楕円の長径・短径が 10，7 のとき，正六角形の一辺の長さは sqrt((10/2)^2 + 3(7/2)^2) + sqrt((10/2)^2 + (7/2)^2/3) = 13.251013385205335 である。

sqrt((10/2)^2 + 3(7/2)^2) + sqrt((10/2)^2 + (7/2)^2/3) |> println

   13.251013385205335

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

   a = 5;  b = 3.5;  x0 = 7.85812;  R = 13.251;  x1 = 12.4938;  x2 = 4.67669

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (10, 7)./2
   t = a^2 + 3b^2
   u = 3a^2 + b^2
   x0 = sqrt(t)
   R = x0 + sqrt(u/3)
   R = sqrt(a^2 + 3b^2) + sqrt(a^2 + b^2/3)
   x1 = (a^2*sqrt(3t*u) + t*u)/(x0*u)
   x2 = 3b^2/x0
   y1 = √3(R - x1)
   y2 = √3x2/3
   @printf("楕円の長径・短径が %g, %g のとき，正六角形の一辺の長さは %g である。\n", 2a, 2b, R)
   @printf("a = %g;  b = %g;  x0 = %g;  R = %g;  x1 = %g;  x2 = %g\n", a, b, x0, R, x1, x2)
   θ = 0:60:420
   x = @. R*cosd(θ)
   y = @. R*sind(θ)
   plot()
   circle(0, 0, R)
   for i = 1:6
       ellipse(x0*cosd(60(i)), x0*sind(60(i)), a, b, φ=60(i), color=:green)
       segment(x[i], y[i], x[i + 1], y[i + 1], :blue)
   end
   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, y1, "(x1,y1) ", :green, :right, :vcenter)
       point(x2, y2, "(x2,y2)", :green, :left, delta=-delta/2)
       point(x0, 0, "楕円:a, b,(x0, 0)", :green, :center, delta=-delta/2)
   end
end;