Return to 積分：いくつかの応用

参考：SymPy で楕円の面積・周，回転楕円体の表面積・体積

必要なモジュールの import

In [1]:
from sympy.abc import *
from sympy import *
# SymPy Plotting Backends (SPB) 
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

mathtext.fontset の設定

In [2]:
# グラフを描くためではなくデフォルト設定のため
import matplotlib.pyplot as plt

plt.rcParams['mathtext.fontset'] = 'cm'

楕円の面積

SymPy Plotting Backends で長半径 $a$，単半径 $b$ の楕円を描く例。（一部が微妙にカクカクするので，n を大きめに。）

In [3]:
a = 2
b = 1
p = plot_geometry(
    Ellipse(Point(0, 0), hradius=a, vradius=b),
    {"facecolor": "yellow", "edgecolor":"blue", "lw":2}, 
    aspect = 'equal', size=(6.4,6.4), n=500, show=False);

fig, ax = p.fig, p.ax
fig.tight_layout()

ax.grid(linestyle = 'dotted')
# 表示範囲
ax.set_xlim(-a-0.5, a+0.5)
ax.set_ylim(-b-0.5, b+0.5)
# 目盛設定
ax.set_xticks([-a, 0, a])
ax.set_xticklabels(labels = ["$-a$", "$0$", "$a$"], fontsize=14)
ax.set_yticks([-b, 0, b])
ax.set_yticklabels(labels = ["$-b$", "$0$", "$b$"], fontsize=14)

# x 軸，y 軸
ax.axhline(0, color='lightgray', dashes=(3, 3), linewidth=0.6)
ax.axvline(0, color='lightgray', dashes=(3, 3), linewidth=0.6);

$\displaystyle \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$ の楕円の面積。

$a > b > 0$ として $a$ は長半径，$b$ は短半径と呼ぶ。離心率 $e$ を使って書くと，$b = a \sqrt{1-e^2}$。

$\displaystyle \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$ より $\displaystyle y = \pm b \sqrt{1-\frac{x^2}{a^2}}$。上半分の面積，つまり $\displaystyle y = b \sqrt{1-\frac{x^2}{a^2}}$ と $x$ 軸で囲まれる部分の面積を求めて2倍すればよいから，

\begin{eqnarray}
S &=& 2 \int_{-a}^a b \sqrt{1 – \frac{x^2}{a^2}} dx \\
&&\qquad (x \equiv a \sin\theta, \ \sqrt{1 – \frac{x^2}{a^2}} = \cos\theta, \ dx = a \cos\theta d\theta)\\
&=& 2 a b \int_{-\pi/2}^{\pi/2} \cos^2\theta d\theta\\
&=& a b \int_{-\pi/2}^{\pi/2} (1 + \cos 2\theta) d\theta \\
&=&a b \left[ \theta + \frac{1}{2}\sin 2\theta\right]_{-\pi/2}^{\pi/2} \\
&=& \pi a b
\end{eqnarray}

SymPy で直接積分してみると，（$a>0, b>0$ を仮定して）

In [4]:
var('a b e', positive = True)
var('x y')

def f(x):
    return b*sqrt(1-x**2/a**2)

Eq(
    2*Integral(f(x), (x, -a, a)),
    2*integrate(f(x), (x, -a, a))
)
Out[4]:
$\displaystyle 2 \int\limits_{- a}^{a} b \sqrt{1 – \frac{x^{2}}{a^{2}}}\, dx = \pi a b$

楕円の周

In [5]:
a = 2
b = 1
p = plot_geometry(
    Ellipse(Point(0, 0), hradius=a, vradius=b),
    {"color":"blue", "lw":2}, is_filled = False,
    aspect = 'equal', size=(6.4,6.4), n=500, show=False);

fig, ax = p.fig, p.ax
fig.tight_layout()

