ヘルムホルツ分解 — 任意のベクトル場は回転なし成分と発散なし成分に分かれる

天気図の風の矢印を眺めていると、そこには二種類の「模様」が混ざっていることに気づきます。低気圧の中心に向かって渦を巻きながら回り込む流れと、高気圧から四方八方へ吹き出していく流れです。前者は回転が主役で、後者は湧き出しが主役です。では、任意の風の場を「純粋に渦だけの成分」と「純粋に湧き出しだけの成分」にきれいに仕分けることはできるのでしょうか。

できます。しかも、その仕分けは必ず存在し、しかも一通りしかない。これがヘルムホルツ分解(Helmholtz decomposition、ヘルムホルツの定理)です。式で書けば、適当な減衰条件を満たすベクトル場 $\bm{F}$ は必ず

$$ \bm{F} = -\nabla \phi + \nabla \times \bm{A} $$

の形に書けます。第1項は回転がゼロ、第2項は発散がゼロ。つまり「渦を持たない部分」と「湧き出しを持たない部分」への完全な仕分けです。

らせん状に広がる流れの場が、中心から放射状に外向きの湧き出し成分と、中心のまわりを回る渦成分の和に分解される様子

左のパネルは渦を巻きながら外へ広がる「らせん状」の流れで、渦と湧き出しのどちらが主役なのか目では判別できません。しかしこれは、中央の「放射状に外向きな成分」と右の「反時計回りに回るだけの成分」をぴったり足し合わせたものです。ヘルムホルツ分解が主張しているのは、こうした仕分けがどんな場に対しても必ずでき、しかも仕分け方は一通りしかない、ということです。

この定理は、抽象的なベクトル解析の一命題にとどまりません。応用先を二つだけ挙げます。ひとつは電磁気学です。マクスウェル方程式を解くときに登場するスカラーポテンシャル $\phi$ とベクトルポテンシャル $\bm{A}$、そして「クーロンゲージ $\nabla \cdot \bm{A} = 0$」という一見天下り的な条件は、実はヘルムホルツ分解そのものです。もうひとつは非圧縮流体のシミュレーションです。CG映画の煙も、航空機まわりの流れ解析も、その中核には「速度場から発散成分を引き算する」という操作があります。これは圧力射影法と呼ばれ、まさにヘルムホルツ分解の横波成分を取り出す操作に他なりません。さらに、画像処理での勾配場の積分、地磁気データの内部起源・外部起源の分離、弾性波動の P 波と S 波の分離も、すべて同じ数学の上に立っています。

本記事の内容

  • 「回転なし成分」「発散なし成分」とは何かを直感的につかむ
  • ポアソン方程式のグリーン関数から、分解を構成的に導出する(存在証明)
  • 無限遠での減衰条件のもとで分解が一意であることを証明する
  • フーリエ変換によって、分解が波数空間での射影行列 $\bm{P}_L = \bm{k}\bm{k}^{\mathsf T}/|\bm{k}|^2$、$\bm{P}_T = \bm{I} – \bm{P}_L$ に対応することを導く
  • 電磁場のポテンシャルと、流体の圧力射影法への応用
  • Python(FFT)で 2 次元場を実際に分解し、直交性と解析解との一致を数値検証する

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

特に $\nabla \times (\nabla \phi) = \bm{0}$ と $\nabla \cdot (\nabla \times \bm{A}) = 0$ という二つの恒等式、そして発散定理は本記事の全編で使います。

ヘルムホルツ分解とは — 場を「湧き出し」と「渦」に仕分ける

ベクトル場を調べたいとき、私たちは二つの微分演算を場に当てます。発散 $\nabla \cdot \bm{F}$ と回転 $\nabla \times \bm{F}$ です。発散は「その点で場がどれだけ湧き出しているか」、回転は「その点で場がどれだけ渦を巻いているか」を測ります。水面に浮かべた小さなゴミを想像すると分かりやすいでしょう。ゴミが周囲から離れていくように広がるなら発散が正、ゴミがくるくる自転するなら回転がゼロでない、というイメージです。

ここで自然な問いが生まれます。発散と回転の二つだけで、場の情報は尽くされているのでしょうか。 言い換えると、「湧き出しの分布」と「渦の分布」を全部教えてもらったら、元の場を復元できるのでしょうか。

ヘルムホルツの定理はこの問いに「イエス」と答えます。しかもその答え方が美しい。場そのものを二つのパーツに分け、片方は発散だけを担当し、もう片方は回転だけを担当する、という分業体制を作れると主張するのです。

言葉で書き下すと次のようになります。

  • 回転なし成分(縦成分、longitudinal part) $\bm{F}_L$:$\nabla \times \bm{F}_L = \bm{0}$ を満たす。あるスカラー関数 $\phi$ を使って $\bm{F}_L = -\nabla \phi$ と書ける。この成分だけが場の発散を担う
  • 発散なし成分(横成分、transverse part) $\bm{F}_T$:$\nabla \cdot \bm{F}_T = 0$ を満たす。あるベクトル場 $\bm{A}$ を使って $\bm{F}_T = \nabla \times \bm{A}$ と書ける。この成分だけが場の回転を担う

そして $\bm{F} = \bm{F}_L + \bm{F}_T$。これがヘルムホルツ分解です。$\phi$ をスカラーポテンシャル、$\bm{A}$ をベクトルポテンシャルと呼びます。

なぜ分業が成立するかは、恒等式から一瞬で見えます。$\bm{F} = -\nabla \phi + \nabla \times \bm{A}$ の発散を取ると、$\nabla \cdot (\nabla \times \bm{A}) = 0$ なので

$$ \nabla \cdot \bm{F} = -\nabla^2 \phi $$

となり、$\bm{A}$ は完全に消えます。同じように回転を取ると、$\nabla \times (\nabla \phi) = \bm{0}$ なので

$$ \nabla \times \bm{F} = \nabla \times (\nabla \times \bm{A}) $$

となり、今度は $\phi$ が完全に消えます。発散を測ると $\phi$ しか見えず、回転を測ると $\bm{A}$ しか見えない。 この「見えなくなる」性質こそが、分業を可能にしている仕掛けです。

スカラーポテンシャル $\phi$ に付いているマイナス符号は数学的には本質ではありませんが、物理では標準です。静電場が $\bm{E} = -\nabla \phi$ と書かれるのと同じ流儀で、「ポテンシャルの高いところから低いところへ場が向く」という描像に合わせています。

ここまでは「そう分解できたら嬉しいね」という願望の話です。次に、その願望が本当に叶うことを示さなければなりません。しかも、都合よく $\phi$ と $\bm{A}$ が存在することを祈るのではなく、与えられた $\bm{F}$ から実際に $\phi$ と $\bm{A}$ を作るレシピを書き下したい。そのための道具が、次のセクションで用意する二つの恒等式とポアソン方程式です。

準備:勾配場と回転場が持つ「消える」性質

分解の構成に入る前に、上で使った二つの恒等式を自分の手で確かめておきましょう。証明は成分計算だけです。

まず $\nabla \times (\nabla \phi) = \bm{0}$ から。回転の $x$ 成分を書き下すと

$$ [\nabla \times (\nabla \phi)]_x = \frac{\partial}{\partial y}\left(\frac{\partial \phi}{\partial z}\right) – \frac{\partial}{\partial z}\left(\frac{\partial \phi}{\partial y}\right) $$

です。$\phi$ が 2 回連続微分可能なら偏微分の順序を交換できます(シュワルツの定理)。すると二つの項は同じものになり、差し引きゼロです。$y$ 成分・$z$ 成分も添字を回すだけで同じ計算になるので、全成分がゼロです。

次に $\nabla \cdot (\nabla \times \bm{A}) = 0$。定義どおり展開すると

$$ \nabla \cdot (\nabla \times \bm{A}) = \frac{\partial}{\partial x}\left(\frac{\partial A_z}{\partial y} – \frac{\partial A_y}{\partial z}\right) + \frac{\partial}{\partial y}\left(\frac{\partial A_x}{\partial z} – \frac{\partial A_z}{\partial x}\right) + \frac{\partial}{\partial z}\left(\frac{\partial A_y}{\partial x} – \frac{\partial A_x}{\partial y}\right) $$

となります。ここで各 2 階偏導関数がどう現れるかを追いかけると、たとえば $\partial^2 A_z/\partial x \partial y$ は第 1 項に $+$、第 2 項に $-$ で現れます。6 個ある項がすべてこのようにペアを組んで打ち消し合うので、全体はゼロです。

勾配場では閉曲線に沿った循環がゼロ、回転場では閉曲線を横切る流束がゼロになることを数値積分で確かめた図

この二つの恒等式は、積分の言葉に直すと目で見えるようになります。左の勾配場では、破線の閉曲線に沿ってぐるりと一周したときの循環(接線方向成分の線積分)が $-1.6\times10^{-17}$、つまり丸め誤差レベルのゼロです。ストークスの定理により循環は内部の回転の総和ですから、これは「回転がゼロ」の別表現になっています。右の回転場では、同じ閉曲線を横切る流束(法線方向成分の線積分)が $1.6\times10^{-17}$ でやはりゼロ。発散定理により流束は内部の湧き出しの総和なので、「発散がゼロ」を意味します。閉曲線をどこにどう置いても結果は変わりません。

