3Dで学ぶ電磁気学

G-7 相補性と保証付き誤差評価

有限要素解が得られても、その誤差は真の解を知らなければ直接計算できません。Bossavitの相補性は、curl条件を満たす場とdiv条件を満たす場を独立に構成し、材料則の不一致から誤差を測る考え方です。この回では、大域的に許容な二つの場からエネルギー誤差と磁気抵抗の上下界を作り、「保証付き評価」と呼べる条件を明確にします。

研究室所蔵の亀有資料「Adaptive Meshing with Bossavit's Error Estimator」とも照合しています。保証付きの相補評価と、局所補間による実用的な回復型指標は成立条件が異なります。後者を適応細分へ使う手順はG-8で扱います。

まず30秒でつかむ

真の解を表示できなくても、Maxwell方程式の異なる半分を厳密に満たす二つの近似場は作れます。一方は $\nabla\times\vec H_h=\vec0$、もう一方は $\nabla\cdot\vec B_h=0$ を満たします。二つが同じ材料則 $\vec B=\mu\vec H$ でも結ばれれば正解です。したがって、構成則のずれが「二つの解がまだ一致していない量」になります。

\nabla\times\vec H_h=\vec0, \qquad \nabla\cdot\vec B_h=0, \qquad \vec e_h=\vec B_h-\mu\vec H_h

大域的な境界条件、循環、総磁束まで満たす許容場を使えば、$\vec e_h$ のエネルギーノルムから保証付き誤差と上下界を得られます。近傍要素の平均だけで場を滑らかにした場合は、細分位置を選ぶ指標にはなっても、それだけで保証付きとは呼べません。

式を読む順番

1. 二つの許容集合を定める。 $\vec H_h$ 側にはcurl条件と循環を、$\vec B_h$ 側にはdiv条件と総磁束を課します。

2. 各集合でエネルギーを最小化する。 一方の保存則を構成上満たす二つの近似場を独立に求めます。

3. 材料則残差を作る。 $\vec e_h=\vec B_h-\mu\vec H_h$ を計算し、どこで二つの場が食い違うかを見ます。

4. 大域許容性を確認してから誤差へ読む。 条件がそろったときだけ相補恒等式と磁気抵抗の上下界を使います。

5. 局所指標と保証を区別する。 局所回復場はG-8のMARKには使えますが、境界・位相条件の検査なしに厳密な上下界とはしません。

この回の到達点

1. curl側とdiv側の許容場が、それぞれどのMaxwell方程式を厳密に満たすか説明できる。

2. 構成則残差から相補誤差を作り、エネルギー誤差との関係を示せる。

3. 磁気抵抗の下界と上界を、二つの変分問題から読み取れる。

4. 大域的に許容な相補場と、局所回復だけで作る誤差指標を区別できる。

修士課程の電磁気学演習

真値を知らずに、解を上下から挟めるか

電磁気の問い: 同じ磁気回路について、curl条件を満たす磁界とdiv条件を満たす磁束を別々に構成する。真の解を表示せずに、計算結果の信頼性をどこまで保証できるだろうか。

観察ポイント: メッシュを細かくしたとき、磁気抵抗の下界、上界、構成則残差がそれぞれどちらへ動くかを書く。材料コントラストを上げた場合も確かめる。

操作と観測: 可視化 G-7-2で試行空間と透磁率比を変え、最初は真値を隠したまま二つの許容解と上下界だけで収束を判断し、最後に真値を表示して判定を検証する。

判定基準と微分幾何による説明

curl側とdiv側の許容集合で得た二つの場が大域境界条件・循環・総磁束を満たせば、構成則のずれ $\vec e_h=\vec B_h-\mu\vec H_h$ が相補ギャップを与え、真の磁気抵抗 $\mathcal R$ を $\mathcal R_h^-\le\mathcal R\le\mathcal R_h^+$ と上下から挟める。

診断: 下界が上界を超えた場合は、どちらかの場が大域的に許容でない。保存則残差が零でない場を使った値は局所指標にはなっても保証付き上下界ではない。材料則が非線形なら線形の二次エネルギー式をそのまま使わない。

研究へ: 商用・自作FEMの解を真値なしで検証するため、保存則を独立に満たす相補場と、上下界幅 $\mathcal R_h^+-\mathcal R_h^-$ を停止条件にした適応解析を構築する。

