算額（その1801） - 算額あれこれ

46 岩手県一関市舞川字龍ヶ沢観福寺観音堂明治34年(1901)

安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：3次元，球，回転楕円体，長立円
#Julia #SymPy #算額 #和算 #数学

長立円（楕円の長軸を中心とした回転楕円体）の周りに 6 個 2 組の等球が外接している。等球は隣同士外接している。
長立円の長径，短径が与えられたとき，等球の直径を得る術を述べよ。

長立円の長半径，短半径，中心座標を \(a, b, (0, 0, 0)\)
等円の半径と中心座標を \(r, (r, 0, y), -(r, 0, y), \dots\)
長立円と等球の接点座標を \( (x_0, 0, y_0)\)
とおき，以下の連立方程式を解く。

3元連立方程式であるが，SymPy では適切に解けない。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms a::positive, b::positive,
r::positive, y::positive,
x0::positive, y0::positive
r = y/2
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (x0 - r)/(y0 - y)
eq3 = (x0 - r)^2 + (y0 - y)^2 - r^2;
# res = solve([eq1, eq2, eq3], (y, x0, y0))

数値解

数値解は以下のようになる。

function H(u)
(y, x0, y0) = u
r = y/2
return [
x0^2/a^2 + y0^2/b^2 - 1, # eq1
-b^2*x0/(a^2*y0) + (x0 - r)/(y0 - y), # eq2
(x0 - r)^2 + (y0 - y)^2 - r^2, # eq3
]
end;

(a, b) = (5, 2)

iniv = BigFloat[3.75, 1.62, 1.89]
res = nls(H, ini=iniv)

([3.7494560410294966, 1.6202440532790994, 1.892080725880953], true)

y = 3.7494560410294966 なので，r = y/2 より等球の直径は 3.7494560410294966 である。

解析解

解析解は一度に得られないので，逐次求めていく。
まず eq1 を解いて \(y_0\) を求める。

@syms a::positive, b::positive,
r::positive, y::positive,
x0::positive, y0::positive
r = y/2
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (x0 - r)/(y0 - y)
eq3 = (x0 - r)^2 + (y0 - y)^2 - r^2;

# y0
ans_y0 = solve(eq1, y0)[1]

　\(\displaystyle \frac{b \sqrt{a^{2} - x_{0}^{2}}}{a}\)

ついで，eq2 に代入し，\(y\) を求める。

eq12 = eq2(y0 => ans_y0) |> simplify

　\(\displaystyle - \frac{x_{0} - \frac{y}{2}}{y - \frac{b \sqrt{a^{2} - x_{0}^{2}}}{a}} - \frac{b x_{0}}{a \sqrt{a^{2} - x_{0}^{2}}}\)

# y
ans_y = solve(eq12, y)[1]

　\(\displaystyle \frac{2 x_{0} \left(a^{2} - b^{2}\right) \sqrt{a^{2} - x_{0}^{2}}}{a \left(a \sqrt{a^{2} - x_{0}^{2}} - 2 b x_{0}\right)}\)

最後に eq3 に \(y, y_0\) を代入し \(x_0\) を求める。

eq13 = eq3(y => ans_y, y0 => ans_y0) |> expand |> simplify |> numerator

　\(\displaystyle - 4 a^{7} b x_{0} + a^{6} b^{2} \sqrt{a^{2} - x_{0}^{2}} + 3 a^{6} x_{0}^{2} \sqrt{a^{2} - x_{0}^{2}} - 2 a^{5} b^{3} x_{0} + 2 a^{5} b x_{0}^{3} + 8 a^{4} b^{2} x_{0}^{2} \sqrt{a^{2} - x_{0}^{2}} - 3 a^{4} x_{0}^{4} \sqrt{a^{2} - x_{0}^{2}} - 4 a^{3} b^{3} x_{0}^{3} + 2 a^{3} b x_{0}^{5} - 5 a^{2} b^{2} x_{0}^{4} \sqrt{a^{2} - x_{0}^{2}} - 2 a b^{3} x_{0}^{5} + 8 b^{4} x_{0}^{4} \sqrt{a^{2} - x_{0}^{2}}\)

# x0
ans_x0 = solve(eq13, x0)[1]

　\(\displaystyle \frac{a^{2}}{\sqrt{a^{2} + 4 b^{2}}}\)

しかし，得られた解は不適切解である。

ans_x0(a => 5, b=> 2).evalf() |> println

3.90434404721515

\(a=5, b=2\) のとき，eq13 を図に描くと以下のようになる。
解は 3 個あるのに，右端の解しか得られない。

using LaTeXStrings
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(eq13(a => 5, b => 2), xlims=(0, 4), xlabel=L"x_0", ylabel=L"eq13")
hline!([0])
vline!([0, 0.697663251751854, 1.6202440532790994, 3.90434404721515], color=:gray80)

res = solve(eq13(a => 5, b => 2), x0)

　\(\displaystyle \left[\begin{smallmatrix}\frac{25 \sqrt{41}}{41}\\5 \operatorname{CRootOf} {\left(127449 x^{8} - 331674 x^{6} + 235225 x^{4} - 25000 x^{2} + 400, 2\right)}\\5 \operatorname{CRootOf} {\left(127449 x^{8} - 331674 x^{6} + 235225 x^{4} - 25000 x^{2} + 400, 3\right)}\end{smallmatrix}\right]\)

res[1].evalf() |> println
res[2].evalf() |> println
res[3].evalf() |> println

3.90434404721515
0.697663251751854
1.62024405327910

\(x_0 = 1.62024405327910\)が適解である。
不適切な \(x)0\) からは，不適切な解しか得られない。

\(x_0 = 1.62024405327910, a = 5, b = 2\) のとき，
\(y = 2r = 3.74945604102950\) となり，等球の直径は 3.7494560410294966 である。

ans_y(x0 => 1.62024405327910, a => 5, b => 2) |> println

3.74945604102950

描画関数プログラムのソースを見る

function draw(a, b, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b) = (5, 2)
delta=a/50
p1 = plot(xlabel="x-axis", ylabel="z-axis")
(y, x0, y0) = [3.7494560410294966, 1.6202440532790994, 1.892080725880953]
r = y/2
ellipse(0, 0, a, b, color=:blue)
circle4(r, y, r)
circle4(r, r, r)
point(a, 0, "a", :blue, :left, :bottom, delta=delta/2)
point(0, b, "b", :blue, :center, :bottom, delta=delta/2)
point(x0, y0, "(x0,y0)", :green, :center, delta=-delta)
point(r, y, "(r,y)", :red, :center, delta=-delta)
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
p2 = plot(xlabel="y-axis", ylabel="z-axis")
circle(0, 0, y, :gray80)
circle(0, 0, b, :blue)
rotate(0, y, r, angle=60)
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
plot(p1, p2)
end;