もう一つ、後で主役になる恒等式を用意します。ベクトル場の二重回転の公式です。

$$ \nabla \times (\nabla \times \bm{W}) = \nabla (\nabla \cdot \bm{W}) – \nabla^2 \bm{W} $$

ここで $\nabla^2 \bm{W}$ は各デカルト成分にラプラシアンを当てたものです(デカルト座標でのみこの単純な形になります)。$x$ 成分で確認しておくと、左辺は

$$ [\nabla \times (\nabla \times \bm{W})]_x = \frac{\partial}{\partial y}\left(\frac{\partial W_y}{\partial x} – \frac{\partial W_x}{\partial y}\right) – \frac{\partial}{\partial z}\left(\frac{\partial W_x}{\partial z} – \frac{\partial W_z}{\partial x}\right) $$

です。ここに $\partial^2 W_x/\partial x^2$ を足して引くという細工をします。足す側を第 1 項に回すと $\partial(\partial W_x/\partial x + \partial W_y/\partial y + \partial W_z/\partial z)/\partial x = \partial(\nabla\cdot\bm{W})/\partial x$ が組み上がり、引く側は残りの 2 階微分と合わさって $-\nabla^2 W_x$ になります。整理すると

$$ [\nabla \times (\nabla \times \bm{W})]_x = \frac{\partial}{\partial x}(\nabla \cdot \bm{W}) – \nabla^2 W_x $$

となり、確かに右辺の $x$ 成分と一致します。

この三つの恒等式を並べてみると、分解の設計図が見えてきます。二重回転の公式を移項すると

$$ -\nabla^2 \bm{W} = -\nabla (\nabla \cdot \bm{W}) + \nabla \times (\nabla \times \bm{W}) $$

です。右辺の第 1 項は「あるスカラーの勾配のマイナス」、第 2 項は「あるベクトルの回転」——まさにヘルムホルツ分解の形をしています。ということは、$-\nabla^2 \bm{W} = \bm{F}$ となる補助場 $\bm{W}$ さえ見つかれば、分解は自動的に完成するわけです。次のセクションでその $\bm{W}$ を実際に作ります。

分解を作る:ポアソン方程式のグリーン関数

$-\nabla^2 \bm{W} = \bm{F}$ を解くというのは、各成分についてポアソン方程式 $\nabla^2 W_x = -F_x$ を解くということです。電磁気学で「電荷分布からポテンシャルを求める」のとまったく同じ問題ですから、答えはよく知られたクーロン積分の形になります。それを丁寧に確認しましょう。

3 次元ラプラシアンの基本解

出発点は次の関係です。

$$ \nabla^2 \frac{1}{|\bm{r}|} = -4\pi \delta^3(\bm{r}) $$

「点電荷が作るポテンシャルは $1/r$ に比例する」という物理の事実を、数学の言葉で書いたものです。二段階で示します。

段階 1:$\bm{r} \neq \bm{0}$ ではゼロ。 球対称な関数に対するラプラシアンは、球座標で

$$ \nabla^2 f(r) = \frac{1}{r^2}\frac{d}{dr}\left(r^2 \frac{df}{dr}\right) $$

と書けます。$f = 1/r$ を入れると $df/dr = -1/r^2$、したがって $r^2 \, df/dr = -1$ という定数です。定数を $r$ で微分すればゼロなので、原点以外では $\nabla^2 (1/r) = 0$ が成り立ちます。

段階 2:原点まわりの積分値が $-4\pi$。 原点を含む半径 $\varepsilon$ の球 $B_\varepsilon$ で積分し、発散定理で表面積分に変えます。

$$ \int_{B_\varepsilon} \nabla^2 \frac{1}{r} \, dV = \oint_{\partial B_\varepsilon} \nabla \frac{1}{r} \cdot d\bm{S} $$

$\nabla (1/r) = -\hat{\bm{r}}/r^2$ で、球面上の面素は $d\bm{S} = \hat{\bm{r}}\, dS$ ですから、内積は $-1/\varepsilon^2$ という定数になります。球面の面積 $4\pi\varepsilon^2$ を掛けると

$$ \oint_{\partial B_\varepsilon} \nabla \frac{1}{r} \cdot d\bm{S} = -\frac{1}{\varepsilon^2} \cdot 4\pi \varepsilon^2 = -4\pi $$

半径 $\varepsilon$ に依存しない値になりました。「原点以外ではゼロ、原点まわりの積分は $-4\pi$」——これはまさにデルタ関数の $-4\pi$ 倍の振る舞いです。

1/r のグラフ、原点以外でラプラシアンがゼロになること、球面を貫く流束が半径によらず -4π になることを示した三面図

三つのパネルが、この二段階の証明をそのまま数値で追いかけています。中央のパネルは 3 次元の中心差分でラプラシアンを計算したもので、残差は最大でも $1.4\times10^{-5}$、$r$ が大きいところでは $10^{-8}$ 以下まで落ちます(残っているのは差分の打ち切り誤差で、真の値はゼロです)。右のパネルが決定的で、球の半径を $0.2$ から $8$ まで 40 倍変えても流束は $-4\pi = -12.566$ のまま動きません。半径を縮めても値が減らないということは、源が原点の一点に集中している——すなわちデルタ関数だということです。

補助場 $\bm{W}$ の構成

基本解が手に入ったので、$\bm{W}$ を積分で書けます。

$$ \begin{equation} \bm{W}(\bm{r}) = \frac{1}{4\pi}\int \frac{\bm{F}(\bm{r}’)}{|\bm{r} – \bm{r}’|} \, d^3 r’ \end{equation} $$

これが $-\nabla^2 \bm{W} = \bm{F}$ を満たすことを確認します。ラプラシアンは観測点 $\bm{r}$ についての微分なので、積分の中に入れられます。

$$ \nabla^2 \bm{W}(\bm{r}) = \frac{1}{4\pi}\int \bm{F}(\bm{r}’) \, \nabla_{\bm{r}}^2 \frac{1}{|\bm{r} – \bm{r}’|} \, d^3 r’ $$

ここで基本解の関係を $\bm{r} \to \bm{r} – \bm{r}’$ とずらして使うと $\nabla_{\bm{r}}^2 |\bm{r}-\bm{r}’|^{-1} = -4\pi \delta^3(\bm{r}-\bm{r}’)$ ですから、デルタ関数が積分を潰して

$$ \nabla^2 \bm{W}(\bm{r}) = \frac{1}{4\pi}\int \bm{F}(\bm{r}’) \cdot \left(-4\pi \delta^3(\bm{r}-\bm{r}’)\right) d^3 r’ = -\bm{F}(\bm{r}) $$

となります。狙いどおりです。

あとは前セクションの設計図に乗せるだけです。二重回転の公式に $\bm{W}$ を入れると

$$ \bm{F} = -\nabla^2 \bm{W} = -\nabla (\nabla \cdot \bm{W}) + \nabla \times (\nabla \times \bm{W}) $$

したがって

$$ \begin{equation} \phi = \nabla \cdot \bm{W}, \qquad \bm{A} = \nabla \times \bm{W} \end{equation} $$

と定義すれば $\bm{F} = -\nabla\phi + \nabla\times\bm{A}$ が成立します。分解の存在が、具体的な積分公式つきで示せました。 「存在するはず」ではなく「これが答えだ」と手渡せるのが、この構成的証明の強みです。

与えられた場からポアソン方程式で補助場を作り、二重回転の公式に代入してスカラーポテンシャルとベクトルポテンシャルを得るまでの流れ図

この流れ図が導出全体の見取り図です。上段の右向き矢印がポアソン方程式を解くステップで、ここだけがグリーン関数($1/r$)を必要とします。そこから下へ降りる矢印は二重回転の公式を当てるだけの代数操作、左向きの矢印はこのあと行う「源だけで書き直す」部分積分に対応します。ポイントは、どのステップにも恣意的な選択が入っていないことです。だからこそ次のセクションで一意性まで示せます。

ただしこの形の $\phi$ と $\bm{A}$ は、いったん $\bm{W}$ を経由しているぶん見通しが悪い。次のセクションで、これらを「$\bm{F}$ の発散」「$\bm{F}$ の回転」という源だけで書き直します。そうすれば、冒頭で述べた「湧き出しと渦の分布さえ分かれば場は復元できる」という主張が、文字どおりの形で現れます。

$\phi$ と $\bm{A}$ を源だけで書き直す

スカラーポテンシャル

$\phi = \nabla \cdot \bm{W}$ に $\bm{W}$ の積分表示を入れます。発散は $\bm{r}$ についての微分で、$\bm{F}(\bm{r}’)$ は $\bm{r}$ に依存しないので、微分は $1/|\bm{r}-\bm{r}’|$ にだけ効きます。

$$ \phi(\bm{r}) = \frac{1}{4\pi}\int \bm{F}(\bm{r}’) \cdot \nabla_{\bm{r}} \frac{1}{|\bm{r}-\bm{r}’|} \, d^3 r’ $$