可視化 G-7-1 二つの許容場

この回の相補恒等式と磁気抵抗の上下界は、線形で対称・正定値な透磁率 $\mu$ を仮定します。非線形材料では、$\mu|\vec H|^2/2$ ではなく凸な磁気エネルギーとその双対を使う別の定式化が必要です。磁気静解析の正解 $(\vec H,\vec B)$ は、位相条件と境界条件を含めて

\nabla\times\vec H=\vec0, \qquad \nabla\cdot\vec B=0, \qquad \vec B=\mu\vec H

を満たします。電流のない領域を想定していますが、G-6の独立循環は許容集合へ含めます。スカラーポテンシャル側では

\vec H_h=-\nabla\Phi_h, \qquad \nabla\times\vec H_h=\vec0

が構成上成り立ちます。intrinsicには $\widetilde\Phi_h$ は磁界 $\widetilde{\mathcal H}_h$ と同じtwisted系列に属するtwisted 0-formです。本教材では右手系を固定し、ordinaryな節点スカラー $\Phi_h$ と1-form $\mathcal H_h$ として離散化します。領域内部ではこの通常の有限要素量だけを使い、向きを反転する写像で別の右手系領域と対応させる場合だけ $\Phi_h=s_g g^*\Phi'_h$ の符号表を参照します。これに対し、磁気ベクトルポテンシャル $\mathcal A_h$ と磁束密度 $\mathcal B_h=d\mathcal A_h$ はintrinsicにもstraight系列です。

ベクトルポテンシャル側では

\vec B_h=\nabla\times\vec A_h, \qquad \nabla\cdot\vec B_h=0

が構成上成り立ちます。Bossavitの記法に合わせ、前者をdiv side、後者をcurl sideと呼びます。どちらも一方のMaxwell方程式には適合していますが、有限次元空間では一般に $\vec B_h\ne\mu\vec H_h$ です。

許容集合を $\mathbb H^I$ と $\mathbb B^F$ とすれば、二つの場は相補的なエネルギー最小化として書けます。

\vec H_h =\operatorname*{arg\,min}_{\vec h\in\mathbb H_h^I} \dfrac12\int_D\mu|\vec h|^2\,dV

\vec B_h =\operatorname*{arg\,min}_{\vec b\in\mathbb B_h^F} \dfrac12\int_D\mu^{-1}|\vec b|^2\,dV

可視化 G-7-1では、青がcurl-freeな $\vec H_h$、オレンジがdiv-freeな $\nu\vec B_h$ です。重要なのは、近似段階を進めてもこの二つの保存則は最初から厳密に満たされており、減るのは第三の量、構成則残差 $\vec r_h=\nu\vec B_h-\vec H_h$ だという点です。

機構を数値で閉じるため、Canvasでは厚さ方向に一様な正規化モデルの $xy$ 代表断面 $D_2=[-1,1]^2$ を使います。$a=0.72$、$\varepsilon_m=0.115/(1+0.72m)$ として

\vec H_h^{(m)} =\begin{pmatrix} 1-\varepsilon_m\pi\cos(\pi x)\sin(\pi y)\\ -\varepsilon_m\pi\sin(\pi x)\cos(\pi y) \end{pmatrix}, \qquad \nu\vec B_h^{(m)} =\begin{pmatrix} 1+a\varepsilon_m\pi\cos(\pi x)\cos(\pi y)\\ a\varepsilon_m\pi\sin(\pi x)\sin(\pi y) \end{pmatrix}

と置きます。この族では、すべての $m$ で保存則が成り立ちます。

\partial_x H_{h,y}^{(m)}-\partial_y H_{h,x}^{(m)}=0, \qquad \partial_x B_{h,x}^{(m)}+\partial_y B_{h,y}^{(m)}=0

一方、断面平均RMS構成則残差は解析的に

\rho_m =\left(\dfrac1{|D_2|}\int_{D_2}|\vec r_h^{(m)}|^2\,dA\right)^{1/2} =\pi\varepsilon_m\sqrt{\dfrac{1+a^2}{2}}

です。固定比較 $m:1\to5$ では $\rho_m:0.1830\to0.0684$、基準比は $0.374$ になります。これは完全なFEM問題の次元付き $\eta_{\mathrm{comp}}$ そのものではなく、保存則と構成則残差を分けて読むための正規化診断モデルです。

