算額（その731）

算額（その731）

福島県耶麻郡猪苗代町中小松西浜 猪苗代小平潟天満宮 明治14年(1881)
http://www.wasan.jp/fukusima/kohiragata3.html

街角の数学 Street Wasan　～落書き帳「○△□」～ 330．小平潟天満宮（その３） 明治14年
http://streetwasan.web.fc2.com/math17.11.1.html

等脚台形の中に大円が内接しており，大円に 3 個の楕円が内接している。楕円は隣の楕円と大円に一点で接する。台形の頂点近辺にある 6 個の小円の直径が 1 寸のとき，真ん中にある内円の直径はいかほどか。

「問」には「楕円は隣の楕円と大円に一点で接する」としか書いていないが，「大円と一点で接する『短径が最も長い楕円（曲率楕円）』」である。短径が短い（すなわち内円の直径が大きい）場合には解が定まらない。
『算法助術の公式86』で，直径 d の円に 1 点で内接する曲率楕円の長径 a と短径 b の関係が述べられている d = a^2/b

大円の中心を原点に置く。
大円の半径と中心座標を R, (0, 0)
台形の右上と右下の頂点の座標を (xb, R), (xa, -R)
上方の楕円の長径と短径と中心座標を a, b, (0, r1 + b)
内円の半径と中心座標を r1, (0, 0)
小円の半径と中心座標を r2, (x1, R - r2), (x2, y2), (x3, r2)
とおき以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R, a, b, r1, x0, y0, xa::positive, xb::positive, r2::positive, x1::positive, x2::positive, y2::negative, x3::positive
@syms d
eq1 = 2b + r1 - R
eq2 = R - a^2/b # b = (R - r1)/2
eq3 = x0^2/a^2 + (y0 - r1 - b)^2/b^2 - 1
eq4 = -b^2*x0/(a^2*(y0 - r1 - b)) - 1/sqrt(Sym(3))
eq5 = x0/y0 - sqrt(Sym(3));
eq6 = dist(xa, -R, xb, R, x1, R - r2) - r2^2
eq7 = dist(xa, -R, xb, R, x2, y2) - r2^2
eq8 = dist(0, 0, xa, - R, x2, y2) - r2^2
eq9 = x3^2 + (r2 - R)^2 - (R + r2)^2;
eq10 = x2^2 + y2^2 - (R + r2)^2
eq11 = xa*xb - R^2;
eq12 = x1^2 + (R - r2)^2 - (R + r2)^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)
   (R, r1, a, b, x0, y0, xa, xb, x1, x2, y2, x3) = u
   return [
       -R + 2*b + r1,  # eq1
       R - a^2/b,  # eq2
       -1 + (-b - r1 + y0)^2/b^2 + x0^2/a^2,  # eq3
       -sqrt(3)/3 - b^2*x0/(a^2*(-b - r1 + y0)),  # eq4
       x0/y0 - sqrt(3),  # eq5
       -r2^2 + (2*R - 2*R*(2*R*(2*R - r2) + (x1 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2) - r2)^2 + (x1 - xa - (-xa + xb)*(2*R*(2*R - r2) + (x1 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2))^2,  # eq6
       -r2^2 + (R - 2*R*(2*R*(R + y2) + (x2 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2) + y2)^2 + (x2 - xa - (-xa + xb)*(2*R*(R + y2) + (x2 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2))^2,  # eq7
       -r2^2 + (x2 - xa*(-R*y2 + x2*xa)/(R^2 + xa^2))^2 + (R*(-R*y2 + x2*xa)/(R^2 + xa^2) + y2)^2,  # eq8
       x3^2 + (-R + r2)^2 - (R + r2)^2,  # eq9
       x2^2 + y2^2 - (R + r2)^2,  # eq10
       -R^2 + xa*xb,  # eq11
       x1^2 + (R - r2)^2 - (R + r2)^2,  # eq12
   ]
end;

r2 = 1/2
iniv = BigFloat[7.1, 1.1, 4.6, 3.1, 3.2, 1.8, 8, 6.1, 5.3, 6.6, -4.5, 3]
res = nls(H, ini=iniv)

   ([3.5, 0.5, 2.29128784747792, 1.5, 1.5155444566227676, 0.875, 3.968626966596886, 3.086709862908689, 2.6457513110645907, 3.307189138830738, -2.25, 2.6457513110645907], true)

小円の直径が 1 寸のとき，内円の直径も 1 寸である。
なお，大円の直径は小円の直径の 7 倍，楕円の短径は小円の直径の 3 倍である。

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

R = 3.5;  r1 = 0.5;  a = 2.29129;  b = 1.5;  x0 = 1.51554;  y0 = 0.875;  xa = 3.96863;  xb = 3.08671;  x1 = 2.64575;  x2 = 3.30719;  y2 = -2.25;  x3 = 2.64575

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (R, r1, a, b, x0, y0, xa, xb, x1, x2, y2, x3) = res[1]
   @printf("R = %g;  r1 = %g;  a = %g;  b = %g;  x0 = %g;  y0 = %g;  xa = %g;  xb = %g;  x1 = %g;  x2 = %g;  y2 = %g;  x3 = %g\n", R, r1, a, b, x0, y0, xa, xb, x1, x2, y2, x3)
   plot()
   circle(0, 0, R)
   circle(0, 0, r1, :green)
   l = r1 + b
   ellipse(0, r1 + b, a, b, color=:blue)
   x = l*cos(pi/6)
   y = -l*sin(pi/6)
   ellipse(x, y, a, b, color=:blue, φ= 240)
   ellipse(-x, y, a, b, color=:blue, φ= 120)
   circle(x3, r2 - R, r2)
   circle(-x3, r2 - R, r2)
   segment(-xa, -R, xa, -R)
   segment(-xb, R, xb, R)
   segment(xa, -R, xb, R)
   segment(-xa, -R, -xb, R)
   circle(x1, R - r2, r2)
   circle(-x1, R - r2, r2)
   circle(x2, y2, r2)
   circle(-x2, y2, r2)
   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(x0, y0, "(x0,y0)", :blue, :left, delta=-delta/2)
       point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
       point(xa, -R, "(xa,-R)", :black, :right, :bottom, delta=delta/2)
       point(xb, R, "(xb,R)", :black, :right, :bottom, delta=delta/2)
       point(0, r1, " r1", :blue, :center, :bottom, delta=delta/2)
       point(0, R - b, "R-b", :blue, :center, :bottom, delta=delta/2)
       point(a, R - b, "(a,R-b) ", :blue, :right, :vcenter)
       point(x1, R - r2, "(x1,R-r2)", :red, :right, :bottom, delta=delta)
       point(x2, y2, "小円:r2,(x2,y2)", :red, :right, :vcenter)
       point(x3, r2 - R, " (x3,r2-R)", :red, :left, :vcenter)
       point(0, 0, "内円:r1", :green, :center, :vcenter, mark=false)
   end
end;

 

算額（その730）

算額（その730）

福島県福島市 『黒岩虚空蔵堂算額』第9問 明治26年
街角の数学 Street Wasan　～落書き帳「○△□」～ 347．楕円に内外接する円
http://streetwasan.web.fc2.com/math17.12.5.html

外円内に互いに接している楕円と円 2 個が入っている。楕円の長径と短径，乙円の直径が与えられたとき，甲円の直径を求めよ。

まさに「算法助術の公式85」が該当する問題であるが，図形を描くためのパラメータをすべて得るためには公式85は力不足である。

基本的には，以下の 6 元連立方程式を解けばよいのだが，SymPy の数式処理能力に限界があるので，逐次的に解いていく。

eq1 = (x1 - x2)^2 + y1^2 - (r1 + r2)^2
eq2 = x1^2 + y1^2 - (a - r1)^2
eq3 = (x0 - x1)^2 + (y0 - y1)^2 - r1^2
eq4 = (x0 - x2)^2 + y0^2 - r2^2
eq5 = x0^2/a^2 + y0^2/b^2 - 1
eq6 = -b^2*x0/(a^2*y0) - (x2 - x0)/y0

原点を中心とする楕円と，(x, 0) を中心とする円の接点を (x0, y0) とする。
x0，y0 は同時に求めるとかえってややこしい式になるので，x0 を求める。y0 は円周上の点であることから決定できる。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, y1::positive, r2::positive, x2::positive,
     a::positive, b::positive, x0::positive, y0::positive
a, b, r2 が既知
y0 = sqrt(r2^2 - (x0 - x2)^2)
eq4 = (x0 - x2)^2 + y0^2 - r2^2
eq6 = -b^2*x0/(a^2*y0) - (x2 - x0)/y0
res1 = solve(eq6, x0)

   1-element Vector{Sym{PyCall.PyObject}}:
    a^2*x2/(a^2 - b^2)

x0, y0 がわかったとして(x2 を含む式であるが），これを方程式に組み込み，x2 を求める。

@syms r1::positive, x1::positive, y1::positive, r2::positive, x2::positive,
     a::positive, b::positive, x0::positive, y0::positive

x0 = a^2*x2/(a^2 - b^2)
y0 = sqrt(r2^2 - (x0 - x2)^2)
eq5 = x0^2/a^2 + y0^2/b^2 - 1
res2 = solve(eq5, x2)

   1-element Vector{Sym{PyCall.PyObject}}:
    sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)/b

更に x2 を方程式に組み込み，r1, x1, y1 を求める。

@syms r1::positive, x1::positive, y1::positive, r2::positive, x2::positive,
     a::positive, b::positive, x0::positive, y0::positive

x2 = sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)/b
x0 = a^2*x2/(a^2 - b^2)
y0 = sqrt(r2^2 - (x0 - x2)^2)
eq1 = (x1 - x2)^2 + y1^2 - (r1 + r2)^2
eq2 = x1^2 + y1^2 - (a - r1)^2
eq3 = (x0 - x1)^2 + (y0 - y1)^2 - r1^2;
res3 = solve([eq1, eq2, eq3], (r1, x1, y1))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (r2*(a*r2 - b^2)/(2*b^2), (2*a^2*b^2 - 2*a^2*r2^2 + a*b^2*r2 - a*r2^3 - b^4 + b^2*r2^2)/(2*b*sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)), sqrt((a^6*r2^2 - a^4*b^4 - 2*a^4*b^2*r2^2 + 2*a^2*b^6 + a^2*b^4*r2^2 - b^8)/(a^6 - 3*a^4*b^2 + 3*a^2*b^4 - b^6))*(a*r2 + b^2)/(2*b^2))

以上をまとめると，
x2 = sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)/b
x0 = a^2*x2/(a^2 - b^2)
y0 = sqrt(r2^2 - (x0 - x2)^2)
r1 = r2*(a*r2 - b^2)/(2b^2)
x1 = (2a^2*b^2 - 2a^2*r2^2 + a*b^2*r2 - a*r2^3 - b^4 + b^2*r2^2)/(2b*sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2))
y1 = sqrt((a^6*r2^2 - a^4*b^4 - 2a^4*b^2*r2^2 + 2a^2*b^6 + a^2*b^4*r2^2 - b^8)/(a^6 - 3a^4*b^2 + 3a^2*b^4 - b^6))*(a*r2 + b^2)/(2b^2)
となる。