ax.grid(linestyle = 'dotted')
# 表示範囲
ax.set_xlim(-a-0.5, a+0.5)
ax.set_ylim(-b-0.5, b+0.5)
# 目盛設定
ax.set_xticks([-a, 0, a])
ax.set_xticklabels(labels = ["$-a$", "$0$", "$a$"], fontsize=14)
ax.set_yticks([-b, 0, b])
ax.set_yticklabels(labels = ["$-b$", "$0$", "$b$"], fontsize=14)

# x 軸，y 軸
ax.axhline(0, color='lightgray', dashes=(3, 3), linewidth=0.6)
ax.axvline(0, color='lightgray', dashes=(3, 3), linewidth=0.6);

楕円の式を媒介変数表示すると，（$x = a \cos\theta, \ y = b \sin\theta$ とするところを趣向を変えて） $x = a \sin\theta, \ y = b \cos\theta$。$0 \le \theta \le \pi/2$ の長さを4倍すればよいから，

\begin{eqnarray}
L &=& 4 \int_0^{\frac{\pi}{2}} \sqrt{\left(\frac{dx}{d\theta}\right)^2 + \left(\frac{dy}{d\theta}\right)^2} \,d\theta\\
&=& 4 \int_0^{\frac{\pi}{2}} \sqrt{a^2 \cos^2\theta + b^2 \sin^2\theta} \,d\theta \\
&=& 4 a \int_0^{\frac{\pi}{2}} \sqrt{\cos^2\theta + (1-e^2) \sin^2\theta} \,d\theta \\
&=& 4 a \int_0^{\frac{\pi}{2}} \sqrt{1 -e^2 \sin^2\theta} \,d\theta \\
&=& 4 a E(e)
\end{eqnarray}

ここで
$\displaystyle E(k) \equiv \int_0^{\frac{\pi}{2}} \sqrt{1 -k^2 \sin^2\theta} \,d\theta$ は第二種完全楕円積分である。

要するに，楕円の周を求める積分は初等関数で表すことができない。

SymPy では，第二種完全楕円積分は elliptic_e(m) として組み込まれている。

In [6]:
def E(e):
    return elliptic_e(e**2)

離心率 $e$ は $0 \le e < 1$ であるから…

In [7]:
E(0)
Out[7]:
$\displaystyle \frac{\pi}{2}$
In [8]:
E(1)
Out[8]:
$\displaystyle 1$

$0 \le e \le 1$ の範囲で関数 $E(e)$ のグラフを描いておく。

In [9]:
# plot(E(e), (e, 0, 1)) と書きたいところだが...
# plot(E(x), (x, 0, 1)) として (x) を省略して...
plot(E, (x, 0, 1), xlabel="$e$", ylabel="$E(e)$");

In [10]:
plot(E(e), (e, 0, 1), xlabel="$e$", ylabel="$E(e)$");
# と書くと，以下のような UserWarning が最初に出ますが，
# 文句を言いつつもグラフを描いてくれます。