近似段階 $m$ 基準 m=1 / RMS||r_h||=0.1830 / curl H_h=0 / div B_h=0

確認：$m:1\to5$ で何が小さくなる？二つの保存則と、二つを結ぶ材料則を別々に考えてください。

動かして確かめる順序（操作して答える問い）。 初期状態では $m=1$ に固定します。左の青矢印は勾配から作った $\vec H_h$、中央の橙矢印は回転から作った $\nu\vec B_h$、右の赤矢印と濃淡は差 $\vec r_h$ です。$m=5$ の結果を図を動かして判定を開き、右欄に残る灰色の $m=1$ 基準と赤い現在値を同じ尺度で比較してください。$\rho:0.1830\to0.0684$ と減っても、$\nabla\times\vec H_h=0$ と $\nabla\cdot\vec B_h=0$ はどちらも $0\to0$ です。操作後に $m=0,6$ も動かし、近づく理由は保存則を緩めたからではなく、各許容集合の中で構成則を同時に満たす点へ近づいたからだと説明してください。

可視化 G-7-2 相補性と磁気抵抗の上下界

二つの許容場の構成則残差を、相補誤差と定義します。

\eta_{\mathrm{comp}}^2 =\int_D\mu^{-1} \left|\vec B_h-\mu\vec H_h\right|^2\,dV

境界条件、独立循環、総磁束を含めて両方が大域的に許容なら、交差項が消え、真の解に対する二つのエネルギー誤差へ分解できます。

\eta_{\mathrm{comp}}^2 =\left\|\vec B_h-\vec B\right\|_{\mu^{-1}}^2 +\left\|\vec H_h-\vec H\right\|_{\mu}^2

ここで

\left\|\vec v\right\|_{\mu}^2 =\int_D\mu|\vec v|^2\,dV, \qquad \left\|\vec w\right\|_{\mu^{-1}}^2 =\int_D\mu^{-1}|\vec w|^2\,dV

です。起磁力 $I$ を与えるdiv sideと、総磁束 $F$ を与えるcurl sideを同じ磁気回路へ適用すると、磁気抵抗 $R=I/F$ を両側から挟めます。

R_{\mathrm{div},h} =\dfrac{I^2}{\displaystyle\int_D\mu|\vec H_h|^2\,dV}, \qquad R_{\mathrm{curl},h} =\dfrac{\displaystyle\int_D\mu^{-1}|\vec B_h|^2\,dV}{F^2}

R_{\mathrm{div},h}\le R\le R_{\mathrm{curl},h}

可視化 G-7-2では、単に幅が縮む演出は使いません。3×3節点・12枝の磁気回路をその場で解き、磁位側と保存磁束側のエネルギーから下界と上界を計算します。近似を改善して、実際に計算された区間が狭くなったときだけ、真値を知らなくても計算の確からしさが上がったと判断します。

二つの変分原理で磁気抵抗を挟む

最初のイメージ。 答えそのものが分からなくても、必ず小さすぎる見積りと必ず大きすぎる見積りを別々に作り、両側から挟むことができます。磁気回路では、磁位を試すDirichlet側と、保存則を満たす磁束を試すThomson側がその二本の物差しです。二本が近づけば、厳密解を表示しなくても誤差幅を判断できます。

枝と節点のつながりを接続行列 $\boldsymbol D$、枝磁気抵抗を対角行列 $\boldsymbol{\mathcal R}$、枝パーミアンスを $\boldsymbol{\mathcal P}=\boldsymbol{\mathcal R}^{-1}$ とします。節点磁位 $\boldsymbol\Phi$ から得る枝起磁力差と、枝磁束 $\boldsymbol f$ の保存則は

\boldsymbol g=\boldsymbol D^{\mathsf T}\boldsymbol\Phi, \qquad \boldsymbol D\boldsymbol f=\boldsymbol q

です。$\boldsymbol D$ が「どことどこがつながるか」、$\boldsymbol{\mathcal R}$ と $\boldsymbol{\mathcal P}$ が「材料と形状によって磁束がどれだけ通りやすいか」を担います。起磁力差を1に固定したDirichlet原理では、有効パーミアンス $C$ が

C =\min_{\Phi_s-\Phi_t=1} (\boldsymbol D^{\mathsf T}\boldsymbol\Phi)^{\mathsf T} \boldsymbol{\mathcal P} (\boldsymbol D^{\mathsf T}\boldsymbol\Phi), \qquad R=\dfrac1C

