ベクトル解析の微分演算子には勾配（grad）、発散（div）、回転（curl）がありますが、これらを組み合わせたもう一つの重要な演算子が ラプラシアン（Laplacian）$\nabla^2$ です。ラプラシアンは「勾配の発散」として定義され、イメージとしては「ある点の値が周囲の平均からどれだけずれているか」を測る量です。
ラプラシアンは物理学の根幹をなす偏微分方程式に登場します。熱伝導方程式、波動方程式、シュレーディンガー方程式、静電場のポアソン方程式など、理工学の至る所で $\nabla^2$ が現れます。ラプラシアンを理解することは、これらの偏微分方程式を学ぶための必須の準備です。
本記事の内容
- ラプラシアンの定義（$\nabla^2 f = \nabla \cdot (\nabla f)$）
- ラプラシアンの物理的意味（周囲との差）
- 調和関数の定義と性質
- ラプラス方程式とポアソン方程式
- 各座標系でのラプラシアンの表現
- Pythonでのラプラシアンの計算と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ラプラシアンとは
直感的な理解
ラプラシアンの物理的意味を一言で言うと、「ある点の値が、その周囲の平均値からどれだけずれているか」 です。
温度場 $T(x, y, z)$ で考えてみましょう。ある点の温度が周囲の平均温度より高ければ、その点の温度は時間とともに下がるはずです（熱が周囲に流れ出す）。逆に周囲より低ければ温度は上がります。この「周囲との差」を表すのがラプラシアン $\nabla^2 T$ です。
- $\nabla^2 T > 0$：周囲の平均より値が 低い → 温度は上昇する方向
- $\nabla^2 T < 0$：周囲の平均より値が 高い → 温度は下降する方向
- $\nabla^2 T = 0$：周囲の平均と 一致 → 平衡状態
1次元での直感
1次元では $\nabla^2 f = f”(x)$（2階微分）です。$f”(x) > 0$ は $f$ が下に凸（谷の底）、$f”(x) < 0$ は上に凸（山の頂上）を意味します。多次元のラプラシアンはこの概念をすべての方向に拡張したものです。
ラプラシアンの数学的定義
スカラー場のラプラシアン
スカラー場 $f(x, y, z)$ の ラプラシアン（Laplacian）は、勾配の発散として定義されます。
$$ \begin{equation} \nabla^2 f = \Delta f = \nabla \cdot (\nabla f) = \mathrm{div}(\mathrm{grad}\, f) \end{equation} $$
具体的に書くと：
$$ \begin{equation} \nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2} \end{equation} $$
すなわち、各変数の 2階偏微分の和 です。
導出
勾配は：
$$ \nabla f = \left(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z}\right) $$
この勾配の発散を取ると：
$$ \begin{align} \nabla \cdot (\nabla f) &= \frac{\partial}{\partial x}\frac{\partial f}{\partial x} + \frac{\partial}{\partial y}\frac{\partial f}{\partial y} + \frac{\partial}{\partial z}\frac{\partial f}{\partial z} \\ &= \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2} \end{align} $$
ベクトル場のラプラシアン
ベクトル場 $\bm{F} = (F_x, F_y, F_z)$ のラプラシアンは、各成分にスカラーラプラシアンを適用したものです。
$$ \begin{equation} \nabla^2 \bm{F} = (\nabla^2 F_x, \nabla^2 F_y, \nabla^2 F_z) \end{equation} $$
ベクトル恒等式を用いると：
$$ \nabla^2 \bm{F} = \nabla(\nabla \cdot \bm{F}) – \nabla \times (\nabla \times \bm{F}) $$
周囲の平均との関係
ラプラシアンが「周囲との差」を表すことを数学的に確認します。点 $(x_0, y_0)$ を中心とする半径 $\epsilon$ の円上での $f$ の平均値を $\bar{f}_\epsilon$ とすると：
$$ \begin{align} \bar{f}_\epsilon &= \frac{1}{2\pi} \int_0^{2\pi} f(x_0 + \epsilon \cos\theta, y_0 + \epsilon \sin\theta) \, d\theta \end{align} $$
テイラー展開を用いると：
$$ \begin{align} \bar{f}_\epsilon &\approx f(x_0, y_0) + \frac{\epsilon^2}{4} \left(\frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}\right) + O(\epsilon^4) \\ &= f(x_0, y_0) + \frac{\epsilon^2}{4} \nabla^2 f + O(\epsilon^4) \end{align} $$
したがって：
$$ \nabla^2 f = \lim_{\epsilon \to 0} \frac{4}{\epsilon^2} (\bar{f}_\epsilon – f) $$
ラプラシアンは確かに「周囲の平均と中心値の差」に比例します。
調和関数
定義
スカラー場 $f$ が 調和関数（harmonic function）であるとは：
$$ \begin{equation} \nabla^2 f = 0 \end{equation} $$
を満たすことです。この方程式を ラプラス方程式（Laplace equation）と呼びます。
調和関数の性質
(1) 平均値の性質: 調和関数の任意の点での値は、その点を中心とする球面（2Dでは円）上での平均値に等しいです。
$$ f(\bm{r}_0) = \frac{1}{4\pi R^2} \oint_{|\bm{r} – \bm{r}_0| = R} f(\bm{r}) \, dS $$
(2) 最大値原理: 調和関数は（定数でない限り）領域の内部で最大値・最小値を取りません。極値は必ず境界上にあります。
(3) 一意性: ラプラス方程式の解は、境界条件を指定すれば一意に定まります。
調和関数の例
| 関数 | 次元 | 確認 |
|---|---|---|
| $f = ax + by + c$ | 2D | $\nabla^2 f = 0 + 0 = 0$ |
| $f = x^2 – y^2$ | 2D | $\nabla^2 f = 2 + (-2) = 0$ |
| $f = e^x \cos y$ | 2D | $\nabla^2 f = e^x \cos y + (-e^x \cos y) = 0$ |
| $f = \frac{1}{r}$ | 3D | $\nabla^2 f = 0$ ($r \neq 0$) |
| $f = \ln r$ | 2D | $\nabla^2 f = 0$ ($r \neq 0$) |
$f = x^2 – y^2$ の確認:
$$ \begin{align} \frac{\partial^2 f}{\partial x^2} &= \frac{\partial}{\partial x}(2x) = 2 \\ \frac{\partial^2 f}{\partial y^2} &= \frac{\partial}{\partial y}(-2y) = -2 \\ \nabla^2 f &= 2 + (-2) = 0 \quad \checkmark \end{align} $$
$f = e^x \cos y$ の確認:
$$ \begin{align} \frac{\partial^2 f}{\partial x^2} &= e^x \cos y \\ \frac{\partial^2 f}{\partial y^2} &= -e^x \cos y \\ \nabla^2 f &= e^x \cos y – e^x \cos y = 0 \quad \checkmark \end{align} $$
ラプラス方程式とポアソン方程式
ラプラス方程式
$$ \begin{equation} \nabla^2 f = 0 \end{equation} $$
源がない領域（電荷がない領域の電位、熱源がない定常温度場など）を記述します。
ポアソン方程式
$$ \begin{equation} \nabla^2 f = g(\bm{r}) \end{equation} $$
源 $g(\bm{r})$ がある場合の方程式です。例えば：
- 静電場: $\nabla^2 \phi = -\frac{\rho}{\varepsilon_0}$（$\rho$ は電荷密度）
- 重力場: $\nabla^2 U = 4\pi G \rho$（$\rho$ は質量密度）
- 定常熱伝導: $\nabla^2 T = -\frac{Q}{k}$（$Q$ は発熱率）
関連する偏微分方程式
ラプラシアンはさまざまな物理方程式に登場します。
| 方程式 | 式 | 物理 |
|---|---|---|
| ラプラス方程式 | $\nabla^2 f = 0$ | 定常状態（平衡） |
| ポアソン方程式 | $\nabla^2 f = g$ | 源を持つ定常状態 |
| 熱伝導方程式 | $\frac{\partial T}{\partial t} = \alpha \nabla^2 T$ | 温度の時間発展 |
| 波動方程式 | $\frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u$ | 波の伝播 |
| シュレーディンガー方程式 | $i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m}\nabla^2 \psi + V\psi$ | 量子力学 |
各座標系でのラプラシアン
直交座標系（デカルト座標）
$$ \nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2} $$
円柱座標系 $(r, \theta, z)$
$x = r\cos\theta$, $y = r\sin\theta$, $z = z$ のとき：
$$ \begin{equation} \nabla^2 f = \frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial f}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 f}{\partial \theta^2} + \frac{\partial^2 f}{\partial z^2} \end{equation} $$
展開すると：
$$ \nabla^2 f = \frac{\partial^2 f}{\partial r^2} + \frac{1}{r}\frac{\partial f}{\partial r} + \frac{1}{r^2}\frac{\partial^2 f}{\partial \theta^2} + \frac{\partial^2 f}{\partial z^2} $$
球座標系 $(r, \theta, \phi)$
$x = r\sin\theta\cos\phi$, $y = r\sin\theta\sin\phi$, $z = r\cos\theta$ のとき：
$$ \begin{equation} \nabla^2 f = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2 \frac{\partial f}{\partial r}\right) + \frac{1}{r^2 \sin\theta}\frac{\partial}{\partial \theta}\left(\sin\theta \frac{\partial f}{\partial \theta}\right) + \frac{1}{r^2 \sin^2\theta}\frac{\partial^2 f}{\partial \phi^2} \end{equation} $$
球対称関数 $f = f(r)$ の場合は大幅に簡略化されます：
$$ \nabla^2 f = \frac{1}{r^2}\frac{d}{dr}\left(r^2 \frac{df}{dr}\right) = \frac{d^2 f}{dr^2} + \frac{2}{r}\frac{df}{dr} $$
具体例
例1: $f(x, y) = x^2 + y^2$
$$ \nabla^2 f = \frac{\partial^2}{\partial x^2}(x^2 + y^2) + \frac{\partial^2}{\partial y^2}(x^2 + y^2) = 2 + 2 = 4 $$
$\nabla^2 f = 4 > 0$ なので、各点で値が周囲の平均より低い（放物面の谷の部分）。
例2: $f(r) = \frac{1}{r}$（$r = \sqrt{x^2 + y^2 + z^2}$, $r \neq 0$）
球座標で：
$$ \begin{align} \nabla^2 \frac{1}{r} &= \frac{1}{r^2}\frac{d}{dr}\left(r^2 \frac{d}{dr}\frac{1}{r}\right) \\ &= \frac{1}{r^2}\frac{d}{dr}\left(r^2 \cdot \left(-\frac{1}{r^2}\right)\right) \\ &= \frac{1}{r^2}\frac{d}{dr}(-1) \\ &= 0 \quad (r \neq 0) \end{align} $$
$f = 1/r$ は $r \neq 0$ で調和関数です。
例3: 熱伝導方程式
1次元の熱伝導方程式 $\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}$ は、$\nabla^2 T > 0$ の場所で温度が上昇し、$\nabla^2 T < 0$ の場所で温度が下降することを意味します。
Pythonでの実装
ラプラシアンの数値計算
import numpy as np
def laplacian_2d(f, dx, dy):
"""
2次元スカラー場のラプラシアンを中心差分で数値計算する
Parameters
----------
f : ndarray, shape (ny, nx)
スカラー場
dx, dy : float
格子間隔
Returns
-------
ndarray, shape (ny, nx)
ラプラシアン
"""
# 2階偏微分を中心差分で計算
d2f_dx2 = np.zeros_like(f)
d2f_dy2 = np.zeros_like(f)
# x方向の2階微分
d2f_dx2[:, 1:-1] = (f[:, 2:] - 2*f[:, 1:-1] + f[:, :-2]) / dx**2
# y方向の2階微分
d2f_dy2[1:-1, :] = (f[2:, :] - 2*f[1:-1, :] + f[:-2, :]) / dy**2
return d2f_dx2 + d2f_dy2
# テスト: f = x^2 + y^2 → ∇²f = 4
x = np.linspace(-3, 3, 200)
y = np.linspace(-3, 3, 200)
X, Y = np.meshgrid(x, y)
dx = x[1] - x[0]
dy = y[1] - y[0]
f = X**2 + Y**2
lap_f = laplacian_2d(f, dx, dy)
# 境界を除いた内部での評価
interior = lap_f[5:-5, 5:-5]
print(f"f = x^2 + y^2")
print(f"解析的な ∇²f = 4")
print(f"数値的な ∇²f (mean) = {interior.mean():.6f}")
print(f"数値的な ∇²f (std) = {interior.std():.6f}")
調和関数と非調和関数の可視化
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-3, 3, 200)
y = np.linspace(-3, 3, 200)
X, Y = np.meshgrid(x, y)
dx = x[1] - x[0]
dy = y[1] - y[0]
fig, axes = plt.subplots(2, 3, figsize=(15, 9))
# --- 行1: 関数のプロット ---
# f1 = x^2 - y^2（調和関数）
f1 = X**2 - Y**2
ax = axes[0, 0]
contour = ax.contourf(X, Y, f1, levels=20, cmap='RdBu_r')
ax.contour(X, Y, f1, levels=10, colors='black', linewidths=0.3)
plt.colorbar(contour, ax=ax)
ax.set_title('$f = x^2 - y^2$ (Harmonic)')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
# f2 = x^2 + y^2（非調和関数）
f2 = X**2 + Y**2
ax = axes[0, 1]
contour = ax.contourf(X, Y, f2, levels=20, cmap='viridis')
ax.contour(X, Y, f2, levels=10, colors='white', linewidths=0.3)
plt.colorbar(contour, ax=ax)
ax.set_title('$f = x^2 + y^2$ (Not harmonic)')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
# f3 = e^x cos(y)（調和関数）
f3 = np.exp(X) * np.cos(Y)
ax = axes[0, 2]
contour = ax.contourf(X, Y, f3, levels=20, cmap='RdBu_r')
ax.contour(X, Y, f3, levels=10, colors='black', linewidths=0.3)
plt.colorbar(contour, ax=ax)
ax.set_title('$f = e^x \\cos y$ (Harmonic)')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
# --- 行2: ラプラシアンのプロット ---
lap1 = laplacian_2d(f1, dx, dy)
ax = axes[1, 0]
contour = ax.contourf(X, Y, lap1, levels=20, cmap='RdBu_r')
plt.colorbar(contour, ax=ax)
ax.set_title(f'$\\nabla^2 f \\approx {lap1[50:-50, 50:-50].mean():.2f}$ (≈ 0)')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
lap2 = laplacian_2d(f2, dx, dy)
ax = axes[1, 1]
contour = ax.contourf(X, Y, lap2, levels=20, cmap='RdBu_r')
plt.colorbar(contour, ax=ax)
ax.set_title(f'$\\nabla^2 f \\approx {lap2[50:-50, 50:-50].mean():.2f}$ (≈ 4)')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
lap3 = laplacian_2d(f3, dx, dy)
ax = axes[1, 2]
contour = ax.contourf(X, Y, lap3, levels=20, cmap='RdBu_r')
plt.colorbar(contour, ax=ax)
ax.set_title(f'$\\nabla^2 f \\approx {lap3[50:-50, 50:-50].mean():.4f}$ (≈ 0)')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig("laplacian_harmonic.png", dpi=150, bbox_inches='tight')
plt.show()
平均値の性質の数値検証
import numpy as np
import matplotlib.pyplot as plt
def check_mean_value_property(f_func, x0, y0, R, N_theta=1000):
"""
調和関数の平均値の性質を数値検証する
Parameters
----------
f_func : callable
スカラー場の関数
x0, y0 : float
中心点
R : float
円の半径
N_theta : int
分割数
Returns
-------
tuple
(中心での値, 円上の平均値)
"""
# 中心での値
f_center = f_func(x0, y0)
# 円上での平均値
theta = np.linspace(0, 2*np.pi, N_theta, endpoint=False)
x_circle = x0 + R * np.cos(theta)
y_circle = y0 + R * np.sin(theta)
f_circle = np.array([f_func(xi, yi) for xi, yi in zip(x_circle, y_circle)])
f_mean = f_circle.mean()
return f_center, f_mean
# 調和関数 f = x^2 - y^2
print("=== 調和関数: f = x^2 - y^2 ===")
f_harmonic = lambda x, y: x**2 - y**2
for R in [0.5, 1.0, 2.0]:
center, mean = check_mean_value_property(f_harmonic, 1.0, 1.0, R)
print(f"R = {R:.1f}: f(center) = {center:.4f}, mean = {mean:.4f}, diff = {abs(center-mean):.6f}")
# 非調和関数 f = x^2 + y^2
print("\n=== 非調和関数: f = x^2 + y^2 ===")
f_not_harmonic = lambda x, y: x**2 + y**2
for R in [0.5, 1.0, 2.0]:
center, mean = check_mean_value_property(f_not_harmonic, 1.0, 1.0, R)
print(f"R = {R:.1f}: f(center) = {center:.4f}, mean = {mean:.4f}, diff = {abs(center-mean):.6f}")
ラプラス方程式の数値解法（ヤコビ法）
import numpy as np
import matplotlib.pyplot as plt
# 2Dラプラス方程式を正方形領域でヤコビ法により数値的に解く
# ∇²u = 0、境界条件: 上辺 u=100, 他の3辺 u=0
N = 50 # 格子点数
u = np.zeros((N, N))
# 境界条件
u[0, :] = 100.0 # 上辺（y=ymax）
u[-1, :] = 0.0 # 下辺
u[:, 0] = 0.0 # 左辺
u[:, -1] = 0.0 # 右辺
# ヤコビ反復法
max_iter = 5000
tolerance = 1e-5
for iteration in range(max_iter):
u_old = u.copy()
# 内部点の更新: u(i,j) = [u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1)] / 4
u[1:-1, 1:-1] = 0.25 * (u_old[:-2, 1:-1] + u_old[2:, 1:-1] +
u_old[1:-1, :-2] + u_old[1:-1, 2:])
# 収束判定
diff = np.max(np.abs(u - u_old))
if diff < tolerance:
print(f"収束: {iteration + 1} 回で tolerance {tolerance} 以下")
break
# 可視化
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
x = np.linspace(0, 1, N)
y = np.linspace(0, 1, N)
X, Y = np.meshgrid(x, y)
contour = ax.contourf(X, Y, u, levels=30, cmap='hot')
ax.contour(X, Y, u, levels=15, colors='white', linewidths=0.3)
plt.colorbar(contour, ax=ax, label='$u(x, y)$')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_title('Solution of $\\nabla^2 u = 0$\n(Top: $u=100$, Other sides: $u=0$)')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig("laplace_equation_solution.png", dpi=150, bbox_inches='tight')
plt.show()
# ラプラシアンが0に近いことを確認
dx = x[1] - x[0]
dy = y[1] - y[0]
lap_u = np.zeros_like(u)
lap_u[1:-1, 1:-1] = ((u[:-2, 1:-1] - 2*u[1:-1, 1:-1] + u[2:, 1:-1]) / dx**2 +
(u[1:-1, :-2] - 2*u[1:-1, 1:-1] + u[1:-1, 2:]) / dy**2)
print(f"内部の ∇²u の最大値: {np.max(np.abs(lap_u[2:-2, 2:-2])):.6f}")
各座標系でのラプラシアンの比較
import numpy as np
import matplotlib.pyplot as plt
# f = 1/r のラプラシアンを直交座標と球座標で計算・比較
# 直交座標で数値的に計算
x = np.linspace(-3, 3, 300)
y = np.linspace(-3, 3, 300)
X, Y = np.meshgrid(x, y)
dx = x[1] - x[0]
dy = y[1] - y[0]
# 2Dの f = 1/r（対数ポテンシャルの方が2Dでは調和関数）
R = np.sqrt(X**2 + Y**2)
R_safe = np.where(R < 0.3, 0.3, R) # 特異点を避ける
f = np.log(R_safe) # 2Dでは ln(r) が調和関数
# ラプラシアンの計算
lap_f = laplacian_2d(f, dx, dy)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# 関数 f = ln(r)
ax = axes[0]
contour = ax.contourf(X, Y, f, levels=20, cmap='viridis')
ax.contour(X, Y, f, levels=10, colors='white', linewidths=0.3)
plt.colorbar(contour, ax=ax, label='$f$')
ax.set_title('$f = \\ln r$ (2D harmonic for $r > 0$)')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
# ラプラシアン
ax = axes[1]
# 特異点付近を除外して表示
lap_clipped = np.clip(lap_f, -5, 5)
contour = ax.contourf(X, Y, lap_clipped, levels=20, cmap='RdBu_r')
plt.colorbar(contour, ax=ax, label='$\\nabla^2 f$')
ax.set_title('$\\nabla^2 (\\ln r) \\approx 0$ (away from origin)')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig("laplacian_log_r.png", dpi=150, bbox_inches='tight')
plt.show()
# r > 1 の領域でのラプラシアンの値を確認
mask = R > 1.0
print(f"r > 1 の領域での ∇²(ln r) の平均: {lap_f[mask].mean():.6f}")
print(f"r > 1 の領域での ∇²(ln r) の標準偏差: {lap_f[mask].std():.6f}")
まとめ
本記事では、ラプラシアンと調和関数について解説しました。
- ラプラシアン $\nabla^2 f = \mathrm{div}(\mathrm{grad}\, f)$ は各変数の 2階偏微分の和 である
- 物理的意味: ある点の値と周囲の平均値との差 に比例する量である
- $\nabla^2 f = 0$ を満たす関数を 調和関数 と呼び、平均値の性質と最大値原理を持つ
- ラプラス方程式 $\nabla^2 f = 0$ は平衡状態、ポアソン方程式 $\nabla^2 f = g$ は源を持つ場を記述する
- 熱伝導方程式、波動方程式、シュレーディンガー方程式などの核心にラプラシアンが存在する
- 円柱座標・球座標では追加の項が現れ、特に球対称の場合は大幅に簡略化される
次のステップとして、以下の記事も参考にしてください。