1変数関数の定積分 $\int_a^b f(x)\,dx$ は「面積」を求める操作でした。これを2変数、3変数に拡張したものが重積分です。二重積分は曲面の下の「体積」を、三重積分はさらに高次元の領域での「質量」や「総量」を計算します。
重積分は物理学のあらゆる場面で登場します。質量分布の計算、慣性モーメント、電荷分布、流体の流量など、空間に広がった量を扱うには重積分が不可欠です。
本記事の内容
- 二重積分の定義
- 逐次積分（累次積分）による計算
- 積分順序の交換
- 三重積分
- 体積・質量・重心の計算
- Pythonでの数値積分と可視化
前提知識
二重積分の定義
リーマン和による定義
$xy$ 平面上の閉領域 $D$ 上で定義された関数 $f(x, y)$ の二重積分を定義します。
領域 $D$ を $n$ 個の小領域 $\Delta D_1, \Delta D_2, \dots, \Delta D_n$ に分割し、各小領域の面積を $\Delta A_i$、小領域内の代表点を $(x_i^*, y_i^*)$ とします。リーマン和を作ると：
$$ S_n = \sum_{i=1}^{n} f(x_i^*, y_i^*)\,\Delta A_i $$
分割を無限に細かくしたときの極限が二重積分です：
$$ \boxed{\iint_D f(x, y)\,dA = \lim_{n \to \infty}\sum_{i=1}^{n} f(x_i^*, y_i^*)\,\Delta A_i} $$
ここで $dA = dx\,dy$ は面積要素です。
幾何学的意味
$f(x, y) \geq 0$ のとき、$\iint_D f(x, y)\,dA$ は曲面 $z = f(x,y)$ と $xy$ 平面の間の体積を表します。
逐次積分（累次積分）
矩形領域の場合
領域 $D = [a, b] \times [c, d]$（長方形）の場合、二重積分は1変数の積分を2回繰り返す逐次積分として計算できます：
$$ \boxed{\iint_D f(x, y)\,dA = \int_a^b \left[\int_c^d f(x, y)\,dy\right]dx = \int_c^d \left[\int_a^b f(x, y)\,dx\right]dy} $$
内側の積分では、もう一方の変数を定数として扱います。
具体例
$f(x, y) = x^2 y$ を $D = [0, 2] \times [0, 3]$ 上で積分します。
$$ \iint_D x^2 y\,dA = \int_0^2 \left[\int_0^3 x^2 y\,dy\right]dx $$
まず内側の積分（$x$ を定数として $y$ で積分）：
$$ \int_0^3 x^2 y\,dy = x^2 \left[\frac{y^2}{2}\right]_0^3 = x^2 \cdot \frac{9}{2} = \frac{9x^2}{2} $$
次に外側の積分：
$$ \int_0^2 \frac{9x^2}{2}\,dx = \frac{9}{2}\left[\frac{x^3}{3}\right]_0^2 = \frac{9}{2} \cdot \frac{8}{3} = 12 $$
$$ \therefore \iint_D x^2 y\,dA = 12 $$
一般領域の場合
領域が矩形でない場合、積分の上下限が変数の関数になります。
縦線先型（Type I）: $D = \{(x, y) \mid a \leq x \leq b,\; g_1(x) \leq y \leq g_2(x)\}$
$$ \iint_D f(x, y)\,dA = \int_a^b \int_{g_1(x)}^{g_2(x)} f(x, y)\,dy\,dx $$
横線先型（Type II）: $D = \{(x, y) \mid c \leq y \leq d,\; h_1(y) \leq x \leq h_2(y)\}$
$$ \iint_D f(x, y)\,dA = \int_c^d \int_{h_1(y)}^{h_2(y)} f(x, y)\,dx\,dy $$
具体例：三角形領域
$f(x, y) = xy$ を三角形 $D = \{(x,y) \mid 0 \leq x \leq 1,\; 0 \leq y \leq x\}$ 上で積分します。
$$ \iint_D xy\,dA = \int_0^1 \int_0^x xy\,dy\,dx $$
内側の積分：
$$ \int_0^x xy\,dy = x\left[\frac{y^2}{2}\right]_0^x = \frac{x^3}{2} $$
外側の積分：
$$ \int_0^1 \frac{x^3}{2}\,dx = \frac{1}{2}\left[\frac{x^4}{4}\right]_0^1 = \frac{1}{8} $$
積分順序の交換
同じ二重積分でも、$x$ を先に積分するか $y$ を先に積分するかで計算の難易度が変わることがあります。積分順序を交換するには、積分領域を正しく書き直すことが必要です。
フビニの定理
$f(x, y)$ が領域 $D$ 上で連続（またはルベーグ可積分）ならば：
$$ \int_a^b \int_{g_1(x)}^{g_2(x)} f(x,y)\,dy\,dx = \int_c^d \int_{h_1(y)}^{h_2(y)} f(x,y)\,dx\,dy $$
ここで両辺の積分限は同じ領域 $D$ を異なる方法で記述したものです。
具体例
$$ I = \int_0^1 \int_x^1 e^{y^2}\,dy\,dx $$
内側の積分 $\int_x^1 e^{y^2}\,dy$ は初等関数で計算できません。積分順序を交換しましょう。
領域は $D = \{(x, y) \mid 0 \leq x \leq 1,\; x \leq y \leq 1\}$ です。これを $y$ を先に固定する形に書き直すと：
$$ D = \{(x, y) \mid 0 \leq y \leq 1,\; 0 \leq x \leq y\} $$
順序を交換すると：
$$ I = \int_0^1 \int_0^y e^{y^2}\,dx\,dy = \int_0^1 y e^{y^2}\,dy $$
$u = y^2$ とおくと $du = 2y\,dy$ より：
$$ I = \frac{1}{2}\int_0^1 e^u\,du = \frac{1}{2}(e – 1) $$
積分順序の交換により、計算不可能だった積分が容易に求まりました。
三重積分
定義
3次元領域 $V$ 上の関数 $f(x, y, z)$ の三重積分は：
$$ \boxed{\iiint_V f(x, y, z)\,dV = \iiint_V f(x, y, z)\,dx\,dy\,dz} $$
計算は逐次積分を3回行います。
具体例：直方体
$f(x, y, z) = xyz$ を $V = [0, 1] \times [0, 2] \times [0, 3]$ 上で積分します。
$$ \iiint_V xyz\,dV = \int_0^1 \int_0^2 \int_0^3 xyz\,dz\,dy\,dx $$
最内側の積分：
$$ \int_0^3 xyz\,dz = xy\left[\frac{z^2}{2}\right]_0^3 = \frac{9xy}{2} $$
中間の積分：
$$ \int_0^2 \frac{9xy}{2}\,dy = \frac{9x}{2}\left[\frac{y^2}{2}\right]_0^2 = 9x $$
最外側の積分：
$$ \int_0^1 9x\,dx = 9\left[\frac{x^2}{2}\right]_0^1 = \frac{9}{2} $$
体積・質量・重心の計算
体積
領域 $V$ の体積は $f = 1$ として：
$$ \text{Volume} = \iiint_V dV $$
質量
密度分布 $\rho(x, y, z)$ をもつ物体の全質量：
$$ \boxed{M = \iiint_V \rho(x, y, z)\,dV} $$
重心
$$ \boxed{\bar{x} = \frac{1}{M}\iiint_V x\,\rho\,dV, \quad \bar{y} = \frac{1}{M}\iiint_V y\,\rho\,dV, \quad \bar{z} = \frac{1}{M}\iiint_V z\,\rho\,dV} $$
具体例：半球の重心
半径 $R$ の上半球（一様密度 $\rho_0$）の重心の $z$ 座標を求めます。
$$ V = \{(x, y, z) \mid x^2 + y^2 + z^2 \leq R^2,\; z \geq 0\} $$
全質量は $M = \rho_0 \cdot \frac{2}{3}\pi R^3$ です。
$z$ 座標の1次モーメント（球座標で計算すると効率的ですが、ここでは直交座標で立式します）：
$$ \iiint_V z\,\rho_0\,dV = \rho_0 \int_{-R}^{R}\int_{-\sqrt{R^2-x^2}}^{\sqrt{R^2-x^2}}\int_0^{\sqrt{R^2-x^2-y^2}} z\,dz\,dy\,dx $$
最内側の積分：
$$ \int_0^{\sqrt{R^2-x^2-y^2}} z\,dz = \frac{R^2-x^2-y^2}{2} $$
残りの積分は極座標 $x = r\cos\theta$, $y = r\sin\theta$ で：
$$ \int_0^{2\pi}\int_0^R \frac{R^2 – r^2}{2}\,r\,dr\,d\theta = \pi\int_0^R (R^2 r – r^3)\,dr = \pi\left[\frac{R^2 r^2}{2} – \frac{r^4}{4}\right]_0^R = \frac{\pi R^4}{4} $$
したがって：
$$ \bar{z} = \frac{\rho_0 \cdot \pi R^4/4}{\rho_0 \cdot 2\pi R^3/3} = \frac{3R}{8} $$
$$ \boxed{\bar{z} = \frac{3R}{8}} $$
半球の重心は底面から $3R/8$ の高さにあります。
Pythonでの可視化
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import integrate
# ============================
# 重積分の可視化と数値計算
# ============================
fig = plt.figure(figsize=(18, 6))
# === (1) 二重積分のリーマン和による可視化 ===
ax1 = fig.add_subplot(131, projection='3d')
def f_2d(x, y):
return np.sin(x) * np.cos(y) + 2
x = np.linspace(0, np.pi, 100)
y = np.linspace(0, np.pi, 100)
X, Y = np.meshgrid(x, y)
Z = f_2d(X, Y)
# 曲面
ax1.plot_surface(X, Y, Z, alpha=0.3, cmap='viridis')
# リーマン和の直方体
n_boxes = 8
x_edges = np.linspace(0, np.pi, n_boxes + 1)
y_edges = np.linspace(0, np.pi, n_boxes + 1)
for i in range(n_boxes):
for j in range(n_boxes):
xc = (x_edges[i] + x_edges[i+1]) / 2
yc = (y_edges[j] + y_edges[j+1]) / 2
zc = f_2d(xc, yc)
dx = x_edges[i+1] - x_edges[i]
dy = y_edges[j+1] - y_edges[j]
# 上面
xx = [x_edges[i], x_edges[i+1], x_edges[i+1], x_edges[i]]
yy = [y_edges[j], y_edges[j], y_edges[j+1], y_edges[j+1]]
zz = [zc, zc, zc, zc]
ax1.plot_surface(np.array([[xx[0], xx[1]], [xx[3], xx[2]]]),
np.array([[yy[0], yy[1]], [yy[3], yy[2]]]),
np.array([[zc, zc], [zc, zc]]),
alpha=0.2, color='orange')
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
ax1.set_title('Double integral\nas Riemann sum', fontsize=12)
ax1.view_init(elev=20, azim=-60)
# === (2) 積分順序の交換の可視化 ===
ax2 = fig.add_subplot(132)
# D = {(x,y) | 0<=x<=1, x<=y<=1}
x_fill = np.linspace(0, 1, 100)
ax2.fill_between(x_fill, x_fill, 1, alpha=0.3, color='blue',
label='Integration region $D$')
ax2.plot(x_fill, x_fill, 'b-', lw=2)
ax2.axhline(1, color='blue', lw=2)
ax2.axvline(0, color='blue', lw=2, ymin=0, ymax=1)
# Type I: 外側x, 内側y
ax2.annotate('Type I:\n$\\int_0^1\\int_x^1 f\\,dy\\,dx$',
xy=(0.15, 0.7), fontsize=11, color='red',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
# Type II 方向の矢印
ax2.annotate('Type II:\n$\\int_0^1\\int_0^y f\\,dx\\,dy$',
xy=(0.5, 0.3), fontsize=11, color='green',
bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.5))
# 縦方向のスライス
for xi in [0.2, 0.5, 0.8]:
ax2.plot([xi, xi], [xi, 1], 'r-', lw=1.5, alpha=0.6)
ax2.plot(xi, xi, 'ro', ms=4)
ax2.plot(xi, 1, 'ro', ms=4)
# 横方向のスライス
for yi in [0.3, 0.6, 0.9]:
ax2.plot([0, yi], [yi, yi], 'g-', lw=1.5, alpha=0.6)
ax2.plot(0, yi, 'go', ms=4)
ax2.plot(yi, yi, 'go', ms=4)
ax2.set_xlabel('$x$', fontsize=12)
ax2.set_ylabel('$y$', fontsize=12)
ax2.set_title('Changing integration order', fontsize=12)
ax2.set_xlim(-0.1, 1.2)
ax2.set_ylim(-0.1, 1.2)
ax2.set_aspect('equal')
ax2.grid(True, alpha=0.3)
# === (3) 数値積分の収束 ===
ax3 = fig.add_subplot(133)
# 二重積分 ∫∫_D x^2*y dA, D=[0,2]x[0,3] の真の値は12
exact_val = 12.0
# リーマン和の収束
n_values = np.arange(2, 101)
riemann_sums = []
for n in n_values:
dx = 2.0 / n
dy = 3.0 / n
total = 0.0
for i in range(n):
for j in range(n):
xc = (i + 0.5) * dx
yc = (j + 0.5) * dy
total += xc**2 * yc * dx * dy
riemann_sums.append(total)
ax3.plot(n_values, riemann_sums, 'b-', lw=1.5, label='Riemann sum')
ax3.axhline(exact_val, color='r', ls='--', lw=2, label=f'Exact = {exact_val}')
ax3.set_xlabel('Number of subdivisions $n$', fontsize=12)
ax3.set_ylabel('Integral value', fontsize=12)
ax3.set_title('Convergence of Riemann sum\n'
'$\\iint x^2 y\\, dA$ on $[0,2]\\times[0,3]$', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('multiple_integrals.png', dpi=150, bbox_inches='tight')
plt.show()
# === (4) 体積と重心の計算 ===
fig2, axes2 = plt.subplots(1, 2, figsize=(14, 6))
# SciPyで数値積分
# 半球 x^2+y^2+z^2 <= R^2, z >= 0 の重心z座標
R = 1.0
# 体積の数値計算
def vol_integrand_z(z, y, x):
return 1.0
def y_lower(x):
return -np.sqrt(max(R**2 - x**2, 0))
def y_upper(x):
return np.sqrt(max(R**2 - x**2, 0))
def z_lower(y, x):
return 0
def z_upper(y, x):
val = R**2 - x**2 - y**2
return np.sqrt(max(val, 0))
vol_num, _ = integrate.tplquad(vol_integrand_z,
-R, R,
y_lower, y_upper,
z_lower, z_upper)
# z座標の1次モーメント
def moment_z(z, y, x):
return z
moment_z_num, _ = integrate.tplquad(moment_z,
-R, R,
y_lower, y_upper,
z_lower, z_upper)
z_bar_num = moment_z_num / vol_num
z_bar_exact = 3*R/8
# 結果の可視化
# 半球の断面図と重心位置
theta = np.linspace(0, np.pi, 100)
x_circle = R * np.cos(theta)
z_circle = R * np.sin(theta)
axes2[0].fill_between(x_circle, 0, z_circle, alpha=0.3, color='blue')
axes2[0].plot(x_circle, z_circle, 'b-', lw=2)
axes2[0].axhline(0, color='k', lw=1)
# 重心位置
axes2[0].plot(0, z_bar_exact, 'r*', ms=20, label=f'$\\bar{{z}}$ = {z_bar_exact:.4f} (exact)')
axes2[0].plot(0.05, z_bar_num, 'go', ms=10, label=f'$\\bar{{z}}$ = {z_bar_num:.4f} (numerical)')
axes2[0].axhline(z_bar_exact, color='r', ls='--', alpha=0.5)
axes2[0].set_xlabel('$x$', fontsize=12)
axes2[0].set_ylabel('$z$', fontsize=12)
axes2[0].set_title(f'Hemisphere centroid\n$R = {R}$, $\\bar{{z}} = 3R/8$', fontsize=12)
axes2[0].legend(fontsize=10)
axes2[0].set_aspect('equal')
axes2[0].grid(True, alpha=0.3)
# 2D: 密度が変化する板の重心
# ρ(x,y) = x + y, D: 単位正方形 [0,1]^2
x_2d = np.linspace(0, 1, 100)
y_2d = np.linspace(0, 1, 100)
X2, Y2 = np.meshgrid(x_2d, y_2d)
rho = X2 + Y2
M_2d, _ = integrate.dblquad(lambda y, x: x + y, 0, 1, 0, 1)
Mx_2d, _ = integrate.dblquad(lambda y, x: x*(x + y), 0, 1, 0, 1)
My_2d, _ = integrate.dblquad(lambda y, x: y*(x + y), 0, 1, 0, 1)
xbar_2d = Mx_2d / M_2d
ybar_2d = My_2d / M_2d
cs = axes2[1].contourf(X2, Y2, rho, levels=20, cmap='YlOrRd')
plt.colorbar(cs, ax=axes2[1], label='$\\rho(x,y) = x + y$')
axes2[1].plot(xbar_2d, ybar_2d, 'k*', ms=20,
label=f'Centroid ({xbar_2d:.3f}, {ybar_2d:.3f})')
axes2[1].set_xlabel('$x$', fontsize=12)
axes2[1].set_ylabel('$y$', fontsize=12)
axes2[1].set_title('Centroid of non-uniform plate\n$\\rho(x,y) = x + y$', fontsize=12)
axes2[1].legend(fontsize=10, loc='lower right')
axes2[1].set_aspect('equal')
axes2[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('centroid_calculation.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"半球の体積（数値）: {vol_num:.6f}")
print(f"半球の体積（解析）: {2/3*np.pi*R**3:.6f}")
print(f"半球の重心z（数値）: {z_bar_num:.6f}")
print(f"半球の重心z（解析）: {z_bar_exact:.6f}")
print(f"2D板の重心: ({xbar_2d:.4f}, {ybar_2d:.4f})")
上段左の図は、曲面の下にリーマン和の直方体を描いたもので、二重積分が「曲面と底面の間の体積」を細かな直方体の和で近似する操作であることを視覚的に示しています。中央の図は、三角形領域 $D$ での積分順序の交換を図示しています。赤い縦線がType I（$x$ を固定して $y$ で積分）、緑の横線がType II（$y$ を固定して $x$ で積分）に対応します。右の図は、分割数を増やすにつれてリーマン和が真の値 12 に収束する様子を示しています。
下段の左の図は、半球の断面と重心位置 $\bar{z} = 3R/8$ を示しています。数値積分の結果と解析解が一致していることを確認できます。右の図は、密度 $\rho(x,y) = x + y$ の非一様な板の密度分布と重心位置を示しており、密度が大きい方向に重心が偏っていることがわかります。
まとめ
本記事では、重積分の定義と計算方法について解説しました。
- 二重積分はリーマン和の極限として定義され、曲面の下の体積を表す
- 逐次積分により、二重積分を1変数の積分の繰り返しとして計算できる
- 積分順序の交換（フビニの定理）により、計算が大幅に簡単になる場合がある
- 三重積分は同様に3回の逐次積分として計算する
- 質量・重心は密度関数を重積分して求まる
- 半球の重心は底面から $3R/8$ の高さにある
次のステップとして、以下の記事も参考にしてください。