で決まります。有限個の試行関数だけを使うと最小値を十分に下げられないので、試行エネルギー $E_{H,h}$ は $C$ 以上です。したがって、その逆数は磁気抵抗の下界になります。

E_{H,h}\ge C, \qquad R_{\mathrm{lower},h}=\dfrac1{E_{H,h}}\le R

一方、単位磁束を端子間へ流すThomson原理では、節点保存則を満たすすべての枝磁束から損失を最小にします。

R =\min_{\boldsymbol D\boldsymbol f=\boldsymbol q} \boldsymbol f^{\mathsf T}\boldsymbol{\mathcal R}\boldsymbol f

保存則を壊さない閉路基底を一部だけ使うと、流れの選択肢が少ないためエネルギーは真値以上です。これが上界です。

R_{\mathrm{upper},h} =\boldsymbol f_h^{\mathsf T}\boldsymbol{\mathcal R}\boldsymbol f_h \ge R, \qquad R_{\mathrm{lower},h}\le R\le R_{\mathrm{upper},h}

近似段階 $m$ を上げると、左の節点磁位空間と右の保存磁束空間が広がり、二つの境界が近づきます。$\mu_r$ を上げると緑の枝の磁気抵抗が下がり、オレンジの磁束がそこへ再配分されます。「厳密参照値」は答え合わせ専用です。まず非表示のまま区間幅と保存則残差だけで結果を判断し、最後に表示して確認してください。

試行空間の段階 $m$ 緑の枝の比透磁率 $\mu_r$ 答え合わせ厳密参照値を表示離散磁気回路の下界と上界

$m$ を1段上げると区間幅は？真値を表示せず、二つの試行空間が広がる効果を図を動かして確かめてください。

節点色 $\Phi_h$ → 磁位エネルギー $E_H$ → 緑の下界 $R_L=1/E_H$ 橙矢印の保存磁束 $\boldsymbol f_h$ → 磁束エネルギー $E_B$ → 橙の上界 $R_U=E_B$

動かして確かめる順序。 まず厳密参照値を隠します。色付き節点 $\Phi_h$ から得る $E_H$ が緑の下界 $R_L=1/E_H$ を作り、保存則を満たす橙矢印 $\boldsymbol f_h$ から得る $E_B$ が橙の上界 $R_U=E_B$ を作ることを図中の対応で追ってください。次に、$m$ を1段上げたときの区間幅を選択欄で図を動かして計算します。緑枝の $\mu_r$ を上げる前にも磁束がどちらへ再配分されるかを確認し、枝色と上下界で照合します。最後だけ参照値を表示し、真値を知らなくても区間幅から品質を判断できたかを確認してください。

右手系代表を使ってもFEMの物理解は変わらない

intrinsicなtwisted formを、各領域で右手系のordinary formとして表すと、向きを反転する領域写像ではtwisted量の係数に符号が付きます。しかし、これは未知量の基底を符号付きで取り替えたものであり、界面の対応を同じ規約で組み立てれば物理場は変わりません。

符号付き自由度対応を対角行列 $\boldsymbol S$ で表し、元の係数を $\boldsymbol u=\boldsymbol S\widehat{\boldsymbol u}$ と書くと、変換後の有限要素方程式は

\widehat{\boldsymbol K} =\boldsymbol S^{\mathsf T}\boldsymbol K\boldsymbol S, \qquad \widehat{\boldsymbol f} =\boldsymbol S^{\mathsf T}\boldsymbol f

です。$\widehat{\boldsymbol K}\widehat{\boldsymbol u}=\widehat{\boldsymbol f}$ を解いて $\boldsymbol u=\boldsymbol S\widehat{\boldsymbol u}$ と戻せば、元の有限要素解と一致します。特に

\dfrac12 \widehat{\boldsymbol u}^{\mathsf T} \widehat{\boldsymbol K} \widehat{\boldsymbol u} = \dfrac12 \boldsymbol u^{\mathsf T} \boldsymbol K \boldsymbol u

なので、磁気エネルギー、磁束密度、磁界、力などの物理量は同じです。異なって見える可能性があるのは、twisted量を右手系代表で保存した生の係数の符号だけです。外部領域と変換領域を同じソルバーで解く場合は、界面に $\boldsymbol u_{\mathrm{ext}}=\boldsymbol S\boldsymbol u_{\mathrm{int}}$ を課すか、変換領域の未知量へ符号を吸収して通常の等値拘束にします。

