多変数関数の積分である多重積分のうち，もっとも簡単な2変数関数の積分つまり2重積分についてまとめる。
2重積分は累次積分で
Maxima における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}
\begin{eqnarray}
\iint_D y \,dx dy &=& \int_{x_1}^{x_2} \left\{\int_{y_1}^{y_2} f(x, y) dy \right\} dx\\
&=& \cdots
\end{eqnarray}
f(x, y):= y$
x1: 0$
x2: 2$
y1: 0$
y2: 2$
I = integrate(
integrate(f(x,y), y, y1, y2),
x, x1, x2);
参考：領域 $D$ の塗りつぶし
ついでに，領域 $D$ の部分を塗りつぶして描いてみます。
f(x, y):= y$
x1: 0$
x2: 2$
y1: 0$
y2: 2$
draw2d(
title = "領域 D",
/* x 軸 y 軸の表示範囲の設定。*/
xrange = [-0.5, 2.5], xaxis = true, xlabel="x",
yrange = [-0.5, 2.5], yaxis = true, ylabel="y",
user_preamble="set grid front",
/* 縦横比 */
proportional_axes = xy,
/* 塗りつぶす色の指定。*/
fill_color = yellow,
/* 上の線。y = y2 を filled_func に代入。*/
filled_func = y2,
/* 下の線。y = y1 を x1 < x < x2 の範囲で描く。*/
explicit(y1, x, x1, x2)
)$
ついでに，3次元的なグラフを描いてみます。
f(x, y):= y$
x1: 0$
x2: 2$
y1: 0$
y2: 2$
draw3d(
/* 表示範囲 */
xrange = [-0.5, 2.5],
yrange = [-0.5, 2.5],
zrange = [0, 2.5],
/* 縦横比 */
proportional_axes = xyz, xyplane=0,
/* 領域 D*/
color=yellow,
parametric_surface(x, y, 0,
x, x1, x2, y, y1, y2),
/* z = f(x, y) */
color=blue,
explicit(f(x, y), x, x1, x2, y, y1, y2)
)$
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(x, y):= 1$
phi1(x):= -sqrt(1-x**2)$
phi2(x):= sqrt(1-x**2)$
x1: -1$
x2: 1$
I = integrate(
integrate(f(x, y), y, phi1(x), phi2(x)),
x, x1, x2);
参考：領域 $D$ の塗りつぶし
ついでに，領域 $D$ の部分を塗りつぶして描いてみます。
phi1(x):= -sqrt(1-x**2)$
phi2(x):= sqrt(1-x**2)$
x1: -1$
x2: 1$
draw2d(
title="領域 D",
/* x 軸 y 軸の表示範囲の設定。*/
xrange = [-1.5, 1.5], xaxis = true, xlabel="x",
yrange = [-1.5, 1.5], yaxis = true, ylabel="y",
proportional_axes = xy, /* 円をまんまるく表示したいときに */
user_preamble="set grid front",
fill_color = yellow,
/* 縁をなめらかにするには適宜大きい値を */
x_voxel = 50, y_voxel = 50,
/* 不等式で示された領域を塗りつぶす例 */
region(x**2 + y**2 < 1, x, -1.5, 1.5, y, -1.5, 1.5)
)$
ガウス積分の準備としての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}$$ をヤコビアンという。
Maxima でヤコビアンを計算する例：
x: r * cos(theta);
y: r * sin(theta);
jacobian([x, y], [r, theta]);
Maxima の jacobian() 関数は「関数行列」を与えるので，「関数行列式」にするには，さらに determinant() をとる。
determinant(jacobian([x, y], [r, theta]));
trigsimp(%);
ヤコビアンが $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
\end{eqnarray}
2*%pi* 'integrate(exp(-r**2) * r, r, 0, inf) =
2*%pi* integrate(exp(-r**2) * r, r, 0, inf);
/* 次に進む前に */
/* x や y の定義をなしにしておく。*/
kill(x, y)$
ガウス積分
$$ 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}$$
ということになる。Maxima がガウス積分 $I$ を知っていることを確認する。
'integrate(exp(-x**2), x, -inf, inf) =
integrate(exp(-x**2), x, -inf, inf);
被積分関数は偶関数であるから，積分範囲を $0$ からにすると答えは半分になるはずで，これも確認。
'integrate(exp(-x**2), x, 0, inf) =
integrate(exp(-x**2), x, 0, inf);
ガウス積分は，積分の上限が $\infty$ のときにだけ解析的に積分できる例。
積分範囲が $0$ から一般に $x$ までだと積分できない（結果が初等関数で表せない）。実際にやってみると…
assume(x > 0)$
'integrate(exp(-t**2), t, 0, x) =
integrate(exp(-t**2), t, 0, x);
参考：誤差関数
$\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$$
I(m):= integrate( x**m * exp(-x**2), x, 0, inf);
I(0);
I(2);
I(4);
I(6);
$m = 2 n$ ($n$ は整数) の場合を推測するために，$I(0)$ の値で割って，factor() 関数で分母分子を素因数分解してみる。
I(2)/I(0), factor;
I(4)/I(0), factor;
I(6)/I(0), factor;
ということは，
\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}$$
と推測される。実際確かめてみると…
/* m >= 2 の偶数に対して，上記の 2 n を m とおいて */
EvenI(m):= (m-1)!!/2**(m/2) * I(0);
/* たとえば m = 8 のとき */
EvenI(8) - I(8);
これで $m \geq 2$ が偶数の場合の $I(m)$ は求まった。
/* m = 2 から 10 まで MyEvenI(m) の値を表示させる。*/
for m: 2 thru 10 step 2
do(print("m = ", m, " のとき: ",
'integrate( x**m * exp(-x**2), x, 0, inf) = EvenI(m)))$
次に，$m$ が奇数（$m = 2 n +1$）の場合の $I(m)$ の一般形を予想してみる。
I(1);
I(3);
I(5);
I(7);
I(3)/I(1), factor;
I(5)/I(1), factor;
I(7)/I(1), factor;
以上のことから次のような推測ができる。
$$ I(2n+1) = n! \times I(1) $$