ベクトル解析の微分演算子には勾配（grad）、発散（div）、回転（curl）がありますが、これらを組み合わせたもう一つの重要な演算子が ラプラシアン（Laplacian）$\nabla^2$ です。ラプラシアンは「勾配の発散」として定義され、イメージとしては「ある点の値が周囲の平均からどれだけずれているか」を測る量です。
ラプラシアンは物理学の根幹をなす偏微分方程式に登場します。熱伝導方程式、波動方程式、シュレーディンガー方程式、静電場のポアソン方程式など、理工学の至る所で $\nabla^2$ が現れます。ラプラシアンを理解することは、これらの偏微分方程式を学ぶための必須の準備です。
本記事の内容
- ラプラシアンの定義（$\nabla^2 f = \nabla \cdot (\nabla f)$）
- ラプラシアンの物理的意味（周囲との差）
- 調和関数の定義と性質
- ラプラス方程式とポアソン方程式
- 各座標系でのラプラシアンの表現
- Pythonでのラプラシアンの計算と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ラプラシアンとは
ラプラシアンの核心は図の一枚に尽きます。ある点の値 $f_0$ と、その周囲の点の平均との差——これがラプラシアンです。中心が周囲の平均より低ければ $\nabla^2 f > 0$、高ければ $\nabla^2 f < 0$ になります。
直感的な理解
ラプラシアンの物理的意味を一言で言うと、「ある点の値が、その周囲の平均値からどれだけずれているか」 です。
温度場 $T(x, y, z)$ で考えてみましょう。ある点の温度が周囲の平均温度より高ければ、その点の温度は時間とともに下がるはずです（熱が周囲に流れ出す）。逆に周囲より低ければ温度は上がります。この「周囲との差」を表すのがラプラシアン $\nabla^2 T$ です。
- $\nabla^2 T > 0$：周囲の平均より値が 低い → 温度は上昇する方向
- $\nabla^2 T < 0$：周囲の平均より値が 高い → 温度は下降する方向
- $\nabla^2 T = 0$：周囲の平均と 一致 → 平衡状態
1次元での直感
1次元で考えると、ラプラシアンは2階微分 $f”$ そのものです。図のように、下に凸（谷、$f”>0$）なら周囲より低く、上に凸（山、$f”<0$）なら周囲より高い。$f''$ の符号が「周囲との高低」を表しているのが見て取れます。
1次元では $\nabla^2 f = f”(x)$（2階微分）です。$f”(x) > 0$ は $f$ が下に凸（谷の底）、$f”(x) < 0$ は上に凸（山の頂上）を意味します。多次元のラプラシアンはこの概念をすべての方向に拡張したものです。
ラプラシアンの数学的定義
定義に入る前に、符号の意味を整理しておきます。図のように、$\nabla^2 f > 0$ は「周囲より低い（谷）」、$\nabla^2 f = 0$ は「周囲と同じ」、$\nabla^2 f < 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}) $$
周囲の平均との関係
2次元のスカラー場で確認しましょう。左がガウス型の山 $f$、右がそのラプラシアン $\nabla^2 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) $$
ラプラシアンは確かに「周囲の平均と中心値の差」に比例します。
調和関数
ラプラシアンがどこでもゼロ（$\nabla^2 f = 0$）になる関数を調和関数と呼びます。図は調和関数 $f = x^2 – y^2$ の等高線で、$\nabla^2 f = 2 – 2 = 0$ です。「どの点も周囲の平均に等しい」ため、極大も極小もなく、鞍点（峠）だけが現れるのが調和関数の特徴です。
定義
スカラー場 $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} $$
ラプラス方程式とポアソン方程式
ラプラシアンを使った2つの基本方程式を整理します。図のように、源（湧き出し）がない場は $\nabla^2 f = 0$（ラプラス方程式、例：真空中の静電ポテンシャル）、源 $\rho$ がある場は $\nabla^2 f = \rho$（ポアソン方程式、例：電荷分布のある静電場）で記述されます。
ラプラス方程式
$$ \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$ の役割がどの方程式でも同じであることが見えてきます。ラプラシアンは「その点の値が周囲の平均からどれだけずれているか」を測る量でした。それが何と等しいと置くかで、方程式の性格が決まります。
- ゼロと置く（ラプラス方程式）→ どの点も周囲の平均に一致している。もう動かない平衡状態
- 源と置く（ポアソン方程式）→ 電荷や質量のあるところだけ平均からずれる
- 時間の1階微分と置く（熱伝導方程式）→ 周囲より低ければ温まり、高ければ冷める。ずれが解消される向きに変化する
- 時間の2階微分と置く（波動方程式）→ ずれが「加速度」になる。行き過ぎて振動し、波として伝わる
- エネルギーと置く（シュレディンガー方程式）→ $-\frac{\hbar^2}{2m}\nabla^2$ が運動エネルギー演算子。波動関数の曲がり具合が運動エネルギーに対応する
熱伝導と波動の違いが、時間微分の階数だけである点は味わい深いところです。1階なら「ずれが減る方向にじわじわ進む」拡散、2階なら「ずれが加速度になって行き過ぎる」振動になります。同じ $\nabla^2$ が、拡散も波動も生み出しているわけです。
各座標系でのラプラシアン
ラプラシアンは座標系によって形が変わります。図のように、問題の対称性に合わせて選ぶのがコツです。平面・箱型なら直交座標、円筒対称なら円柱座標、球対称なら球座標を使うと、方程式が解きやすくなります。
直交座標系（デカルト座標）
$$ \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$
図のように $f = x^2 + y^2$ はお椀型で、ラプラシアンは $\nabla^2 f = 2 + 2 = 4$（どこでも定数）です。お椀のどの点も「周囲の平均より低い谷」なので、ラプラシアンが正の一定値になる、と直感的に理解できます。
$$ \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$）
図の $f = 1/r$ は、3次元では原点以外でラプラシアンがゼロ（調和関数）になります。これは点電荷がつくるクーロンポテンシャルの形そのもので、電荷のない空間で静電ポテンシャルがラプラス方程式を満たすことに対応します。
球座標で：
$$ \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: 熱伝導方程式
ラプラシアンの意味が最もよく現れるのが熱伝導方程式 $\partial u/\partial t = \nabla^2 u$ です。図は、最初に箱型だった温度分布が時間とともに平滑化されていく様子です。「周囲より高い点は下がり、低い点は上がる」というラプラシアンの作用が、凸凹をならして拡散させているのが見て取れます。
1次元の熱伝導方程式 $\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}$ は、$\nabla^2 T > 0$ の場所で温度が上昇し、$\nabla^2 T < 0$ の場所で温度が下降することを意味します。
離散ラプラシアン — 差分と画像処理
ここまでは連続の関数を扱ってきましたが、実際にコンピュータで計算するときは値が格子上にしかありません。そこで微分を差分に置き換えた離散ラプラシアンを使います。
1次元では、2階微分の中心差分から
$$ \frac{d^2 f}{dx^2} \approx \frac{f(x+h) – 2f(x) + f(x-h)}{h^2} $$
が得られます。2次元では $x$ 方向と $y$ 方向の両方を足して、格子幅を $h = 1$ とすれば
$$ \nabla^2 f_{i,j} \approx f_{i+1,j} + f_{i-1,j} + f_{i,j+1} + f_{i,j-1} – 4 f_{i,j} $$
という式になります。「上下左右の4つを足して、中心の4倍を引く」という形です。これはまさに「周囲の平均との差」の離散版で、$\frac{1}{4}(f_{i+1,j}+f_{i-1,j}+f_{i,j+1}+f_{i,j-1}) – f_{i,j}$ の4倍にあたります。
この係数を並べたものを5点ステンシルと呼びます。
2枚目がそのステンシルで、中心が $-4$、上下左右が $+1$、斜めは $0$ です。これを画像の各画素に適用した結果が3枚目、その絶対値が4枚目です。
読み取れることが2つあります。ひとつは、四角形や円の内部（値が一定な平坦部）では結果がぴったり $0$ になっていることです。実際に四角形の内部だけを取り出して測ると、最大値は厳密に $0$ でした。周囲と同じ値なら「平均とのずれ」はゼロ、というラプラシアンの意味そのものです。
もうひとつは、輪郭のところだけが大きな値を持つことです。3枚目では境界の内側が青（負）、外側が赤（正）になっており、値が急に変わる場所でラプラシアンが符号を変えながら大きくなっています。この符号が変わる点（ゼロ交差）が、エッジの正確な位置にあたります。
これが画像処理のラプラシアンフィルタの原理です。ノイズがあると2階微分は敏感に反応してしまうため、実務ではガウシアンフィルタで平滑化してからラプラシアンを掛ける LoG（Laplacian of Gaussian） がよく使われます。この記事の図でも、先に $\sigma = 1.4$ のガウシアンで平滑化してからラプラシアンを適用しています。
グラフラプラシアン — 格子からグラフへ
離散ラプラシアンをさらに一般化すると、格子である必要すらなくなります。「隣は誰か」という接続関係さえ決まっていれば、同じ発想が使えるからです。これがグラフラプラシアンです。
グラフの隣接行列を $A$、各ノードの次数（つながっている辺の本数）を並べた対角行列を $D$ とすると、グラフラプラシアンは
$$ L = D – A $$
で定義されます。
左は7個のノードを1列に並べ、隣どうしだけを結んだグラフです。中央がその $L = D – A$ で、内部のノードでは対角が $+2$、両隣が $-1$ になっています。これは1次元の離散ラプラシアン $f_{i+1} – 2f_i + f_{i-1}$ の符号を反転したものです。
符号の違いは慣習の問題です。連続のラプラシアンは平坦な点から凹むと正になりますが、グラフラプラシアンは $L$ の固有値がすべて $0$ 以上になるよう（半正定値になるよう）符号を取っています。両者を行き来するときに混乱しやすい点なので、意識しておくと安全です。
右は $L$ の固有ベクトルをプロットしたものです。第2固有ベクトルは端から端へ単調に変化し、第3は1回、第4は2回符号を変えます。これは弦の振動の基本モード・2倍音・3倍音とまったく同じ構造で、離散版の正弦波になっています。連続のラプラシアンの固有関数が三角関数だったことの、そのままの対応です。
固有値も特徴的で、最小固有値は必ず $0$ になります（この例でも厳密に $0$）。対応する固有ベクトルは「全ノードが同じ値」という定数ベクトルで、連続の世界で「定数関数のラプラシアンは $0$」だったことに対応します。2番目に小さい固有値はグラフの「つながりにくさ」を測る量として知られ、スペクトラルクラスタリングでグラフを2つに分ける基準に使われます。
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$ は源を持つ場を記述する
- 熱伝導方程式、波動方程式、シュレーディンガー方程式などの核心にラプラシアンが存在する
- 円柱座標・球座標では追加の項が現れ、特に球対称の場合は大幅に簡略化される
- 微分を差分に置き換えた 離散ラプラシアン（5点ステンシル）は、平坦部で $0$、輪郭で大きくなる。画像処理のエッジ検出はこの性質を使っている
- さらに一般化した グラフラプラシアン $L = D – A$ は固有値が非負で、第2固有値がグラフの分けやすさを表す。スペクトラルクラスタリングの土台になっている
次のステップとして、以下の記事も参考にしてください。