保証付き相補評価と局所指標を区別する

完全な相補解を毎回求める代わりに、計算済みの場を別の適合空間へ局所補間して構成則のずれを測ることもできます。ただし、局所回復場が境界条件、全体の循環、総磁束まで満たすとは限りません。

\text{global admissibility} \quad\Longrightarrow\quad \text{guaranteed complementary bound}

\text{local recovered field only} \quad\Longrightarrow\quad \text{refinement indicator}

したがって「相補誤差」と「適応細分に使える指標」は同義ではありません。前者は大域的な許容性を検査して初めて保証をもち、後者は誤差が大きそうな要素を選ぶ実用的な手掛かりです。局所補間、MARK、h細分、p次数上昇は次のG-8で実際に動かします。

保存則・検算と次への接続

各近似段階で、接続行列の恒等式と磁気抵抗の挟み込みを独立に検査します。

\boldsymbol C_h\boldsymbol G_h=\boldsymbol0, \qquad \boldsymbol D_h\boldsymbol C_h=\boldsymbol0

R_{\mathrm{div},h} \le R\le R_{\mathrm{curl},h}

非入れ子空間、反復解法誤差、不適切な境界条件・位相制約があると、期待する単調な挟み込みは崩れます。G-3のWhitney形式と離散Hodge、G-6の大域循環を正しく組み立てて初めて、相補性が誤差保証として働きます。

可視化 G-7-3 相補誤差の局所分布

curl側とdiv側から得た許容場の不一致を、要素 $K$ ごとに

\eta_K^2 = \int_K \mu\left| \vec H_h-\nu\vec B_h \right|^2dV

と測ります。大域相補誤差は

\eta_{\mathrm{comp}}^2 = \sum_K\eta_K^2

です。不一致量を増減し、同じ大域値でも誤差が局在する要素を優先して細分すべき理由を3Dで確認します。

局所不一致 $\delta$

可視化 G-7-4 局所回復指標と保証付き相補量を分ける

回復場の局所差だけを測る指標を

\eta_{\mathrm{loc}}^2 =\sum_{K\in\mathcal T_h} \|\vec H_h-\nu\vec B_h^*\|_{K}^{2}

とします。境界条件、切断面、位相制約まで満たす大域許容場でなければ、これは保証付き上界とは呼べません。概念的には、大域許容性の欠陥 $\eta_{\mathrm{bdry}}$ も含めて

\eta_{\mathrm{comp}}^2 =\eta_{\mathrm{loc}}^2+\eta_{\mathrm{bdry}}^2

と監視する必要があります。局所差を小さくしたまま境界欠陥を増やし、「図の中で場が近い」ことと「真の解を挟み込む」ことが別であると確認してください。

局所不一致 $\eta_{\mathrm{loc}}$ 大域許容性欠陥 $\eta_{\mathrm{bdry}}$

可視化 G-7-5 Dörflerマーキングで上下界を縮める

局所相補誤差を大きい順に選び、細分対象と保証付き上下界の幅を同時に確認します。

\sum_{K\in\mathcal M}\eta_K^2\ge\theta\sum_{K\in\mathcal T_h}\eta_K^2

R_{\mathrm{curl},h}-R_{\mathrm{div},h}\longrightarrow0

マーキング閾値 $\theta$

理解の確認

1. $\vec H_h=-\nabla\Phi_h$ と $\vec B_h=\nabla\times\vec A_h$ は、それぞれ何を厳密に満たすか。

2. $\eta_{\mathrm{comp}}=0$ なら、二つの許容場について何が言えるか。

3. 右手系代表の係数に符号差があっても、物理場とエネルギーが同じになる理由を説明してください。

4. 回復型指標を保証付き誤差上界として使う前に、どの大域条件を確認すべきか。

5. 下界と上界の間隔が狭くても、反復解法残差を別に監視すべきなのはなぜか。

参考資料

Alain Bossavit, *Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements*, Academic Press (1998), Chapter 6 and Appendix C.
Akihisa Kameari, “Adaptive Meshing with Bossavit's Error Estimator” (2022), 研究室所蔵のNGSolve実装資料。