ここで鍵になるのが、$|\bm{r}-\bm{r}’|$ が差にしか依存しないという性質です。$\bm{r}$ で微分するのと $\bm{r}’$ で微分するのは符号だけ違います。

$$ \nabla_{\bm{r}} \frac{1}{|\bm{r}-\bm{r}’|} = -\nabla_{\bm{r}’} \frac{1}{|\bm{r}-\bm{r}’|} $$

これを代入すると、積分変数 $\bm{r}’$ についての微分だけが残ります。

$$ \phi(\bm{r}) = -\frac{1}{4\pi}\int \bm{F}(\bm{r}’) \cdot \nabla’ \frac{1}{|\bm{r}-\bm{r}’|} \, d^3 r’ $$

ここで積の微分法則(スカラー $g$ とベクトル $\bm{F}$ に対する $\nabla’\cdot(g\bm{F}) = g\,\nabla’\cdot\bm{F} + \bm{F}\cdot\nabla’ g$)を $g = 1/|\bm{r}-\bm{r}’|$ に使って、被積分関数を書き換えます。

$$ \bm{F} \cdot \nabla’ \frac{1}{|\bm{r}-\bm{r}’|} = \nabla’ \cdot \left(\frac{\bm{F}(\bm{r}’)}{|\bm{r}-\bm{r}’|}\right) – \frac{\nabla’ \cdot \bm{F}(\bm{r}’)}{|\bm{r}-\bm{r}’|} $$

これを代入すると、$\phi$ は 2 項に分かれます。

$$ \phi(\bm{r}) = -\frac{1}{4\pi}\int \nabla’ \cdot \left(\frac{\bm{F}(\bm{r}’)}{|\bm{r}-\bm{r}’|}\right) d^3 r’ + \frac{1}{4\pi}\int \frac{\nabla’ \cdot \bm{F}(\bm{r}’)}{|\bm{r}-\bm{r}’|} \, d^3 r’ $$

第 1 項は全微分の形をしているので、発散定理で無限遠の球面上の面積分に変換できます。

$$ -\frac{1}{4\pi}\oint_{S_\infty} \frac{\bm{F}(\bm{r}’)}{|\bm{r}-\bm{r}’|} \cdot d\bm{S}’ $$

半径 $R$ の球面では面積が $4\pi R^2$、分母 $|\bm{r}-\bm{r}’| \sim R$ ですから、この面積分は $R \cdot |\bm{F}|$ 程度の大きさになります。したがって、$\bm{F}$ が無限遠で $1/R$ より速く減衰すれば この項はゼロに落ちます(実際には後述のとおり $1/R^{1+\epsilon}$ 程度あれば十分です)。この減衰条件を仮定すると、残るのは第 2 項だけです。

$$ \begin{equation} \phi(\bm{r}) = \frac{1}{4\pi}\int \frac{\nabla’ \cdot \bm{F}(\bm{r}’)}{|\bm{r}-\bm{r}’|} \, d^3 r’ \end{equation} $$

スカラーポテンシャルは、場の発散を「電荷」だと思ったときのクーロンポテンシャルそのものです。実際、基本解の性質からこの $\phi$ は

$$ \nabla^2 \phi = -\nabla \cdot \bm{F} $$

を満たします。冒頭で恒等式から予言したポアソン方程式が、ここで積分解つきで確認できました。

ベクトルポテンシャル

$\bm{A} = \nabla \times \bm{W}$ も同じ手順です。回転を積分の中に入れると、$\nabla_{\bm{r}} \times (g \bm{F}(\bm{r}’)) = (\nabla_{\bm{r}} g) \times \bm{F}(\bm{r}’)$ なので

$$ \bm{A}(\bm{r}) = \frac{1}{4\pi}\int \left(\nabla_{\bm{r}} \frac{1}{|\bm{r}-\bm{r}’|}\right) \times \bm{F}(\bm{r}’) \, d^3 r’ = -\frac{1}{4\pi}\int \left(\nabla’ \frac{1}{|\bm{r}-\bm{r}’|}\right) \times \bm{F}(\bm{r}’) \, d^3 r’ $$

ここでも積の微分法則、今度はベクトル版の $\nabla’ \times (g\bm{F}) = g \, \nabla’\times\bm{F} + (\nabla’ g)\times \bm{F}$ を使います。これを $(\nabla’ g) \times \bm{F}$ について解いて代入すると

$$ \bm{A}(\bm{r}) = -\frac{1}{4\pi}\int \nabla’ \times \left(\frac{\bm{F}(\bm{r}’)}{|\bm{r}-\bm{r}’|}\right) d^3 r’ + \frac{1}{4\pi}\int \frac{\nabla’ \times \bm{F}(\bm{r}’)}{|\bm{r}-\bm{r}’|} \, d^3 r’ $$

第 1 項は回転の体積積分なので、発散定理のベクトル版($\int \nabla\times \bm{G}\, dV = \oint d\bm{S} \times \bm{G}$ の符号違い)で無限遠の面積分になり、同じ減衰条件のもとで消えます。したがって

$$ \begin{equation} \bm{A}(\bm{r}) = \frac{1}{4\pi}\int \frac{\nabla’ \times \bm{F}(\bm{r}’)}{|\bm{r}-\bm{r}’|} \, d^3 r’ \end{equation} $$

こちらは場の回転を「電流」だと思ったときのビオ・サバール型のベクトルポテンシャルです。同様に

$$ \nabla^2 \bm{A} = -\nabla \times \bm{F} $$

を満たします。

自動的に成り立つクーロンゲージ

この $\bm{A}$ には、うれしいおまけが付いてきます。$\nabla \cdot \bm{A} = 0$ が自動的に成り立つのです。実際、上と同じ手順で

$$ \nabla \cdot \bm{A} = \frac{1}{4\pi}\int (\nabla’ \times \bm{F}) \cdot \nabla_{\bm{r}} \frac{1}{|\bm{r}-\bm{r}’|} \, d^3 r’ = -\frac{1}{4\pi}\int (\nabla’ \times \bm{F}) \cdot \nabla’ \frac{1}{|\bm{r}-\bm{r}’|} \, d^3 r’ $$

積の微分法則で書き換えると、面積分の項と $\nabla’ \cdot (\nabla’ \times \bm{F})$ の項に分かれます。前者は減衰条件で消え、後者は「回転の発散はゼロ」という恒等式でそもそもゼロです。よって $\nabla \cdot \bm{A} = 0$。

電磁気学で「クーロンゲージを取る」と言われると、何か人為的な条件を課しているように感じます。しかしヘルムホルツ分解の立場からは、$\bm{A}$ を「回転を源とするビオ・サバール積分」として作れば $\nabla\cdot\bm{A}=0$ は選ぶまでもなく満たされているわけです。ゲージ選択の自由度は、$\bm{A}$ に任意の勾配場 $\nabla \chi$ を足しても $\nabla\times\bm{A}$ が変わらないことから来ますが、$\nabla\cdot\bm{A}=0$ を要求すれば $\nabla^2\chi=0$ となり、減衰条件のもとで $\chi$ は定数に限られます。つまり$\bm{A}$ はゲージまで込めて一意に決まるのです。

ここまでで「分解は存在する」ことが示せました。しかし、もし別のレシピで作った $\phi’, \bm{A}’$ が同じ $\bm{F}$ を再現したらどうでしょう。分解が「複数ある」のなら、それは物理量として使いものになりません。次はこの心配を潰しにいきます。

分解の一意性

主張はこうです。$\bm{F}$ が無限遠で十分速く減衰するなら、$\bm{F} = \bm{F}_L + \bm{F}_T$($\nabla\times\bm{F}_L=\bm{0}$、$\nabla\cdot\bm{F}_T=0$、両者とも無限遠で減衰)という分解は一通りしかない。

証明は背理法ではなく、差を取る素直な議論です。二つの分解

$$ \bm{F} = \bm{F}_L + \bm{F}_T = \bm{F}_L’ + \bm{F}_T’ $$

があったとして、差を $\bm{G} = \bm{F}_L – \bm{F}_L’$ と置きます。移項すると同時に $\bm{G} = \bm{F}_T’ – \bm{F}_T$ でもあります。ここが議論の急所です。$\bm{G}$ は「二つの回転なし場の差」なので回転がゼロ、同時に「二つの発散なし場の差」なので発散もゼロ。つまり

$$ \nabla \times \bm{G} = \bm{0}, \qquad \nabla \cdot \bm{G} = 0 $$

を同時に満たします。ここで二重回転の公式を $\bm{G}$ に当てると

$$ \nabla^2 \bm{G} = \nabla(\nabla\cdot\bm{G}) – \nabla\times(\nabla\times\bm{G}) = \bm{0} $$

となり、$\bm{G}$ の各デカルト成分は $\mathbb{R}^3$ 全体で調和関数だと分かります。

ここでリウヴィルの定理(調和関数版)が効きます。全空間で調和かつ有界な関数は定数に限られる、という定理です。$\bm{G}$ は減衰する場の差なので無限遠でゼロに近づき、有界です。したがって $\bm{G}$ は定数、しかも無限遠でゼロなのでその定数はゼロ。よって $\bm{G} = \bm{0}$、すなわち $\bm{F}_L = \bm{F}_L’$ かつ $\bm{F}_T = \bm{F}_T’$ です。