たとえば，楕円の長径と短径が 100, 50，乙円の直径が 40 のとき，甲円の直径は 12 である。

(a, b, r2) = (50, 25, 20)
r1 = r2*(a*r2 - b^2)/(2*b^2)
println("2r1 = $(2r1)")

   2r1 = 12.0

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

   a = 50;  b = 25;  r2 = 20
   r1 = 6;  x1 = 37.2391;  y1 = 23.4361;  x2 = 25.9808;  x0 = 34.641;  y0 = 18.0278

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r2) = (100, 50, 40) ./ 2
   x2 = sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)/b
   x0 = a^2*x2/(a^2 - b^2)
   y0 = sqrt(r2^2 - (x0 - x2)^2)
   r1 = r2*(a*r2 - b^2)/(2b^2)
   x1 = (2a^2*b^2 - 2a^2*r2^2 + a*b^2*r2 - a*r2^3 - b^4 + b^2*r2^2)/(2b*sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2))
   y1 = sqrt((a^6*r2^2 - a^4*b^4 - 2a^4*b^2*r2^2 + 2a^2*b^6 + a^2*b^4*r2^2 - b^8)/(a^6 - 3a^4*b^2 + 3a^2*b^4 - b^6))*(a*r2 + b^2)/(2b^2)
   @printf("楕円の長径 = %g;  楕円の短径 = %g;  乙円の直径 = %g;  甲円の直径 = %g\n", 2a, 2b, 2r2, 2r1)
   @printf("a = %g;  b = %g;  r2 = %g\n", a, b, r2)
   @printf("r1 = %g;  x1 = %g;  y1 = %g;  x2 = %g;  x0 = %g;  y0 = %g\n", r1, x1, y1, x2, x0, y0)
   plot()
   circle(0, 0, a)
   ellipse(0, 0, a, b, color=:blue)
   circle4(x1, y1, r1, :orange)
   circle(x2, 0, r2, :green)
   circle(-x2, 0, r2, :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(a, 0, " a", :blue, :left, :vcenter)
       point(0, b, "b", :blue, :center, :bottom, delta=delta/2)
       point(x2, 0, "乙円:r2,(x2,0)", :green, :center, delta=-delta/2)
       point(x0, y0, "(x0,y0) ", :black, :right, delta=-delta/2)
       point(x1, y1, "甲円:r1,(x1,y1) ", :black, :right, :bottom, delta=delta/2)
   end
end;