必要なパッケージを import します。
from sympy.abc import *
from sympy import *
# SymPy Plotting Backends (SPB): グラフを描く際に利用
from spb import *
# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']
多変数関数の積分である多重積分のうち，もっとも簡単な2変数関数の積分つまり2重積分についてまとめる。
2重積分は累次積分で
2重積分の計算は，基本的に累次積分でやります。
$\displaystyle \iint_D f(x, y)\, dx dy$ の領域 $D$ が何と何で囲まれているかを明らかにして積分する。
累次積分 1
領域 $D$ が $x_1 \leq x \leq x_2,\ y_1 \leq y \leq y_2$ である場合：
$$\iint_D f(x, y) \,dx\,dy = \int_{x_1}^{x_2} \left\{\int_{y_1}^{y_2} f(x, y)\,dy \right\} dx$$
または
$$\iint_D f(x, y) \,dx\,dy = \int_{y_1}^{y_2} \left\{\int_{x_1}^{x_2} f(x, y)\,dx \right\} dy$$
累次積分 2
領域 $D$ が $x_1 \leq x \leq x_1,\ \phi_1(x) \leq y \leq \phi_2(x)$ である場合：
$$ \iint_D f(x, y) \,dx\,dy = \int_{x_1}^{x_2} \left\{\int_{\phi_1(x)}^{\phi_2(x)} f(x, y)\,dy \right\} dx$$
累次積分 3
領域 $D$ が $\psi_1(y)\leq x \leq \psi_2(y),\ y_1 \leq y \leq y_2$ である場合：
$$ \iint_D f(x, y) \,dx\,dy = \int_{y_1}^{y_2} \left\{\int_{\psi_1(y)}^{\psi_2(y)} f(x, y)\,dx \right\} dy$$
例題
$$I = \iint_D y \,dx dy, \quad D: 0\leq x\leq 2, \ 0 \leq y \leq 2$$
これは，ケース 1。以下のように累次積分にして…
\begin{eqnarray}
f(x, y) &=& y \\
x_1 &=& 0\\
x_2 &=& 2\\
y_1 &=& 0\\
y_2 &=& 2
\end{eqnarray}
SymPy では2重積分は以下のように。
$\displaystyle \int_{y_1}^{y_2} \left\{\int_{x_1}^{x_2} f(x, y) dx \right\} dy =
$
integrate(f(x, y), (x, x1, x2), (y, y1, y2))
f = y
x1 = 0
x2 = 2
y1 = 0
y2 = 2
# 手っ取り早く答えだけを表示させたいなら
integrate(f, (x, x1, x2), (y, y1, y2))
# 何を計算するか表示させてから答えを表示させるなら
Eq(Integral(y, (x, x1, x2), (y, y1, y2)),
Integral(y, (x, x1, x2), (y, y1, y2)).doit())
参考：領域 $D$ の塗りつぶし
ついでに，領域 $D$ の部分を塗りつぶして描いてみます。
$x_1 \leq x \leq x_2$ の範囲で $y = y1$ と $y = y2$ に囲まれた領域を塗りつぶす方法は，SymPy Plotting Backends では以下のようにします。
# NumPy に頼らない例
# x1 ~ x2 を 100 等分した数値リストを作成
x_list = [(x1 + (x2-x1)*i/100) for i in range(101)]
# x_array に対応した y1 の数値リストを作成
y1_list = lambdify(x, y1)(x_list)
# x_array に対応した y2 の数値リストを作成
y2_list = lambdify(x, y2)(x_list)
plot(xlim = (-0.5, 2.5),
ylim = (-0.5, 2.5),
aspect = "equal",
title = "領域 $D$",
xlabel = "x", ylabel = "y",
fill={'x': x_list,
'y1':y1_list,
'y2':y2_list,
'color':'yellow'});
ついでに，3次元的なグラフを描いてみます。
plot3d((0, (x, x1, x2), (y, y1, y2)),
(f, (x, x1, x2), (y, y1, y2)),
xlim = (-0.5, 2.5), ylim = (-0.5, 2.5),
zlim = (0, 3), n = 10,
rendering_kw=[{"color":"yellow"},
{"color":"lightblue"}]);
2重積分の結果は，黄色で塗りつぶされた領域 $D$ と薄青色で塗られた平面で挟まれた部分の体積を表す。
縦横高さともに $2$ の立方体を斜めに2つに切った立体であるので体積は
$$\frac{2\cdot 2\cdot 2}{2} = 4$$
になるはずで，答えも確かにそうなっている。
例題
$$I = \iint_D dx dy, \quad D: x^2 + y^2 \leq 1$$
累次積分の形にする。まず先に $-\sqrt{1-x^2} < y < \sqrt{1-x^2}$ の範囲で $y$ で定積分し，そのあとに $ -1 < x < 1$ の範囲で $x$ で定積分する。
被積分関数は $1$。
ケース 2 の累次積分の形：
\begin{eqnarray}
f(x, y) &=& 1 \\
x_1 &=& -1\\
x_2 &=& 1 \\
\phi_1(x) &=& -\sqrt{1-x^2}\\
\phi_2(x) &=& \sqrt{1-x^2}
\end{eqnarray}
\begin{eqnarray}
\iint_D dx dy &=& \int_{x_1}^{x_2} \left\{\int_{\phi_1(x)}^{\phi_2(x) } f(x, y)\,dy \right\} dx \\
&=& \cdots
\end{eqnarray}
$I$ は半径 $1$ の円の面積であるから，$\pi r^2 = \pi \ (\because r=1)$ となるはず。
f = 1
x1 = -1
x2 = 1
phi1 = -sqrt(1-x**2)
phi2 = sqrt(1-x**2)
Eq(Integral(f, (y, phi1, phi2), (x, x1, x2)),
Integral(f, (y, phi1, phi2), (x, x1, x2)).doit())
参考：領域 $D$ の塗りつぶし
ついでに，領域 $D$ の部分を塗りつぶして描いてみます。
# NumPy を import して使う場合
import numpy as np
x_array = np.linspace(x1, x2, 100)
def y1(x):
return -np.sqrt(1-x**2)
def y2(x):
return np.sqrt(1-x**2)
y1_array = y1(x_array)
y2_array = y2(x_array)
plot(xlim = (-1.5, 1.5),
ylim = (-1.5, 1.5),
aspect = "equal",
title = "領域 $D$",
xlabel = "x", ylabel = "y",
fill={'x': x_array,
'y1':y1_array,
'y2':y2_array,
'color':'yellow'});
ガウス積分の準備としての2重積分
$$ S = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-x^2 – y^2} dx\,dy$$
極座標に変換して積分
極座標に変換して計算する。
デカルト座標 $x, y$ と2次元極座標 $r, \theta$ との間の関係は
\begin{eqnarray}
x &=& r\cos\theta\\
y &=& r\sin\theta
\end{eqnarray}
であり，微小面積要素 $dx\,dy$ は以下のように変換される。
\begin{eqnarray}
dx\,dy &=& \frac{\partial(x, y)}{\partial(r,\theta)} \,dr\, d\theta
\end{eqnarray}
ここで，ここで，$$ \frac{\partial(x,y)}{\partial(r,\theta)} \equiv
\begin{vmatrix}
\frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta}\\
\frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta}\\
\end{vmatrix}
= \frac{\partial x}{\partial r} \frac{\partial y}{\partial \theta} – \frac{\partial y}{\partial r}\frac{\partial x}{\partial \theta}$$ をヤコビアンという。
ヤコビ行例
SymPy でヤコビアンを計算する例：
X = Matrix([r*cos(theta), r*sin(theta)])
Y = Matrix([r, theta])
jaco = X.jacobian(Y)
jaco
ヤコビアン（ヤコビ行列式）
jacobian() はヤコビ行列。この行列の行列式が，我々が求めたいヤコビアンはヤコビ行列の行列式。det() で行列式。
det(jaco).simplify()
ヤコビアンが $r$ となることがわかった。
\begin{eqnarray}
\therefore\ \ S &=& \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-x^2 – y^2} dx\,dy \\
&=& \int_0^{2\pi} \int_0^{\infty} e^{-r^2} r \,dr \,d\theta\\
&=& 2 \pi \int_0^{\infty} e^{-r^2} r \,dr = \cdots
\end{eqnarray}
Eq(Integral(exp(-r**2)*r, (theta, 0, 2*pi), (r, 0, oo)),
Integral(exp(-r**2)*r, (theta, 0, 2*pi), (r, 0, oo)).doit())
ガウス積分
$$ I = \int_{-\infty}^{\infty} e^{-x^2 } dx$$
とすると，上で
\begin{eqnarray}
S &=& \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-x^2 – y^2} dx\,dy \\
&=& I^2 \\
&=& \pi
\end{eqnarray}
ということがわかったので，
$$ I = \int_{-\infty}^{\infty} e^{-x^2 } dx = \sqrt{S} = \sqrt{\pi}$$
ということになる。SymPy がガウス積分 $I$ を知っていることを確認する。
Eq(Integral(exp(-x**2), (x, -oo, oo)),
Integral(exp(-x**2), (x, -oo, oo)).doit())
被積分関数は偶関数であるから，積分範囲を $0$ からにすると答えは半分になるはずで，これも確認。
Eq(Integral(exp(-x**2), (x, 0, oo)),
Integral(exp(-x**2), (x, 0, oo)).doit())
ガウス積分は，積分の上限が $\infty$ のときにだけ解析的に積分できる例。
積分範囲が $0$ から一般に $x$ までだと積分できない（結果が初等関数で表せない）。実際にやってみると…
Eq(Integral(exp(-t**2), (t, 0, x)),
Integral(exp(-t**2), (t, 0, x)).doit())
参考：誤差関数
$\operatorname{erf} (x)$ は「誤差関数」と呼ばれ，これが定義：
\begin{eqnarray}
\operatorname{erf} (x) &\equiv& \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} \, dt
\end{eqnarray}
ガウス積分に関連した積分
$$I(m) \equiv \int_{0}^{\infty} x^m e^{-x^2 } dx$$
def I(m):
return integrate(x**m * exp(-x**2), (x, 0, oo))
I(0)
I(2)
I(4)
I(6)
$m = 2 n$ ($n$ は整数) の場合を推測するために，$I(0)$ の値で割ってみる。
I(2)/I(0)
I(4)/I(0)
I(6)/I(0)
ということは，
\begin{eqnarray}
I(2) &=& 1\times \frac{1}{2} \times I(0), \\
I(4) &=& \frac{3\cdot 1}{2^2} \times I(0), \\
I(6) &=& \frac{5\cdot3\cdot 1}{2^3} \times I(0), \dots
\end{eqnarray}
だから，
$$\displaystyle I(2 n) = \frac{(2n-1)!!}{2^{n}} \frac{\sqrt{\pi}}{2}$$
と推測される。実際確かめてみると…
（SymPy では $n! = $ factorial(n)，$n!! = $ factorial2(n)）
# m >= 2 の偶数に対して，上記の 2n を m とおいて
def EvenI(m):
return factorial2(m-1)/2**(m/2) * I(0)
# 例えば m = 8 のとき
EvenI(8) - I(8)
これで $m \geq 2$ が偶数の場合の $I(m)$ は求まった。
次に，$m$ が奇数（$m = 2 n +1$）の場合の $I(m)$ の一般形を予想してみる。
I(1)
I(3)
I(5)
I(7)
I(3)/I(1)
I(5)/I(1)
I(7)/I(1)
以上のことから次のような推測ができる。
$$ I(2n+1) = n! \times I(1) $$