もう少し手触りのある別証明も紹介しておきます。$\bm{G}$ は回転なしなので $\bm{G} = \nabla \psi$ と書け、発散なしなので $\nabla^2 \psi = 0$。半径 $R$ の球でグリーンの第一恒等式を使うと

$$ \int_{B_R} |\nabla \psi|^2 \, dV = \oint_{S_R} \psi \, \nabla\psi \cdot d\bm{S} – \int_{B_R} \psi \nabla^2\psi \, dV $$

右辺第 2 項は $\nabla^2\psi=0$ でゼロ。第 1 項は、$\psi$ が $1/R$ 程度、$\nabla\psi$ が $1/R^2$ 程度で減衰するなら被積分関数が $1/R^3$、面積が $R^2$ なので $1/R$ で消えます。$R\to\infty$ とすると左辺の $\int |\nabla\psi|^2 dV$ がゼロ、非負量の積分がゼロなので $\nabla\psi = \bm{0}$、つまり $\bm{G}=\bm{0}$ です。こちらは「エネルギーがゼロなら場がゼロ」という物理的に分かりやすい論法になっています。

ここで注意点をひとつ。この一意性は無限領域+減衰条件という設定に強く依存しています。有界領域で同じことをやろうとすると、境界積分が消えないので、$\phi$ や $\bm{A}$ に境界条件(ディリクレ型やノイマン型)を課さないと一意になりません。さらに、穴の空いた領域(トーラス状の領域など)では「回転もゼロ・発散もゼロだが恒等的にゼロでもない」調和場が現れます。これを含めた三分解がホッジ分解で、領域のトポロジーと調和場の次元が結びつく美しい理論に発展します。

一意性証明の論理を段階的に示した流れ図と、穴のあいた環状領域に残る調和場(回転も発散もゼロだがゼロ場ではない)の矢印図

左が証明の流れです。「二つの分解の差」という一手だけで、$\bm{G}$ は回転ゼロと発散ゼロを同時に背負い込み、そこから調和関数へ、さらにリウヴィルの定理でゼロへ追い込まれます。右は一意性が崩れる例で、穴の外側を回る場 $\bm{G} = (-y, x)/\rho^2$ を描いたものです。この場の発散と回転を中心差分で実測すると、環状領域上でどちらも最大 $1.5\times10^{-9}$(差分の丸め誤差レベル)でした。ゼロ場ではないのに回転も発散もゼロ——穴があるせいで「$\bm{G}=\nabla\psi$ と書いて $\psi$ を一価に取る」議論が使えなくなり、この場が第 3 の成分として残るのです。

さて、ここまでの導出は積分公式が主役でした。式は正確ですが、直感はまだ育っていません。実はフーリエ変換に移ると、ヘルムホルツ分解は単なるベクトルの射影になります。次のセクションで、その驚くほど単純な姿を見ましょう。

波数空間で見るヘルムホルツ分解 — 縦波射影と横波射影

微分が掛け算になる

フーリエ変換を

$$ \hat{\bm{F}}(\bm{k}) = \int \bm{F}(\bm{r}) e^{-i \bm{k}\cdot\bm{r}} \, d^3 r $$

と定義します。指数関数の微分は指数関数ですから、$\partial/\partial x_j$ は $i k_j$ の掛け算に化けます。ベクトル記号でまとめると、$\nabla \to i\bm{k}$ という置き換えです。すると発散と回転は

$$ \widehat{\nabla \cdot \bm{F}} = i\bm{k}\cdot\hat{\bm{F}}, \qquad \widehat{\nabla \times \bm{F}} = i\bm{k}\times\hat{\bm{F}} $$