/usr/local/lib/python3.8/dist-packages/spb/series.py:228: UserWarning: The evaluation with NumPy/SciPy failed.
NameError: name 'elliptic_e' is not defined
Trying to evaluate the expression with Sympy, but it might be a slow operation.
  warnings.warn(

回転楕円体の表面積 1.

SymPy Plotting Backends で回転楕円体を描く例。

In [11]:
var('a b e', positive = True)
var('x y')

def f(x):
    return b*sqrt(1-x**2/a**2)
a = 2
b = 1

# 回転楕円体の表面
p1 = plot3d_parametric_surface(
    x, f(x)*cos(t), f(x)*sin(t), (x, -a, a), (t, 0, 2*pi),
    {'cmap':'Blues', "alpha":0.6}, 
    use_cm=True, color_func=-f(x)*cos(t), colorbar=False, 
    axis=False, size=(6.4, 6.4), n=60, 
    xlim=(-2.2,2.2), ylim=(-1.2, 1.2), zlim=(-1.2, 1.2), 
    aspect="equal", show=False
);

p2 = plot3d_parametric_line(
    # x 軸
    (x, 0, 0, (x, -2.5, a), "", {"color":"blue","zorder":0, "lw":2}), 
    (x, 0, 0, (x, a,  2.5), "", {"color":"blue","zorder":10,"lw":2}), 
    # いわゆる「y 軸」，実は z 軸
    (0, 0, z, (z, -1.5, b), "", {"color":"red","zorder":0, "lw":0.7}), 
    (0, 0, z, (z, b,  1.5), "", {"color":"red","zorder":10,"lw":0.7}), 
    # いわゆる「y = y(x)」，実は z = z(x)
    (x, 0, f(x), (x, -a, -a+0.1), "", {"color":"darkblue","zorder":0, "lw":2}),
    (x, 0, f(x), (x, -a+0.1,  a), "", {"color":"darkblue","zorder":10,"lw":2}),
    use_cm=False, show=False
)

(p1+p2).show();

$\displaystyle y = b \sqrt{1-\frac{x^2}{a^2}} = \sqrt{b^2 – (1-e^2) x^2}$ を $x$ 軸のまわりに回転してできる回転楕円体の表面積は

\begin{eqnarray}
S &=& \int_{-a}^a 2\pi y \sqrt{1+ \left(\frac{dy}{dx}\right)^2}\,dx \\
\end{eqnarray}

… ということで，$y$ を定義し，

$$\displaystyle f \equiv 2 \pi y \sqrt{1+\left(\frac{dy}{dx}\right)^2} = 2 \pi \sqrt{y^2 + y^2 \left(\frac{dy}{dx}\right)^2}$$

を計算すると…

In [12]:
var('a b e', positive = True)
var('x y')

y = b*sqrt(1 - (1-e**2)*x**2/b**2)

f = 2*pi*sqrt(y**2 + y**2*diff(y,x)**2)
f = simplify(f)
f
Out[12]:
$\displaystyle 2 \pi \sqrt{b^{2} + e^{4} x^{2} – e^{2} x^{2}}$

以上の結果をじっとみて，

$$f = 2 \pi b\sqrt{1 – \frac{e^2}{a^2} x^2}$$

となることがわかれば，以下のように積分してやると…

In [13]:
f = 2*pi*b*sqrt(1-e**2*x**2/a**2)

ans = integrate(f, (x, -a, a))
simplify(ans)
Out[13]:
$\displaystyle \frac{2 \pi a b \left(e \sqrt{1 – e^{2}} + \operatorname{asin}{\left(e \right)}\right)}{e}$

つまり，この回転楕円体の表面積 $S$ は

\begin{eqnarray}
S &=& 2 \pi ab\left(\sqrt{1-e^2} + \frac{\arcsin e}{e} \right) \\
&=& 2 \pi \left(b^2 + ab \,\frac{\arcsin e}{e} \right)
\end{eqnarray}

確かに定積分できることがわかったので，置換積分の立場から確認してみよう。そもそも計算するのは，

In [14]:
int1 = Integral(f, (x, -a, a))
int1
Out[14]:
$\displaystyle \int\limits_{- a}^{a} 2 \pi b \sqrt{1 – \frac{e^{2} x^{2}}{a^{2}}}\, dx$

$\displaystyle \frac{e x}{a} \rightarrow t$ と置換してやると…

In [15]:
int2 = int1.transform(e*x/a, t)
int2
Out[15]:
$\displaystyle \int\limits_{- e}^{e} \frac{2 \pi a b \sqrt{1 – t^{2}}}{e}\, dt$

つまり，基本的に $\displaystyle \int \sqrt{1-t^2} \, dt$ がわかればよい。これは以下のように積分できる。

In [16]:
Eq(Integral(sqrt(1-t**2), t), 
   integrate(sqrt(1-t**2), t)
  )
Out[16]:
$\displaystyle \int \sqrt{1 – t^{2}}\, dt = \frac{t \sqrt{1 – t^{2}}}{2} + \frac{\operatorname{asin}{\left(t \right)}}{2}$

SymPy では苦もなく出てくるが，人力では以下のように部分積分してみると…

\begin{eqnarray}
I &\equiv& \int \sqrt{1-t^2}\, dt \\
&=& t \sqrt{1-t^2} – \int t \, \frac{d}{dt} \sqrt{1-t^2}\,dt \\
&=& t \sqrt{1-t^2} + \int \frac{t^2}{\sqrt{1-t^2}}\,dt \\
&=& t \sqrt{1-t^2} + \int \frac{1}{\sqrt{1-t^2}} \,dt – \int \frac{1-t^2}{\sqrt{1-t^2}}\, dt \\
&=& t \sqrt{1-t^2} + \int \frac{1}{\sqrt{1-t^2}} \,dt – I \\
\therefore\ \ I &=& \frac{1}{2} \left(t \sqrt{1-t^2} + \int \frac{1}{\sqrt{1-t^2}} \,dt \right) \\
&=& \frac{t\sqrt{1-t^2}}{2} + \frac{\sin^{-1} t}{2}
\end{eqnarray}

定積分も以下のように…

In [17]:
int2.doit()
Out[17]:
$\displaystyle – \frac{2 \pi a b \left(- \frac{e \sqrt{1 – e^{2}}}{2} – \frac{\operatorname{asin}{\left(e \right)}}{2}\right)}{e} + \frac{2 \pi a b \left(\frac{e \sqrt{1 – e^{2}}}{2} + \frac{\operatorname{asin}{\left(e \right)}}{2}\right)}{e}$
In [18]:
simplify(_)
Out[18]:
$\displaystyle \frac{2 \pi a b \left(e \sqrt{1 – e^{2}} + \operatorname{asin}{\left(e \right)}\right)}{e}$

回転楕円体の表面積 2.

In [19]:
var('a b e', positive = True)
var('x y')

def g(z):
    return a*sqrt(1-z**2/b**2)
a = 2
b = 1

# 回転楕円体の表面
p1 = plot3d_parametric_surface(
    g(z)*cos(t), g(z)*sin(t), z, (z, -b, b), (t, 0, 2*pi),
    {'cmap':'Reds', 'alpha':0.5}, 
    use_cm=True, color_func=-g(z)*sin(t), colorbar=False, 
    axis=False, size=(6.4, 6.4), n=80, 
    xlim=(-2.2,2.2), ylim=(-2.2, 2.2), zlim=(-1.2,1.2), 
    aspect="equal", show=False
);

p2 = plot3d_parametric_line(
    # x 軸
    (x, 0, 0, (x, -2.5, a), "", {"color":"blue","zorder":0, "lw":0.7}), 
    (x, 0, 0, (x, a,  2.5), "", {"color":"blue","zorder":10, "lw":0.7}), 
    # いわゆる「y 軸」，実は z 軸
    (0, 0, z, (z, -2, b), "", {"color":"red","zorder":0, "lw":2}), 
    (0, 0, z, (z, b,  2), "", {"color":"red","zorder":10,"lw":2}), 
    # いわゆる「x = x(y)」，実は x = x(z)
    (g(z), 0, z, (z, -2, -0.5), "", {"color":"darkred","zorder":0,"lw":2}), 
    (g(z), 0, z, (z, -0.5,  b), "", {"color":"darkred","zorder":10,"lw":2}), 
    use_cm=False, show=False
)

(p1+p2).show();

$\displaystyle x = a \sqrt{1-\frac{y^2}{b^2}}= \frac{\sqrt{b^2-{y^2}}}{\sqrt{1-e^2}}$ を $y$ 軸のまわりに回転してできる回転楕円体の表面積は

\begin{eqnarray}
S &=& \int_{-b}^b 2\pi x \sqrt{1+ \left(\frac{dx}{dy}\right)^2}\,dy \\
\end{eqnarray}

In [20]:
var('a b e', positive = True)
var('x y')

x = a*sqrt(1 - y**2/b**2)

f = 2*pi*sqrt(x**2 + x**2*diff(x,y)**2)
f = simplify(f)
f
Out[20]:
$\displaystyle \frac{2 \pi a \sqrt{a^{2} y^{2} + b^{2} \left(b^{2} – y^{2}\right)}}{b^{2}}$

以上の結果をじっとみて，

$$f = 2 \pi a\sqrt{1 + \frac{a^2 e^2}{b^4} y^2}$$

となることがわかれば，以下のように積分してやると…

In [21]:
f = 2*pi*a*sqrt(1 + a**2 * e**2*y**2/b**4)

ans = integrate(f, (y, -b, b))
simplify(ans)
Out[21]:
$\displaystyle \frac{2 \pi \left(a e \sqrt{a^{2} e^{2} + b^{2}} + b^{2} \operatorname{asinh}{\left(\frac{a e}{b} \right)}\right)}{e}$

つまり，この回転楕円体の表面積 $S$ は

\begin{eqnarray}
S &=& 2 \pi \left(a \sqrt{a^2 e^2 + a^2 (1-e^2)} +
b^2 \frac{\operatorname{asinh} \frac{a e}{b}}{e} \right) \\
&=& 2 \pi \left(a^2 +
\frac{b^2}{e} \operatorname{asinh} \frac{e}{\sqrt{1-e^2}} \right) \\
&=& 2 \pi \left(a^2 +
b^2 \frac{\operatorname{atanh}e}{e} \right)
\end{eqnarray}

なお，$\operatorname{asinh}$ すなわち $\sinh^{-1}$ と $\operatorname{atanh}$ すなわち $\tanh^{-1}$ の関係については，

\begin{eqnarray}
y &\equiv& \operatorname{asinh} \frac{e}{\sqrt{1-e^2}} \\
\sinh y &=& \frac{e}{\sqrt{1-e^2}} \\
\cosh y &=& \sqrt{1 + \sinh^2 y} = \frac{1}{\sqrt{1-e^2}} \\
\therefore\ \ \tanh y &=& \frac{\sinh y}{\cosh y} = e \\
\therefore\ \ y &=& \operatorname{atanh} e = \operatorname{asinh} \frac{e}{\sqrt{1-e^2}}
\end{eqnarray}

確かに定積分できることがわかったので，置換積分の立場から確認してみよう。そもそも計算するのは，

In [22]:
var('a b e', positive = True)
var('x y')

int1 = Integral(f, (y, -b, b))
int1
Out[22]:
$\displaystyle \int\limits_{- b}^{b} 2 \pi a \sqrt{\frac{a^{2} e^{2} y^{2}}{b^{4}} + 1}\, dy$

$\displaystyle \frac{a e y}{b^2} \rightarrow t$ と置換してやると…

In [23]:
int2 = int1.transform(a*e*y/b**2, t)
int2
Out[23]:
$\displaystyle \int\limits_{- \frac{a e}{b}}^{\frac{a e}{b}} \frac{2 \pi b^{2} \sqrt{t^{2} + 1}}{e}\, dt$

つまり，基本的に $\displaystyle \int \sqrt{1+t^2} \, dt$ がわかればよい。これは以下のように積分できる。

In [24]:
Eq(Integral(sqrt(1+t**2), t), 
   integrate(sqrt(1+t**2), t)
  )
Out[24]:
$\displaystyle \int \sqrt{t^{2} + 1}\, dt = \frac{t \sqrt{t^{2} + 1}}{2} + \frac{\operatorname{asinh}{\left(t \right)}}{2}$

SymPy では苦もなく出てくるが，人力では以下のように部分積分してみると…

\begin{eqnarray}
I &\equiv& \int \sqrt{1+t^2}\, dt \\
&=& t \sqrt{1+t^2} – \int t \, \frac{d}{dt} \sqrt{1+t^2}\,dt \\
&=& t \sqrt{1-t^2} – \int \frac{t^2}{\sqrt{1+t^2}}\,dt \\
&=& t \sqrt{1-t^2} + \int \frac{1}{\sqrt{1+t^2}} \,dt – \int \frac{1+t^2}{\sqrt{1-t^2}}\, dt \\
&=& t \sqrt{1-t^2} + \int \frac{1}{\sqrt{1+t^2}} \,dt – I \\
\therefore\ \ I &=& \frac{1}{2} \left(t \sqrt{1+t^2} + \int \frac{1}{\sqrt{1+t^2}} \,dt \right) \\
&=& \frac{t\sqrt{1+t^2}}{2} + \frac{\sinh^{-1} t}{2}
\end{eqnarray}

定積分も以下のように…

In [25]:
int2.doit()
Out[25]:
$\displaystyle – \frac{2 \pi b^{2} \left(- \frac{a e \sqrt{\frac{a^{2} e^{2}}{b^{2}} + 1}}{2 b} – \frac{\operatorname{asinh}{\left(\frac{a e}{b} \right)}}{2}\right)}{e} + \frac{2 \pi b^{2} \left(\frac{a e \sqrt{\frac{a^{2} e^{2}}{b^{2}} + 1}}{2 b} + \frac{\operatorname{asinh}{\left(\frac{a e}{b} \right)}}{2}\right)}{e}$
In [26]:
simplify(_)
Out[26]:
$\displaystyle \frac{2 \pi \left(a e \sqrt{a^{2} e^{2} + b^{2}} + b^{2} \operatorname{asinh}{\left(\frac{a e}{b} \right)}\right)}{e}$

回転楕円体の体積 1.

$\displaystyle y = b \sqrt{1-\frac{x^2}{a^2}}$ を $x$ 軸のまわりに回転してできる回転楕円体の体積は

\begin{eqnarray}
V &=& \int_{-a}^a \pi y^2 \,dx \\
\end{eqnarray}

$\displaystyle y = b \sqrt{1-\frac{x^2}{a^2}}$ とおいて…

In [27]:
var('a b e', positive = True)
var('x y')

y = b*sqrt(1-x**2/a**2)
y
Out[27]:
$\displaystyle b \sqrt{1 – \frac{x^{2}}{a^{2}}}$
In [28]:
Eq(Integral(pi*y**2, (x, -a, a)),
   integrate(pi*y**2, (x, -a, a)))
Out[28]:
$\displaystyle \int\limits_{- a}^{a} \pi b^{2} \cdot \left(1 – \frac{x^{2}}{a^{2}}\right)\, dx = \frac{4 \pi a b^{2}}{3}$

回転楕円体の体積 2.

$\displaystyle x = a \sqrt{1-\frac{y^2}{b^2}}$ を $y$ 軸のまわりに回転してできる回転楕円体の体積は

\begin{eqnarray}
V &=& \int_{-b}^b \pi x^2 \,dy \\
\end{eqnarray}

$\displaystyle x = a \sqrt{1-\frac{y^2}{b^2}}$ とおいて…

In [29]:
var('a b e', positive = True)
var('x y')

x = a*sqrt(1-y**2/b**2)
x
Out[29]:
$\displaystyle a \sqrt{1 – \frac{y^{2}}{b^{2}}}$
In [30]:
Eq(Integral(pi*x**2, (y, -b, b)),
   integrate(pi*x**2, (y, -b, b)))
Out[30]:
$\displaystyle \int\limits_{- b}^{b} \pi a^{2} \cdot \left(1 – \frac{y^{2}}{b^{2}}\right)\, dy = \frac{4 \pi a^{2} b}{3}$

楕円および楕円体にまつわるエトセトラ

  • 楕円の周は第二種完全楕円積分で表される。（初等関数ではあらわせない。）
  • 回転楕円体の表面積を求める際には，以下のような簡単な無理関数の積分がでてくるので，しっかりやっておく。

$$\int \sqrt{1-t^2} \, dt, \quad \int \sqrt{1+t^2} \, dt$$

答えは逆三角関数や逆双曲線関数を含むからね。