になります。微分方程式が代数方程式に化けるのがフーリエ変換の御利益ですが、ここでの御利益は代数を超えています。上の 2 式を眺めると、回転なし条件と発散なし条件が、$\hat{\bm{F}}(\bm{k})$ と $\bm{k}$ の幾何学的な関係に翻訳されているのです。

  • $\nabla \times \bm{F} = \bm{0} \iff \bm{k}\times\hat{\bm{F}} = \bm{0} \iff \hat{\bm{F}} \parallel \bm{k}$(波数ベクトルに平行
  • $\nabla \cdot \bm{F} = 0 \iff \bm{k}\cdot\hat{\bm{F}} = 0 \iff \hat{\bm{F}} \perp \bm{k}$(波数ベクトルに垂直

つまり波数空間では、「回転なし」=「$\bm{k}$ 方向を向いている」、「発散なし」=「$\bm{k}$ に直交している」という、ただそれだけの話になります。

分解はベクトルの直交分解そのもの

3 次元空間の任意のベクトルは、与えられた方向 $\hat{\bm{k}} = \bm{k}/|\bm{k}|$ に平行な成分と垂直な成分に一意に分けられます。中学校の力の分解と同じです。この当たり前の事実を各波数 $\bm{k}$ ごとに実行すれば、それがヘルムホルツ分解になります。

$$ \hat{\bm{F}}(\bm{k}) = \underbrace{(\hat{\bm{k}}\cdot\hat{\bm{F}})\hat{\bm{k}}}_{\hat{\bm{F}}_L} + \underbrace{\hat{\bm{F}} – (\hat{\bm{k}}\cdot\hat{\bm{F}})\hat{\bm{k}}}_{\hat{\bm{F}}_T} $$

射影行列で書くとさらにすっきりします。

$$ \begin{equation} \bm{P}_L(\bm{k}) = \frac{\bm{k}\bm{k}^{\mathsf T}}{|\bm{k}|^2}, \qquad \bm{P}_T(\bm{k}) = \bm{I} – \frac{\bm{k}\bm{k}^{\mathsf T}}{|\bm{k}|^2} \end{equation} $$

$\bm{P}_L$ は縦波射影(longitudinal projector)、$\bm{P}_T$ は横波射影(transverse projector)と呼ばれます。成分で書けば $(P_L)_{ij} = k_i k_j/|\bm{k}|^2$、$(P_T)_{ij} = \delta_{ij} – k_i k_j/|\bm{k}|^2$ です。

射影行列としての性質もその場で確認できます。$\hat{\bm{k}}$ が単位ベクトルであることから

$$ \bm{P}_L^2 = \hat{\bm{k}}(\hat{\bm{k}}^{\mathsf T}\hat{\bm{k}})\hat{\bm{k}}^{\mathsf T} = \hat{\bm{k}}\hat{\bm{k}}^{\mathsf T} = \bm{P}_L $$

と冪等性が出ます。同様に $\bm{P}_T^2 = \bm{P}_T$、そして

$$ \bm{P}_L \bm{P}_T = \bm{P}_L(\bm{I} – \bm{P}_L) = \bm{P}_L – \bm{P}_L^2 = \bm{0} $$

で二つの射影は互いに直交します。さらに $\bm{P}_L + \bm{P}_T = \bm{I}$、どちらも対称行列。3 次元では $\bm{P}_L$ のランクが 1、$\bm{P}_T$ のランクが 2 で、これは「縦波の偏光は 1 通り、横波の偏光は 2 通り」という物理の事実に対応します。

波数ベクトルに平行な縦成分と垂直な横成分へのベクトル分解の幾何図と、射影行列の冪等性・直交性を乱数で実測した棒グラフ

左の図が波数空間で起きていることのすべてです。$\hat{\bm{F}}(k)$ を $\bm{k}$ に平行な向きと垂直な向きに下ろすだけ——高校で習う力の分解と同じ操作を、波数ごとに実行しているにすぎません。右は射影行列の性質を乱数 2000 個の $\bm{k}$ で実測した結果で、$\bm{P}_L^2 – \bm{P}_L$ や $\bm{P}_L\bm{P}_T$ の成分は最大でも $4.4\times10^{-16}$、$\bm{P}_L + \bm{P}_T – \bm{I}$ に至っては厳密にゼロでした。冪等性・相互直交性・完全性が、どの方向の $\bm{k}$ でも例外なく成り立っています。

積分公式との整合

波数空間の話が、前セクションの積分公式と同じものであることを確かめておきましょう。$\nabla^2 \phi = -\nabla\cdot\bm{F}$ をフーリエ変換すると

$$ -|\bm{k}|^2 \hat{\phi} = -i\bm{k}\cdot\hat{\bm{F}} \quad\Longrightarrow\quad \hat{\phi} = \frac{i\,\bm{k}\cdot\hat{\bm{F}}}{|\bm{k}|^2} $$

です。これを使って $\bm{F}_L = -\nabla\phi$ を変換すると、$-\nabla \to -i\bm{k}$ より

$$ \hat{\bm{F}}_L = -i\bm{k}\hat{\phi} = -i\bm{k}\cdot\frac{i(\bm{k}\cdot\hat{\bm{F}})}{|\bm{k}|^2} = \frac{\bm{k}(\bm{k}\cdot\hat{\bm{F}})}{|\bm{k}|^2} = \bm{P}_L \hat{\bm{F}} $$

($-i \times i = 1$ を使いました)。積分公式で作った縦成分と、波数空間の射影で作った縦成分は同じものです。$1/|\bm{k}|^2$ という因子が、実空間では $1/(4\pi|\bm{r}-\bm{r}’|)$ という遠隔作用の核に対応しています。ここから、ヘルムホルツ分解が非局所的な操作であることも読み取れます。ある点での縦成分を知るには、場全体の情報が必要なのです。

なぜ「縦」「横」なのか、そして直交性

呼び名の由来は平面波を思い浮かべると分かります。$\bm{F}(\bm{r}) = \bm{a}\, e^{i\bm{k}\cdot\bm{r}}$ という波を考えたとき、振幅ベクトル $\bm{a}$ が進行方向 $\bm{k}$ に平行なら、場は進行方向に沿って伸び縮みします。空気の疎密波(音波)がこのタイプで、縦波です。一方 $\bm{a} \perp \bm{k}$ なら、場は進行方向と垂直に振動します。電磁波や弾性体の S 波がこのタイプで、横波です。真空中の電磁波が横波である理由は、$\nabla\cdot\bm{E}=0$ が $\bm{k}\cdot\hat{\bm{E}}=0$ を強制するから——ヘルムホルツ分解の言葉ではそう説明できます。

進行方向に平行に変位する縦波(疎密ができる)と、垂直に変位する横波を並べて比較した図

上段の縦波では、粒子が進行方向に沿って前後するため、粒子が詰まる「密」の場所と離れる「疎」の場所が交互に生まれます。密度が変化するということは体積が変わる、つまり発散を持つということで、これが縦成分の正体です。下段の横波は変位が進行方向と直交するので、粒子の間隔は変わりません。体積を変えずに形だけがずれるので発散はゼロ——横成分が「非圧縮」と同義になる理由がここにあります。

最後に重要な性質を一つ。縦成分と横成分は $L^2$ 内積の意味で直交します。 パーセバルの定理を使うと

$$ \int \bm{F}_L \cdot \bm{F}_T \, d^3 r = \frac{1}{(2\pi)^3}\int \hat{\bm{F}}_L^{*}\cdot\hat{\bm{F}}_T \, d^3 k = \frac{1}{(2\pi)^3}\int (\bm{P}_L\hat{\bm{F}})^{*}\cdot(\bm{P}_T\hat{\bm{F}}) \, d^3 k $$

被積分関数は $\hat{\bm{F}}^{*}\bm{P}_L^{\mathsf T}\bm{P}_T\hat{\bm{F}} = \hat{\bm{F}}^{*}\bm{P}_L\bm{P}_T\hat{\bm{F}} = 0$($\bm{P}_L$ は対称、$\bm{P}_L\bm{P}_T=\bm{0}$)なので、積分は各波数でゼロです。したがってエネルギーがきれいに二分されます。

$$ \int |\bm{F}|^2 \, d^3 r = \int |\bm{F}_L|^2 \, d^3 r + \int |\bm{F}_T|^2 \, d^3 r $$

ベクトル場に対する「ピタゴラスの定理」です。乱流研究で「圧縮成分のエネルギー」と「ソレノイダル成分のエネルギー」を別々に語れるのは、この直交性のおかげです。

理論はここまでで一通り揃いました。次は具体的な場で手を動かし、分解が本当に思ったとおりに働くかを確かめます。

具体例:手で分解できる場

例 1:点源の場は純粋な縦成分

原点に置いた点電荷が作る電場を考えます。

$$ \bm{F}(\bm{r}) = \frac{\bm{r}}{r^3} = \frac{\hat{\bm{r}}}{r^2} $$

この場の回転はゼロです。実際 $\bm{F} = -\nabla(1/r)$ と書けるので、勾配場の回転はゼロという恒等式から直ちに従います。つまり $\phi = 1/r$、$\bm{A} = \bm{0}$ で、純粋な縦成分です。発散を計算すると $\nabla\cdot\bm{F} = 4\pi\delta^3(\bm{r})$ で、原点の点源だけが源になっています。

例 2:剛体回転は純粋な横成分

角速度 $\bm{\omega}$ で回る剛体の速度場を考えます。

$$ \bm{F}(\bm{r}) = \bm{\omega}\times\bm{r} $$

発散を計算すると、$\bm{\omega} = (0,0,\omega)$ の場合 $\bm{F} = (-\omega y, \omega x, 0)$ ですから

$$ \nabla\cdot\bm{F} = \frac{\partial(-\omega y)}{\partial x} + \frac{\partial(\omega x)}{\partial y} + 0 = 0 $$

でゼロです。一方 $\nabla\times\bm{F} = 2\bm{\omega}$ で、一様な渦度を持ちます。純粋な横成分の代表例です。ただしこの場は無限遠で発散するので、そのままでは一意性の議論の枠外にあります。この点は次の例で修正します。

例 3:減衰する混合場(本記事の主役)

減衰条件を満たし、しかも縦成分と横成分の両方を持つ、扱いやすい 2 次元の例を作ります。$\rho^2 = x^2+y^2$ として

$$ \begin{equation} \bm{F}(x,y) = \big((x – y)\,e^{-\rho^2},\; (y + x)\,e^{-\rho^2}\big) \end{equation} $$

ガウス関数で覆いをかけているので、無限遠では指数的に減衰します。この場は次のように分解できます。

$$ \bm{F}_L = (x\,e^{-\rho^2},\; y\,e^{-\rho^2}), \qquad \bm{F}_T = (-y\,e^{-\rho^2},\; x\,e^{-\rho^2}) $$

$\bm{F}_L$ は原点から放射状に外向き(湧き出し型)、$\bm{F}_T$ は原点まわりを反時計回り(渦型)です。それぞれが本当に回転なし・発散なしであることを確かめます。

$\bm{F}_L$ の回転(2 次元なので $z$ 成分だけ)は、積の微分で

$$ (\nabla\times\bm{F}_L)_z = \frac{\partial}{\partial x}\left(y e^{-\rho^2}\right) – \frac{\partial}{\partial y}\left(x e^{-\rho^2}\right) = y(-2x)e^{-\rho^2} – x(-2y)e^{-\rho^2} = 0 $$

$\bm{F}_T$ の発散は

$$ \nabla\cdot\bm{F}_T = \frac{\partial}{\partial x}\left(-y e^{-\rho^2}\right) + \frac{\partial}{\partial y}\left(x e^{-\rho^2}\right) = 2xy\,e^{-\rho^2} – 2xy\,e^{-\rho^2} = 0 $$

どちらもゼロです。

ポテンシャルも書き下せます。$\phi = \tfrac{1}{2}e^{-\rho^2}$ とすると

$$ -\nabla\phi = -\tfrac{1}{2}\left(-2x, -2y\right)e^{-\rho^2} = (x, y)\,e^{-\rho^2} = \bm{F}_L $$

2 次元では、発散なし成分は流れ関数 $\psi$ を使って $\bm{F}_T = (\partial\psi/\partial y,\, -\partial\psi/\partial x)$ と書けます(これは $\bm{A} = \psi\hat{\bm{z}}$ とした $\nabla\times\bm{A}$ に他なりません)。$\psi = \tfrac{1}{2}e^{-\rho^2}$ を入れると

$$ \left(\frac{\partial\psi}{\partial y}, -\frac{\partial\psi}{\partial x}\right) = \left(-y e^{-\rho^2},\; x e^{-\rho^2}\right) = \bm{F}_T $$

きれいに一致します。この例では偶然 $\phi$ と $\psi$ が同じ関数になりました。

源も計算しておきます。$\nabla\cdot\bm{F} = \nabla\cdot\bm{F}_L$ を直接計算すると

$$ \nabla\cdot\bm{F} = (1-2x^2)e^{-\rho^2} + (1-2y^2)e^{-\rho^2} = (2 – 2\rho^2)e^{-\rho^2} $$

同様に $(\nabla\times\bm{F})_z = (2-2\rho^2)e^{-\rho^2}$ です。どちらも原点付近で正、$\rho=1$ でゼロ、それより外では負になります。「中心で湧き出し、周囲のリングで吸い込み、全体としては収支ゼロ」という構造です。

エネルギーの分配も手計算できます。$\int|\bm{F}_L|^2 dA = \int \rho^2 e^{-2\rho^2} \cdot 2\pi\rho \, d\rho = \pi/4$、横成分も同じく $\pi/4$。合計 $\pi/2$ で、これは $\int|\bm{F}|^2 dA$ と一致します。ピタゴラスの定理が成立しているわけです。

この例なら解析解が完全に分かっているので、数値計算の答え合わせに使えます。次のセクションでは、この場を FFT で分解して、上の解析解をどれだけ正確に再現できるか確かめます。

応用 1:電磁気学のスカラー/ベクトルポテンシャル

マクスウェル方程式のうち、源を含まない 2 本

$$ \nabla\cdot\bm{B} = 0, \qquad \nabla\times\bm{E} = -\frac{\partial\bm{B}}{\partial t} $$

を眺めます。1 本目は「磁場は発散なし」と言っています。ヘルムホルツ分解の言葉では、$\bm{B}$ は縦成分を持たない、つまり純粋な横成分です。だから必ず $\bm{B} = \nabla\times\bm{A}$ と書けます。これがベクトルポテンシャルの存在理由です。

これを 2 本目に入れると

$$ \nabla\times\bm{E} = -\frac{\partial}{\partial t}(\nabla\times\bm{A}) = -\nabla\times\frac{\partial\bm{A}}{\partial t} $$

移項して $\nabla\times\left(\bm{E} + \partial\bm{A}/\partial t\right) = \bm{0}$。回転がゼロの場は勾配で書けるので

$$ \bm{E} + \frac{\partial\bm{A}}{\partial t} = -\nabla\phi \quad\Longrightarrow\quad \bm{E} = -\nabla\phi – \frac{\partial\bm{A}}{\partial t} $$

教科書でおなじみの表式が出ました。電場は「静電ポテンシャルの勾配」という縦成分と、「ベクトルポテンシャルの時間変化」という項の和という構造です。

静的な場合はもっと分かりやすくなります。静電場は $\nabla\times\bm{E}=\bm{0}$ なので純粋な縦成分、静磁場は $\nabla\cdot\bm{B}=0$ なので純粋な横成分。電場と磁場は、ヘルムホルツ分解の二つの部屋にそれぞれ住んでいると言えます。そして源の関係は、本記事で導いた積分公式そのものです。$\nabla\cdot\bm{E} = \rho/\varepsilon_0$ を代入すれば

$$ \phi(\bm{r}) = \frac{1}{4\pi\varepsilon_0}\int \frac{\rho(\bm{r}’)}{|\bm{r}-\bm{r}’|}\, d^3 r’ $$

というクーロンポテンシャル、$\nabla\times\bm{B} = \mu_0\bm{j}$ を代入すれば

$$ \bm{A}(\bm{r}) = \frac{\mu_0}{4\pi}\int \frac{\bm{j}(\bm{r}’)}{|\bm{r}-\bm{r}’|}\, d^3 r’ $$

というビオ・サバールのベクトルポテンシャルになります。ヘルムホルツ分解の一般公式に、マクスウェル方程式が指定する源を差し込んだだけなのです。

クーロンゲージ $\nabla\cdot\bm{A}=0$ が自動的に満たされることは既に見ました。電磁場の量子化や、放射場と近接場の分離を議論するときには、電場そのものを $\bm{E} = \bm{E}_L + \bm{E}_T$ と分解し、「縦電場は瞬時のクーロン相互作用を、横電場は光速で伝わる放射を担う」という描像が使われます。ヘルムホルツ分解は、そこでも基本の道具として働いています。

電磁気が「分解の二つの部屋」を静的に使い分ける例だったのに対し、次に見る流体力学では、分解が毎ステップ実行される計算操作として現れます。

応用 2:流体の圧力射影法(非圧縮条件の投影)

非圧縮ナビエ・ストークス方程式は

$$ \frac{\partial\bm{u}}{\partial t} + (\bm{u}\cdot\nabla)\bm{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\bm{u}, \qquad \nabla\cdot\bm{u} = 0 $$

です。この方程式には数値計算者を困らせる特徴があります。圧力 $p$ に時間微分が付いていないのです。$p$ は「時間発展する変数」ではなく、$\nabla\cdot\bm{u}=0$ を守らせるためのラグランジュ未定乗数の役割をしています。では $p$ をどう決めればよいのか。ヘルムホルツ分解が答えを与えます。

Chorin の分数ステップ法(射影法)は、1 ステップを 2 段階に分けます。まず圧力を無視して速度を進め、中間速度 $\bm{u}^{*}$ を作ります。

$$ \bm{u}^{*} = \bm{u}^{n} + \Delta t \left[-(\bm{u}^n\cdot\nabla)\bm{u}^n + \nu\nabla^2\bm{u}^n\right] $$

移流と粘性だけで進めたので、$\bm{u}^{*}$ は一般に $\nabla\cdot\bm{u}^{*}\neq 0$、つまり湧き出しを持ってしまいます。そこで圧力の勾配で補正します。

$$ \bm{u}^{n+1} = \bm{u}^{*} – \frac{\Delta t}{\rho}\nabla p $$

この式を見てください。「元の場」=「発散なし場」+「勾配場」というヘルムホルツ分解そのものです。$\bm{u}^{n+1}$ が横成分、$(\Delta t/\rho)\nabla p$ が縦成分です。

圧力を決めるには、両辺の発散を取ります。$\nabla\cdot\bm{u}^{n+1}=0$ を要求すると

$$ 0 = \nabla\cdot\bm{u}^{*} – \frac{\Delta t}{\rho}\nabla^2 p \quad\Longrightarrow\quad \nabla^2 p = \frac{\rho}{\Delta t}\nabla\cdot\bm{u}^{*} $$

圧力ポアソン方程式です。中間速度の発散が「源」になっています。これはヘルムホルツ分解でスカラーポテンシャルを求めるポアソン方程式 $\nabla^2\phi = -\nabla\cdot\bm{F}$ と(符号の流儀を除いて)同じ形です。

まとめると、射影法の 1 ステップは

$$ \bm{u}^{n+1} = \bm{P}_T \bm{u}^{*} $$

と書けます。横波射影を毎ステップ当てているだけ。周期境界の直接数値シミュレーション(DNS)では、この射影を FFT で $\bm{P}_T = \bm{I} – \bm{k}\bm{k}^{\mathsf T}/|\bm{k}|^2$ の掛け算として実行します。壁のある領域では、圧力ポアソン方程式を有限差分や有限要素で解くことになりますが、やっていることは同じです。

さらに、この射影がエネルギーを減らす方向にしか働かないことも直交性から分かります。$|\bm{u}^{*}|^2 = |\bm{P}_L\bm{u}^{*}|^2 + |\bm{P}_T\bm{u}^{*}|^2$ なので、射影後のエネルギーは必ず元以下です。「圧力は仕事をしない」という非圧縮流体の性質が、射影の幾何から見えてきます。

2 次元では、横成分を流れ関数 $\psi$ で表す渦度・流れ関数定式化もよく使われます。$\bm{u} = (\partial\psi/\partial y, -\partial\psi/\partial x)$ とすれば非圧縮条件は自動で満たされ、渦度 $\omega = -\nabla^2\psi$ の発展方程式だけを解けばよくなります。圧力が方程式から消えるのは、最初から横成分の空間の中だけで議論しているからです。この見通しも、ヘルムホルツ分解を知っていると当たり前に思えてきます。

理論と応用が揃ったので、いよいよ実際に計算してみましょう。

Python で確かめる

ここでは周期境界の 2 次元領域を用意し、FFT で射影行列を掛けるだけの実装で分解を行います。確かめたいのは次の 3 点です。

  1. 分解した縦成分・横成分が解析解と一致する
  2. 縦成分と横成分が数値的に直交するか(内積がゼロ、エネルギーがピタゴラスで分かれるか)
  3. 圧力射影法として使ったときに、発散が本当に消える

準備と分解の実装

まずは格子と波数を作り、解析例の場 $\bm{F} = ((x-y)e^{-\rho^2}, (y+x)e^{-\rho^2})$ を用意します。

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

N = 256          # 格子点数(各方向)
L = 8.0          # 領域の一辺([-4, 4) を周期境界で扱う)
x = -L / 2 + L * np.arange(N) / N
X, Y = np.meshgrid(x, x, indexing="ij")
E = np.exp(-(X**2 + Y**2))       # ガウス覆い

# 解析例:縦成分(x,y)E と 横成分(-y,x)E の和
Fx = (X - Y) * E
Fy = (Y + X) * E

# 波数格子(2π/L 刻み)
kx = 2 * np.pi * np.fft.fftfreq(N, d=L / N)
KX, KY = np.meshgrid(kx, kx, indexing="ij")
K2 = KX**2 + KY**2
K2[0, 0] = 1.0                   # k=0 のゼロ割回避(後で成分を 0 にする)

ガウス関数は $\rho=4$ で $e^{-16}\approx 10^{-7}$ まで落ちるので、周期境界の影響は無視できます。$\bm{k}=\bm{0}$ の点で $\bm{P}_L$ は定義できません。ここは「場の空間平均」に対応する成分で、定数場は回転もゼロ・発散もゼロという調和成分です。今回は縦成分側にゼロを割り当て、平均を横成分に含める流儀を採ります。

分解本体は 10 行足らずです。

def helmholtz_2d(ux, uy, KX, KY, K2):
    """FFT で 2 次元ベクトル場を縦成分(回転なし)と横成分(発散なし)に分解"""
    uxh = np.fft.fft2(ux)
    uyh = np.fft.fft2(uy)
    kdotu = KX * uxh + KY * uyh          # k · u_hat
    lxh = KX * kdotu / K2                # P_L u_hat の x 成分
    lyh = KY * kdotu / K2
    lxh[0, 0] = 0.0                      # 平均成分(調和成分)は縦から外す
    lyh[0, 0] = 0.0
    lx = np.real(np.fft.ifft2(lxh))
    ly = np.real(np.fft.ifft2(lyh))
    return lx, ly, ux - lx, uy - ly      # 縦成分, 横成分

FLx, FLy, FTx, FTy = helmholtz_2d(Fx, Fy, KX, KY, K2)

# 解析解との比較
FLx_a, FLy_a = X * E, Y * E
FTx_a, FTy_a = -Y * E, X * E
err_L = np.max(np.abs(FLx - FLx_a)) / np.max(np.abs(FLx_a))
err_T = np.max(np.abs(FTx - FTx_a)) / np.max(np.abs(FTx_a))
print(f"縦成分の最大相対誤差: {err_L:.3e}")
print(f"横成分の最大相対誤差: {err_T:.3e}")

出力は 縦成分の最大相対誤差: 2.625e-08横成分の最大相対誤差: 2.625e-08 となります。相対誤差が $10^{-8}$ 台に収まっており、FFT による分解が解析解を正しく再現していることが確認できました。倍精度の丸め誤差($10^{-16}$)まで落ちないのは、ガウス関数が周期領域の端で完全にゼロになっていないための切断誤差です。領域 $L$ を広げると誤差はさらに小さくなります。

3 つの場を並べて見る

元の場・縦成分・横成分を quiver プロットで並べます。

s = 12   # 矢印の間引き
fields = [(Fx, Fy, "元の場 $\\bf{F}$(渦と湧き出しが混在)"),
          (FLx, FLy, "縦成分 $\\bf{F}_L$(回転なし・放射状)"),
          (FTx, FTy, "横成分 $\\bf{F}_T$(発散なし・渦)")]

fig, axes = plt.subplots(1, 3, figsize=(15, 5.2))
for ax, (ux, uy, ttl) in zip(axes, fields):
    ax.quiver(X[::s, ::s], Y[::s, ::s], ux[::s, ::s], uy[::s, ::s],
              np.hypot(ux, uy)[::s, ::s], cmap="viridis", scale=12)
    ax.set_title(ttl, fontsize=12)
    ax.set_xlabel("$x$"); ax.set_ylabel("$y$")
    ax.set_xlim(-3, 3); ax.set_ylim(-3, 3)
    ax.set_aspect("equal")
plt.tight_layout()
plt.show()

FFTで分解した元の場・縦成分・横成分の3枚のベクトル場プロット

3 枚の図を見比べると、分解の意味が一目で分かります。左の元の場は矢印が斜めに渦を巻きながら外へ広がる「らせん状」で、渦と湧き出しのどちらが主役なのか目では判別できません。中央の縦成分は原点から放射状にまっすぐ外を向く純粋な湧き出しパターン、右の横成分は原点を中心に反時計回りにきれいに回る純粋な渦パターンです。らせんは「放射」と「回転」の重ね合わせだったことが、数式ではなく絵として理解できます。

発散と回転がどちらの成分に載るか

次に、それぞれの成分の発散と回転を計算して、分業が成立していることを確かめます。微分はスペクトル法($ik$ の掛け算)で行います。

def ddx(f): return np.real(np.fft.ifft2(1j * KX * np.fft.fft2(f)))
def ddy(f): return np.real(np.fft.ifft2(1j * KY * np.fft.fft2(f)))

div_F  = ddx(Fx)  + ddy(Fy)
curl_F = ddx(Fy)  - ddy(Fx)
div_T  = ddx(FTx) + ddy(FTy)     # 0 になるはず
curl_L = ddx(FLy) - ddy(FLx)     # 0 になるはず

src_a = (2 - 2 * (X**2 + Y**2)) * E   # 解析解:div も curl もこの形
print(f"div F の解析解との最大誤差 : {np.max(np.abs(div_F - src_a)):.3e}")
print(f"curl F の解析解との最大誤差: {np.max(np.abs(curl_F - src_a)):.3e}")
print(f"div(F_T) の最大絶対値      : {np.max(np.abs(div_T)):.3e}")
print(f"curl(F_L) の最大絶対値     : {np.max(np.abs(curl_L)):.3e}")

出力は順に 2.018e-052.018e-053.017e-093.017e-09 です。前半 2 つは「発散も回転も $(2-2\rho^2)e^{-\rho^2}$ になる」という手計算の結果を $10^{-5}$ 精度で再現しています。後半 2 つが重要で、横成分の発散と縦成分の回転が $10^{-9}$ 台、すなわち実質ゼロになりました。「縦成分だけが発散を担い、横成分だけが回転を担う」という分業が数値的に確認できたわけです。

源の空間分布も見ておきましょう。

fig, axes = plt.subplots(1, 2, figsize=(11, 4.6))
for ax, (Z, ttl) in zip(axes, [(div_F, "発散 $\\nabla\\cdot\\bf{F}$(縦成分の源)"),
                               (curl_F, "回転 $(\\nabla\\times\\bf{F})_z$(横成分の源)")]):
    im = ax.pcolormesh(X, Y, Z, cmap="RdBu_r", vmin=-2, vmax=2, shading="auto")
    ax.contour(X, Y, Z, levels=[0.0], colors="k", linewidths=1.0)
    ax.set_title(ttl, fontsize=12); ax.set_aspect("equal")
    ax.set_xlim(-3, 3); ax.set_ylim(-3, 3)
    ax.set_xlabel("$x$"); ax.set_ylabel("$y$")
    fig.colorbar(im, ax=ax)
plt.tight_layout()
plt.show()

発散と回転の空間分布のカラーマップと、径方向プロファイルが解析解と一致する様子

どちらのカラーマップも、中心が赤(正)、$\rho=1$ の黒い等高線でゼロ、その外側が青(負)という同じ同心円構造をしています。中心で湧き出し、外側のリングで吸い込む「ドーナツ型の源」です。$\rho \to \infty$ で急速にゼロへ落ちているので、導出で仮定した減衰条件が満たされていることも見て取れます。右端のパネルは $y=0$ 上での断面で、スペクトル微分で求めた発散が解析解 $(2-2\rho^2)e^{-\rho^2}$ の曲線に完全に重なっていること、そして横成分の発散と縦成分の回転がどちらも $3.0\times10^{-9}$ にとどまることを同時に示しています。

直交性とエネルギー分配

理論が予言した直交性を数値で確かめます。

dA = (L / N) ** 2
inner = np.sum(FLx * FTx + FLy * FTy) * dA
nL = np.sqrt(np.sum(FLx**2 + FLy**2) * dA)
nT = np.sqrt(np.sum(FTx**2 + FTy**2) * dA)
nF = np.sqrt(np.sum(Fx**2 + Fy**2) * dA)
print(f"内積 <F_L, F_T>      : {inner:.3e}")
print(f"規格化した内積        : {inner / (nL * nT):.3e}")
print(f"|F_L| = {nL:.6f},  |F_T| = {nT:.6f},  |F| = {nF:.6f}")
print(f"|F|^2 - (|F_L|^2 + |F_T|^2) = {nF**2 - (nL**2 + nT**2):.3e}")

出力は 内積 <F_L, F_T> : -6.385e-18規格化した内積: -8.129e-18|F_L| = 0.886227, |F_T| = 0.886227, |F| = 1.253314|F|^2 - (|F_L|^2 + |F_T|^2) = 4.441e-16 です。規格化した内積が $10^{-18}$ という倍精度の丸め誤差レベルで、縦成分と横成分は完全に直交しています。ノルムも手計算と一致しており、$|\bm{F}_L| = |\bm{F}_T| = \sqrt{\pi}/2 \approx 0.886227$、$|\bm{F}| = \sqrt{\pi/2} \approx 1.253314$ です。最後の行はピタゴラスの定理の残差で、これも丸め誤差レベル。エネルギーが二つの成分にきれいに折半されている様子が数値でも確認できました。

エネルギーが縦成分と横成分にπ/4ずつ折半される棒グラフと、計算領域を広げるほど分解の相対誤差が下がるグラフ

左の棒グラフでは、縦成分と横成分のエネルギーがともに $0.785398 = \pi/4$ で、その和 $1.570796 = \pi/2$ が元の場のエネルギーと小数第 6 位まで一致しています。もし二成分が直交していなければ交差項が残り、和は元の値からずれるはずです。右のグラフは、計算領域の一辺 $L$ を $4$ から $12$ まで広げたときの相対誤差で、$L=4$ での $1.1\times10^{-3}$ から $L=12$ では $2.6\times10^{-16}$(倍精度の限界)まで 13 桁も改善します。誤差の正体はアルゴリズムではなく、ガウス関数が周期領域の端で完全にゼロになっていないことによる切断誤差だ、と結論できます。

圧力射影法のデモ

最後に、応用 2 で説明した射影法を実演します。「非圧縮な速度場」に「勾配場」を人工的に足したものを中間速度 $\bm{u}^{*}$ だと思い、射影で元の非圧縮場を取り戻せるかを試します。

# 真の非圧縮速度場(渦)と、圧力に相当するスカラー
u_true_x, u_true_y = -Y * E, X * E
q = np.exp(-((X - 1.0)**2 + Y**2))          # これが「圧力」の正解
ustar_x = u_true_x + ddx(q)                  # 中間速度 u* = 非圧縮場 + ∇q
ustar_y = u_true_y + ddy(q)

# 射影(Δt/ρ = 1 の単位系)
Lx, Ly, Tx, Ty = helmholtz_2d(ustar_x, ustar_y, KX, KY, K2)

div_before = ddx(ustar_x) + ddy(ustar_y)
div_after  = ddx(Tx) + ddy(Ty)
print(f"射影前の発散 RMS : {np.sqrt(np.mean(div_before**2)):.3e}")
print(f"射影後の発散 RMS : {np.sqrt(np.mean(div_after**2)):.3e}")
print(f"真の非圧縮場との最大誤差: {np.max(np.abs(Tx - u_true_x)):.3e}")

# 圧力の復元:∇²p = ∇·u*  →  p_hat = -div_hat/|k|²
ph = np.fft.fft2(div_before) / (-K2); ph[0, 0] = 0.0
p = np.real(np.fft.ifft2(ph))
print(f"圧力 p と正解 q の差(平均補正後): {np.max(np.abs((p - p.mean()) - (q - q.mean()))):.3e}")

出力は 射影前の発散 RMS : 4.432e-01射影後の発散 RMS : 1.806e-07真の非圧縮場との最大誤差: 7.688e-09圧力 p と正解 q の差(平均補正後): 2.638e-07 です。射影前は発散の RMS が 0.443 あったのに、射影後は $10^{-7}$ 台まで 6 桁以上落ちました。しかも取り出された横成分は真の非圧縮場と $10^{-9}$ 精度で一致し、圧力ポアソン方程式から復元した $p$ は(定数の不定性を除いて)正解の $q$ と $10^{-7}$ 精度で一致しています。流体ソルバが毎ステップやっているのは、まさにこの操作です。

射影の様子を可視化してみます。

fig, axes = plt.subplots(1, 3, figsize=(15, 5.2))
panels = [(ustar_x, ustar_y, "中間速度 $\\bf{u}^*$(発散あり)"),
          (Lx, Ly, "取り除かれる縦成分 $\\nabla p$"),
          (Tx, Ty, "射影後 $\\bf{u}^{n+1}=P_T\\bf{u}^*$(発散ゼロ)")]
for ax, (ux, uy, ttl) in zip(axes, panels):
    ax.quiver(X[::s, ::s], Y[::s, ::s], ux[::s, ::s], uy[::s, ::s],
              np.hypot(ux, uy)[::s, ::s], cmap="plasma", scale=12)
    ax.set_title(ttl, fontsize=12); ax.set_aspect("equal")
    ax.set_xlim(-3, 3); ax.set_ylim(-3, 3)
    ax.set_xlabel("$x$"); ax.set_ylabel("$y$")
plt.tight_layout()
plt.show()

中間速度・取り除かれる勾配成分・射影後の速度場の3枚と、射影前後の発散RMSの比較

左の中間速度は、$x=1$ 付近に吸い寄せの成分が乗って渦が歪んでいます。中央が取り除かれる勾配成分で、$(1,0)$ に向かって集まる形——「圧力の山 $q$ に向かって登る勾配 $\nabla q$」そのものです(山の頂上ではラプラシアンが負、つまり吸い込みになっています)。右の射影後は、左のパネルにあった歪みが消えて、$(0,0)$ を中心とするきれいな渦に戻っています。右端の棒グラフが示すとおり、発散の RMS は $4.43\times10^{-1}$ から $1.81\times10^{-7}$ へ落ちました。射影が中間速度から「圧力が担うべき分」だけを正確に抜き取っていることが、絵でも数値でも分かります。

これで理論・具体例・数値実験がひととおり揃いました。最後に、この分解がいつ成り立ち、いつ注意が必要かを整理しておきます。

分解が成り立つ条件と、破れるとき

ヘルムホルツ分解は「任意のベクトル場に無条件で使える」というほど万能ではありません。使うときに意識すべき点を整理します。

減衰条件。 導出で無限遠の面積分を捨てた箇所がありました。厳密には、$|\bm{r}|\to\infty$ で $|\bm{F}| = O(|\bm{r}|^{-1-\epsilon})$(ある $\epsilon>0$)程度の減衰があれば十分です。一様な場 $\bm{F} = \bm{c}$(定数)や剛体回転 $\bm{\omega}\times\bm{r}$ はこの条件を満たしません。実際、定数場は「回転なし」とも「発散なし」とも言えるので、分解が一意になりません。定数場は $\bm{F}_L$ にも $\bm{F}_T$ にも入れられる調和成分なのです。数値計算で $\bm{k}=\bm{0}$ の扱いに気を使ったのは、まさにこの成分のためです。

有界領域とホッジ分解。 現実の問題は箱や配管の中で解きます。有界領域 $\Omega$ では、境界積分が消えないので、$\phi$ や $\bm{A}$ に境界条件を課さないと一意になりません。典型的には、$\bm{F}_T$ に「壁を貫かない」条件 $\bm{F}_T\cdot\bm{n}=0$ を課し、$\phi$ にノイマン条件 $\partial\phi/\partial n = -\bm{F}\cdot\bm{n}$ を与えます。さらに領域に穴があると、回転もゼロ・発散もゼロだが恒等的にゼロでもない調和場が現れます。これを第 3 の成分として加えた $\bm{F} = \nabla\phi + \nabla\times\bm{A} + \bm{h}$ がホッジ分解で、$\bm{h}$ の張る空間の次元が領域のトポロジー(ベッチ数)で決まる、という深い結果が知られています。数値解析の分野では、この構造を離散レベルで保つ手法が「離散外微分」「有限要素外積解析(FEEC)」として研究されています。

滑らかさ。 ここまでの導出では $\bm{F}$ が 2 回微分できると暗に仮定していました。実用上は超関数の意味で解釈すれば、不連続を含む場(例:異なる媒質の境界での電場)にも適用できます。数値的には、不連続点でギブス振動が出るので、スペクトル法よりも有限体積・有限要素で扱ったほうが安定です。

分解は非局所的。 これは概念として重要です。ある点 $\bm{r}$ での縦成分 $\bm{F}_L(\bm{r})$ を求めるには、積分公式が示すとおり場全体の情報が要ります。近傍だけを見て「ここは縦、ここは横」と決めることはできません。だから並列計算では、射影のたびに全域通信(FFT や大域的なポアソンソルバ)が必要になり、大規模流体シミュレーションのボトルネックになります。「圧力は無限の速さで伝わる」という非圧縮近似の性質は、この非局所性の別表現です。

以上を踏まえれば、ヘルムホルツ分解は非常に広い場面で安心して使える道具です。

まとめ

本記事では、ヘルムホルツ分解について、存在の構成的証明から一意性、波数空間での姿、そして応用と数値検証まで通して解説しました。

  • 主張:適当に減衰するベクトル場は必ず $\bm{F} = -\nabla\phi + \nabla\times\bm{A}$ と、回転なし成分(縦成分)と発散なし成分(横成分)に分解でき、その分解は一通りしかない
  • 仕掛け:$\nabla\times(\nabla\phi)=\bm{0}$ と $\nabla\cdot(\nabla\times\bm{A})=0$ により、発散を測ると $\phi$ しか見えず、回転を測ると $\bm{A}$ しか見えない。だから分業が成立する
  • 構成法:$\nabla^2(1/r)=-4\pi\delta^3(\bm{r})$ を使って $-\nabla^2\bm{W}=\bm{F}$ を解き、二重回転の公式に代入すれば $\phi=\nabla\cdot\bm{W}$、$\bm{A}=\nabla\times\bm{W}$ が得られる。源で書き直すと $\nabla^2\phi=-\nabla\cdot\bm{F}$、$\nabla^2\bm{A}=-\nabla\times\bm{F}$ というポアソン方程式のクーロン型・ビオサバール型の解になる
  • 一意性:二つの分解の差は回転も発散もゼロ、よって各成分が調和関数。減衰条件+リウヴィルの定理でゼロに落ちる
  • 波数空間:$\nabla\to i\bm{k}$ の置き換えで、回転なし=$\hat{\bm{F}}\parallel\bm{k}$、発散なし=$\hat{\bm{F}}\perp\bm{k}$。分解は射影行列 $\bm{P}_L=\bm{k}\bm{k}^{\mathsf T}/|\bm{k}|^2$、$\bm{P}_T=\bm{I}-\bm{P}_L$ の適用に帰着し、$L^2$ 直交性とエネルギーのピタゴラス分解が従う
  • 応用:電磁気学のスカラー/ベクトルポテンシャル(クーロンゲージは自動成立)、流体の圧力射影法($\bm{u}^{n+1}=\bm{P}_T\bm{u}^{*}$ と圧力ポアソン方程式)
  • 数値検証:FFT による分解が解析解を相対誤差 $10^{-8}$ で再現し、縦横成分の内積は $10^{-18}$、射影後の発散 RMS は $4\times10^{-1}$ から $2\times10^{-7}$ へ 6 桁低下した

「場を発散源と渦源に分けて考える」という発想は、いったん身につくと、電磁気・流体・弾性波・地球物理・画像処理のあちこちで同じ骨格が見えるようになります。マクスウェル方程式の 4 本が「$\bm{E}$ と $\bm{B}$ の発散と回転を全部指定している」形をしているのも、ヘルムホルツの定理の視点では「場を一意に決めるのに必要十分な情報を並べている」と読めます。

次のステップとして、以下の記事も参考にしてください。