3Dで学ぶ電磁気学

G-1 有限要素法の基礎

有限要素法は、偏微分方程式をいきなり巨大な行列へ置き換える手法ではありません。解析領域を小さな要素へ分け、各要素で場を補間し、弱形式から得た局所的な関係を全体へ組み立て、境界条件を入れて解く方法です。このページでは、高橋則雄先生の有限要素法の体系に沿って、磁気ベクトルポテンシャル $A_z$ を未知量とする静磁界問題を入口にします。

このページでつなぐ四つの段階

1　分割と近似連続な場を三角形ごとの一次関数で近似し、要素を細かくしたときの誤差を観察します。

2　補間関数面積座標 $N_1,N_2,N_3$ が節点値から要素内部の値を作る様子を動かします。

3　組立と境界要素係数行列を全体行列へ加算し、固定境界と自然境界の違いを確かめます。

4　未知量の置き場所四面体の節点自由度と辺自由度を比べ、電磁界で辺要素が必要になる理由へ進みます。

可視化 G-1-1 メッシュを細かくすると何が改善するか

単位正方形上の既知の場

A_z(x,y)=\sin(\pi x)\sin(\pi y)

を三角形一次要素で補間します。この場が一定の磁気抵抗率 $\nu$ のPoisson型方程式を満たすとき、右辺は

-\nabla\cdot(\nu\nabla A_z)=J_z, \qquad J_z=2\pi^2\nu\sin(\pi x)\sin(\pi y)

です。ここでは解法誤差を混ぜず、節点上では正しい値を与え、要素内部を一次補間したときの誤差だけを見ます。分割数を増やし、折れた近似面が滑らかな場へ近づくことと、計算量が増えることを同時に確認してください。

一辺の分割数 $n$

図を動かして動かしてください。$h=1/n$ を半分にすると、滑らかな解に対する一次補間の $L^2$ 誤差はおおむね $O(h^2)$、勾配の誤差はおおむね $O(h)$ で減ります。

\lVert A_z-I_hA_z\rVert_{L^2(\Omega)}=O(h^2), \qquad \lVert \nabla(A_z-I_hA_z)\rVert_{L^2(\Omega)}=O(h)

可視化 G-1-2 一次三角形要素の中を面積座標で補間する

三角形の節点を $P_1,P_2,P_3$、節点値を $a_1,a_2,a_3$ とすると、要素内部の近似値は

\[ A_h(x,y)=N_1(x,y)a_1+N_2(x,y)a_2+N_3(x,y)a_3 \]

で与えます。一次三角形では $N_i$ は面積座標であり、

N_i(P_j)=\delta_{ij}, \qquad N_1+N_2+N_3=1

を満たします。評価点 $P$ をドラッグするか、$x,y$ を動かしてください。色付きの面積比が $N_1,N_2,N_3$ となり、節点値の重み付き平均が $A_h(P)$ になります。

評価点 $x$ 評価点 $y$ 節点値 $a_1$ 節点値 $a_2$ 節点値 $a_3$

三角形内部で $\nabla N_i$ は一定です。したがって一次要素の $A_h$ は連続でも、その勾配から求める磁束密度

\vec B_h=\nabla\times(A_h\vec e_z) =\dfrac{\partial A_h}{\partial y}\vec e_x -\dfrac{\partial A_h}{\partial x}\vec e_y

は要素ごとに一定となり、要素境界で跳びます。細分によって磁束密度分布が滑らかになる理由はここにあります。

可視化 G-1-3 要素係数行列を全体へ組み立てる

磁気抵抗率 $\nu$ が要素内で一定のとき、一次三角形要素の係数行列は

K_{ij}^{(e)} =\int_{\Omega_e}\nu\nabla N_i\cdot\nabla N_j\,d\Omega =\dfrac{\nu}{4\Delta_e}(b_ib_j+c_ic_j)

です。局所節点番号を全体節点番号へ写すBoolean行列 $\boldsymbol L_e$ を使えば、全体行列は

\boldsymbol K =\sum_e\boldsymbol L_e^{\mathsf T}\boldsymbol K^{(e)}\boldsymbol L_e, \qquad \boldsymbol K\vec a=\vec f

となります。スライダーで要素を一枚ずつ組み立て、同じ全体節点へ複数要素の寄与が足し込まれる様子を見てください。4要素を組み立てた自然境界の行列では、各行の和が0となり、定数ベクトルが零モードです。

組み立て済み要素数表示する系確認：$P_1$―$P_4$を固定すると自由な未知量はいくつ残るか

自然境界だけでは定数を加えても勾配が変わらないため、行列に零空間が残ります。固定境界を入れると、既知の境界値 $a_1$―$a_4$ は自由な未知量から除かれます。この例で残る自由系は中心節点だけです。全要素を組み立て、中心へ単位荷重 $f_5=1$ を与えると、

\boldsymbol K_{ff}=[K_{55}]=[4], \qquad a_5=K_{55}^{-1}f_5=\dfrac14

となります。確認を確定した後、自然境界と固定境界を切り替え、零モードを持つ5自由度系が正則な1自由度系へ変わることを行列の表示で確かめてください。

\boldsymbol K\vec 1=\vec 0 \quad\text{(自然境界のみ)}, \qquad A_z\vert_{\Gamma_D}=0 \quad\text{(固定境界)}

橋渡し G-1-A1 曲面計量から局所剛性行列へ

F-2-C1 測地円と局所曲率では、平面、円筒、球面、一定負曲率の擬球面局所片を、連続な計量と測地円だけで比較しました。ここから先は有限要素法の担当です。同じ曲面 $\Gamma$ 上のPoisson方程式

-\Delta_\Gamma\phi=f

に試験関数 $w$ を掛け、部分積分すると、自然境界項を別にして弱形式

\int_\Gamma\nabla_\Gamma\phi\cdot\nabla_\Gamma w\,dA =\int_\Gamma fw\,dA

を得ます。曲面要素を参照三角形 $\hat T$ の局所座標 $\boldsymbol\xi=(u,v)$ から写すJacobianを $J$、第1基本形式を $G=J^{\mathsf T}J$ とすると、一次形状関数 $\boldsymbol N=(N_1,N_2,N_3)$ の局所剛性行列は

\boldsymbol K_e =\int_{\hat T} (\nabla_{\!\xi}\boldsymbol N)^{\mathsf T} G^{-1} (\nabla_{\!\xi}\boldsymbol N) \sqrt{\det G}\,du\,dv

です。曲面勾配の内積は $G^{-1}$、面積要素は $\sqrt{\det G}$ が担当します。したがって、Gauss曲率 $K$ を行列へ直接足すのではなく、曲率を含む局所距離の変化が計量 $G$ を通って行列へ入ると読むのが正確です。

F-2-C1と同じ座標を使うと、代表的な計量は次のようになります。

| 曲面 | 局所計量 | FEMでの帰結 |

|---|---|---|

| 平面 | $G=I$ | 基準の局所行列 |

| 円筒 | $G=I$ | 外から曲がって見えても平面と同じ局所行列 |

| 半径 $R$ の球面 | $G=\operatorname{diag}(\cos^2(v/R),1)$ | 要素内で係数と面積重みが変化 |

| 擬球面局所片 | $G=\operatorname{diag}\!\left(\tanh^2U/\tanh^2U_0,\operatorname{sech}^2U/\operatorname{sech}^2U_0\right)$ | 負曲率に対応する異方的な局所変化 |

円筒と平面の一致は「曲率を無視した近似」ではなく、両者が局所等長であることのFEM側の表現です。一方、球面と擬球面では中心で $G(0)=I$ にそろえても、有限な要素内で $G$ が変化するため局所行列は平面からずれます。$h\to0$ では一点の計量へ近づきますが、曲面近似誤差、有限要素近似誤差、数値積分誤差は別々に評価しなければなりません。

形作用素 $S$ も通常の曲面Poisson行列へ入れるのか

第1基本形式だけで定まる上のエネルギーには、$S$ を独立な係数として足しません。形作用素が明示的に現れるのは、シェルの曲げ、曲率依存の表面物理、法線を含む境界条件、形状微分などです。連続な $G,S,K$ の区別はF-2、trial/test・要素積分・組立てはこのG-1が正典です。

一次平面三角形で滑らかな曲率を表せるか

一要素内の法線は一定なので、点ごとの曲率は要素内だけでは回復できません。法線ジャンプ、頂点角度欠損、高次曲面写像、または曲率回復が必要です。ただしLaplace--Beltrami局所行列自体は、各要素の接平面上の勾配と面積から通常どおり組み立てられます。

橋渡し可視化 G-1-B1 合同面を折るとFEM行列と法線結合はどう分かれるか

F-2-C2d 多面体版の驚異の定理では、四枚の合同三角形を辺に沿って折るだけなら、頂点角度欠損は変わらないことを見ました。ここでは同じ一頂点starを、そのまま有限要素の局所行列へ渡します。

頂点から出る二辺を単位長、挟角を $\theta=70^\circ$ とする一枚の三角形 $T_e$ では、面積は

|T_e|=\dfrac12\sin\theta=0.469846310393

です。ヒンジ変数 $d$ を変えても三辺と面積は変わりません。一次形状関数 $N_1,N_2,N_3$ の局所質量行列とLaplace--Beltrami剛性行列を

M_{ij}^{(e)}=\int_{T_e}N_iN_j\,dS, \qquad K_{ij}^{(e)}=\int_{T_e}\nabla_\Gamma N_i\cdot\nabla_\Gamma N_j\,dS

とします。質量行列は

\boldsymbol M^{(e)}=\dfrac{|T_e|}{12} \begin{bmatrix}2&1&1\\1&2&1\\1&1&2\end{bmatrix}

であり、剛性行列も辺の内積、すなわち第1基本形式だけで決まります。したがって合同面を折るヒンジ変形では

\boldsymbol M^{(e)}(d)=\boldsymbol M^{(e)}(d_0), \qquad \boldsymbol K^{(e)}(d)=\boldsymbol K^{(e)}(d_0)

です。ところが、固定した外部磁束密度 $\vec B_{\rm ext}=\vec e_z$ と面法線 $\vec n_e(d)$ を結ぶ局所荷重

f_i^{(n)}(d) =\int_{T_e}N_i\,\vec B_{\rm ext}\cdot\vec n_e(d)\,dS =\dfrac{|T_e|}{3}\,\vec B_{\rm ext}\cdot\vec n_e(d)

は変わります。これは行列計算の矛盾ではありません。$\boldsymbol M^{(e)}$ と $\boldsymbol K^{(e)}$ は面内計量だけを見るのに対し、$f_i^{(n)}$ は面が周囲の3次元空間へどう向いているかを見るからです。

ヒンジ変数 $d$ 任意観察：ヒンジを動かしたときの組合せ

基準 $d_0=35^\circ$ から既定の $d=55^\circ$ へ変形すると、一面目では $\vec B_{\rm ext}\cdot\vec n_1$ が $0.554636574843$ から $0.699790071708$ へ、各局所節点の法線荷重係数が $0.086864649433$ から $0.109597927747$ へ変わります。一方、$\boldsymbol M^{(e)}$ と $\boldsymbol K^{(e)}$ の最大成分差は丸め誤差の範囲で0です。

電磁界問題の全応答も不変なのか

不変とは限りません。面内のスカラーLaplace--Beltrami作用素、合同面のP1補間、局所質量・剛性行列は内在計量だけで決まります。しかし法線磁束、表面電流を3次元ベクトルとして見た向き、放射・散乱、外部コイルとの相互結合、曲面境界のtraceは外在配置に依存します。「内在演算子が同じ」と「電磁界境界値問題の右辺・結合・観測量まで同じ」を分けてください。

なぜ一面の局所行列を比較するのか

四面は合同で接続も固定されているため、同じ局所節点順序を使えば各面の $\boldsymbol M^{(e)},\boldsymbol K^{(e)}$ は一致し、その全体組立ても不変です。ここでは原因を見失わないよう一面の3×3行列を表示しています。法線荷重は各面で向きが異なるため、全体へ組み立てる前に面ごとの符号と法線方向を検査します。

橋渡し可視化 G-1-B2 局所3×3行列を5×5全体行列へどう散らすか

G-1-B1では、一枚の合同三角形に限れば局所質量・剛性行列はヒンジ変形で変わらず、法線荷重は変わることを確認しました。次は四面を同じ五つの全体節点へ結び、局所番号から全体番号へのscatterを実行します。中心を $P_0$、境界を $P_1,P_2,P_3,P_4$ とすると、面の局所節点対応は

I_1=(0,1,2),\quad I_2=(0,2,3),\quad I_3=(0,3,4),\quad I_4=(0,4,1)

です。局所値を取り出すBoolean行列を $\boldsymbol L_e\in\mathbb R^{3\times5}$ と書けば、$m$面までの組立ては

\boldsymbol K^{[m]}(d) =\sum_{e=1}^{m}\boldsymbol L_e^{\mathsf T} \boldsymbol K^{(e)}(d)\boldsymbol L_e, \qquad \vec f^{[m]}(d,\beta) =\sum_{e=1}^{m}\boldsymbol L_e^{\mathsf T}\vec f^{(e)}(d,\beta)

となります。四面をすべて組み立てた剛性行列は

\boldsymbol K= \begin{bmatrix} 2.800830&-0.700208&-0.700208&-0.700208&-0.700208\\ -0.700208& 1.064178&-0.181985&0&-0.181985\\ -0.700208&-0.181985& 1.064178&-0.181985&0\\ -0.700208&0&-0.181985& 1.064178&-0.181985\\ -0.700208&-0.181985&0&-0.181985& 1.064178 \end{bmatrix}

です。各面が合同で接続 $I_e$ も変わらないため、ヒンジ変形後も $\boldsymbol K(d)=\boldsymbol K(d_0)$ です。また自然境界だけなので

\boldsymbol K\vec 1=\vec0, \qquad \operatorname{nullity}(\boldsymbol K)=1

となり、定数零モードが残ります。

一方、外部場の向きを

\vec B_{\rm ext}(\beta) =\cos\beta\,\vec e_x+\sin\beta\,\vec e_z, \qquad \vec f^{(e)} =\dfrac{|T_e|}{3}(\vec B_{\rm ext}\cdot\vec n_e) \begin{bmatrix}1\\1\\1\end{bmatrix}

とします。$\beta=0^\circ$では四面の符号が $+,-,-,+$ となります。$q(d)=|T_e|n_{1x}(d)/3$ と置くと、四面をscatterした右辺は

\vec f(d,0)=q(d) \begin{bmatrix}0&2&0&-2&0\end{bmatrix}^{\mathsf T}

となり、中心 $P_0$ と $P_2,P_4$ では正負の面寄与が相殺します。$\beta=90^\circ$では四面が同符号となり、$r(d)=|T_e|n_{1z}(d)/3$ に対して

\vec f(d,90^\circ)=r(d) \begin{bmatrix}4&2&2&2&2\end{bmatrix}^{\mathsf T}

です。全体行列の不変性と、右辺の加算・相殺は別々に検査しなければなりません。

ヒンジ変数 $d$ 外部場の傾き $\beta$ 組立て済み面数 $m$ 任意観察：$\beta=0^\circ,m=4$で荷重が0になる節点

既定の $d=55^\circ,\beta=20^\circ$ では、面別荷重は $(0.109598,-0.034629,-0.034629,0.109598)$、全体右辺は

\vec f= \begin{bmatrix} 0.149939&0.219196&0.074969&-0.069257&0.074969 \end{bmatrix}^{\mathsf T}

です。さらに、右辺の全成分和と組み立てた面上の総法線流束は

\vec1^{\mathsf T}\vec f^{[m]} =\sum_{e=1}^{m}\int_{T_e}\vec B_{\rm ext}\cdot\vec n_e\,dS

で一致します。Canvasはこの流束収支残差も計算し、見た目だけでscatterを判定しないようにしています。

負の荷重は負の材料係数なのか

違います。ここで符号を持つのは $\vec B_{\rm ext}\cdot\vec n_e$ であり、固定した外部場と向き付き面法線の関係です。面積と局所質量・剛性行列は正のままです。向きを反転した面では法線traceの代表値が反転するので、局所右辺の符号も反転します。

既存のG-1-3と何が違うのか

G-1-3は平面メッシュで局所行列を全体行列へ足す一般的な入口です。G-1-B2は同じscatterを曲面starへ適用し、合同面なら全体行列まで不変になる一方、法線依存右辺は面ごとに符号を持って加算・相殺されることへ焦点を絞っています。

辺要素でも同じ符号規則なのか

考え方は同じですが、自由度の向きは面法線ではなく辺の向きと共有traceで管理します。HCurlでは局所辺向きと全体辺向きの符号写像を明示してscatterし、HDivでは面向きの符号写像を使います。このP1スカラー例は、まずlocal-to-global加算と法線荷重の符号を分けて理解する入口です。

橋渡し可視化 G-1-B3 零モードを持つ自然境界系はいつ解けるか

G-1-B2で得た五自由度の剛性行列には、定数ベクトル $\vec1$ が零モードとして残ります。対称行列なので、元の自然境界方程式

\boldsymbol K\vec a=\vec f, \qquad \boldsymbol K\vec1=\vec0

が解を持つための条件は、右辺が零モードに直交すること、すなわち

\boxed{\vec1^{\mathsf T}\vec f=0}

です。これは離散Fredholm可解性条件であり、この例では向き付き総法線流束が0であることと同じです。$\beta=0^\circ$ では正負の面寄与が相殺して条件を満たしますが、既定の $d=55^\circ,\beta=20^\circ$ では

\vec1^{\mathsf T}\vec f=0.449816387475\neq0

なので、元の自然境界問題には厳密解がありません。行列が特異だから適当に一つの節点を0へ置けばよい、という話ではありません。まず右辺の可解性を判定し、その後で零モードのゲージを固定します。

可視化では三つの扱いを比較します。

\begin{array}{ll} \text{自然境界の診断:}& \begin{bmatrix}\boldsymbol K&\vec1\\\vec1^{\mathsf T}&0\end{bmatrix} \begin{bmatrix}\vec a\\\lambda\end{bmatrix} =\begin{bmatrix}\vec f\\0\end{bmatrix},\\[8pt] \text{面積零平均:}& \begin{bmatrix}\boldsymbol K&\vec c\\\vec c^{\mathsf T}&0\end{bmatrix} \begin{bmatrix}\vec a\\\lambda\end{bmatrix} =\begin{bmatrix}\vec f\\0\end{bmatrix}, \qquad \vec c=\boldsymbol M\vec1,\\[8pt] P_0\text{固定:}& a_0=0,\qquad \boldsymbol K_{ff}\vec a_f=\vec f_f. \end{array}

第一式は診断用です。非適合右辺では $\lambda=(\vec1^{\mathsf T}\vec f)/5$ となり、$\boldsymbol K\vec a-\vec f=-\lambda\vec1$ が残るため、元の方程式の解とは呼べません。第二式は一様な分布反力密度 $\lambda$ を面積重み $\vec c$ に載せ、$\vec c^{\mathsf T}\vec a=0$ を課した補完問題です。第三式は固定節点 $P_0$ に集中反力を発生させる別の補完問題です。いずれも反力 $\vec r=\boldsymbol K\vec a-\vec f$ を含めれば

\vec1^{\mathsf T}(\vec f+\vec r)=0

となります。

ヒンジ変数 $d$ 外部場の傾き $\beta$ 可解性・拘束の扱い任意観察：$\beta=0^\circ$で二つの拘束解を比べると

$d=55^\circ,\beta=0^\circ$ の適合状態では、面積零平均解と $P_0$ 固定解は定数だけ異なり、各要素上の勾配と

E_h=\frac12\vec a^{\mathsf T}\boldsymbol K\vec a =0.0221362503083

は同じです。一方、$\beta=20^\circ$ の非適合状態で分布反力と一点反力を選ぶことは、異なる物理的補完です。両解は定数差だけでは結ばれず、勾配エネルギーもそれぞれ $0.0195468009242$ と $0.0356003573819$ になります。星形面の節点色は表示中の解から面積平均を除いた値で塗り、ゲージに依存しない勾配の比較を優先します。

零平均拘束を入れれば、非適合な元問題が解けたことになるのか

なりません。非適合時の鞍点系は、面積重みに沿う分布反力を加えた補完問題を解いています。表示する $\lambda$ と反力ベクトルを省くと、元の自然境界式を解いたように見えるので注意してください。

$P_0=0$ はいつ単なるゲージ固定なのか

$\vec1^{\mathsf T}\vec f=0$ を満たすときです。この場合、零平均解と $P_0$ 固定解は定数差だけで、勾配とエネルギーは一致します。非適合時には $P_0$ の反力が荷重総和を受け持つため、分布反力を選んだ問題とは別物です。

電磁気FEMでは何に対応するのか

純Neumann Poisson問題の電荷・流束収支、磁気スカラーポテンシャルの定数ゲージ、ベクトルポテンシャル系のより大きな核とゲージ条件に対応します。まず核と右辺の直交性を調べ、次に物理的境界条件、ゲージ、Lagrange乗数、nullspace-aware solverのどれで閉じるかを選びます。

可視化 G-1-4 四面体では節点と辺のどちらへ未知量を置くか

スカラーポテンシャルには節点要素が自然です。一方、磁気ベクトルポテンシャル $\vec A$ や電界 $\vec E$ のように接線連続性を持つ1-formは、辺に沿う積分値を自由度とする辺要素で表します。四面体の重心座標を $\lambda_i$ とすると、辺 $i\to j$ の一次Nédélec基底は

\vec W_{ij}=\lambda_i\nabla\lambda_j-\lambda_j\nabla\lambda_i

です。この基底は選んだ辺に沿う線積分だけを1にします。

\int_{e_{kl}}\vec W_{ij}\cdot d\vec\ell =\delta_{ik}\delta_{jl}-\delta_{il}\delta_{jk}

表示を節点要素と辺要素で切り替え、辺要素では選択した辺を変えてください。矢印は選択した辺基底 $\vec W_{ij}$ の要素内部での向きです。

要素の種類選択する辺

なぜ電磁界解析では辺要素がよいのか

静磁界の磁気ベクトルポテンシャル法では、$\vec B=\nabla\times\vec A$ を用いて、代表的には次の弱形式を解きます。

\int_{\Omega}\nu (\nabla\times\vec A)\cdot(\nabla\times\vec v)\,dV =\int_{\Omega}\vec J\cdot\vec v\,dV, \qquad \vec A,\vec v\in H(\operatorname{curl};\Omega)

この問題が要求するのは、要素境界をまたぐ接線成分の整合です。三成分をすべて節点値として連続にするベクトル節点要素は、必要以上に強い連続性を課します。Nédélec辺要素は接線成分を連続に保ちながら、法線成分には材料境界やゲージに応じた不連続を許せるため、curlを含むMaxwell問題の関数空間に適合します。

辺要素の未知量は、点でのベクトル値ではなく、向きを付けた辺に沿う循環です。

a_e=\int_e\vec A\cdot d\vec\ell, \qquad \vec A_h=\sum_{e}a_e\vec W_e

したがって、隣り合う四面体は共有辺の自由度を同じ一個の未知量として使えます。辺の向きを逆にすれば $a_e$ の符号も反転するため、Stokesの定理で面を回る循環を足し合わせる操作が、行列の符号としてそのまま組み立てられます。

\int_f(\nabla\times\vec A_h)\cdot\vec n\,dS =\oint_{\partial f}\vec A_h\cdot d\vec\ell =\sum_{e\subset\partial f}C_{fe}a_e

さらに、節点から辺、辺から面へ進む離散微分は、連続場の恒等式を保ちます。

H^1\xrightarrow{\nabla}H(\operatorname{curl}) \xrightarrow{\nabla\times}H(\operatorname{div}), \qquad \boldsymbol C\boldsymbol G=\boldsymbol0, \qquad \boldsymbol D\boldsymbol C=\boldsymbol0

このため、勾配場をcurlした成分や、curlを発散した成分が離散化だけを原因として現れにくくなります。特にMaxwell固有値問題では、適合するメッシュ・空間・組立てと組み合わせることで、物理的でない偽モードを抑える重要な条件になります。

| 比べる点 | ベクトル節点要素 | Nédélec辺要素 |

|---|---|---|

| 自由度 | 節点での各成分 | 向き付き辺の線積分 |

| 要素間で連続にする量 | ベクトルの全成分 | 接線成分 |

| 適合する代表空間 | $(H^1)^3$ | $H(\operatorname{curl})$ |

| 得意な未知量 | 変位など全成分を連続にしたいベクトル | $\vec E$、$\vec A$ などcurlを使う場 |

| Maxwell構造 | 後から整合を工夫する必要がある | Stokes則と離散複体を直接組み込める |

ただし、辺要素を選ぶだけで解析が自動的に正しくなるわけではありません。磁気ベクトルポテンシャルには $\vec A$ と $\vec A+\nabla\phi$ が同じ $\vec B$ を与えるゲージ自由度が残るため、ゲージ条件、木・余木分解、あるいは適切な正則化が必要です。穴のある領域では切断面やcohomology自由度も別に扱います。逆に、未知量がスカラーポテンシャルのような0-formなら節点要素、磁束のような2-formなら面要素が自然です。大切なのは「辺要素が常に優秀」ではなく、Maxwell方程式が積分する幾何学的対象と自由度の置き場所を一致させることです。

橋渡し可視化 G-1-B4 HCurlの勾配核を反復法はどう扱うか

辺要素の自由度を置いただけでは、curl–curl行列は正定値になりません。最小の例として、一個の向き付き四面体を考えます。四節点から六辺への接続行列を $\boldsymbol G$、六辺から四面への接続行列を $\boldsymbol C$ とすると、

\boldsymbol G= \begin{bmatrix} -1&1&0&0\\ -1&0&1&0\\ -1&0&0&1\\ 0&-1&1&0\\ 0&-1&0&1\\ 0&0&-1&1 \end{bmatrix}, \qquad \boldsymbol C= \begin{bmatrix} 0&0&0&1&-1&1\\ 0&-1&1&0&0&-1\\ 1&0&-1&0&1&0\\ -1&1&0&-1&0&0 \end{bmatrix}

です。境界の境界は零なので $\boldsymbol C\boldsymbol G=\boldsymbol0$ が数値近似ではなく整数行列の恒等式として成立します。面に置いた材料Hodge行列 $\boldsymbol H_2$ を使う最小curl–curlエネルギーは

E_h(\boldsymbol a)=\frac12 (\boldsymbol C\boldsymbol a)^T\boldsymbol H_2 (\boldsymbol C\boldsymbol a), \qquad \boldsymbol K=\boldsymbol C^T\boldsymbol H_2\boldsymbol C

です。したがって任意の節点ポテンシャル $\boldsymbol\phi$ に対し

\boldsymbol K\boldsymbol G\boldsymbol\phi=\boldsymbol0

となります。ここでは定数節点モードを除いた $\boldsymbol G_r\in\mathbb R^{6\times3}$ が三次元の勾配核を張り、$\operatorname{rank}\boldsymbol K=3$、$\operatorname{nullity}\boldsymbol K=3$ です。

下の可視化では、右辺

\boldsymbol f(\eta) =\boldsymbol C^T\boldsymbol s +\eta\boldsymbol G\boldsymbol q

へ勾配核成分を混ぜます。$\eta=0$ なら $\boldsymbol G_r^T\boldsymbol f=\boldsymbol0$ で元式は可解、$\eta>0$ なら非適合です。反復回数を動かし、次の三つを比べてください。

nullspace射影CG：$\boldsymbol P=\boldsymbol I-\boldsymbol G_r(\boldsymbol G_r^T\boldsymbol G_r)^{-1}\boldsymbol G_r^T$ を使い、$\boldsymbol K\boldsymbol a=\boldsymbol P\boldsymbol f$ を解きます。射影残差は零へ行きますが、$\eta>0$ では元式の残差 $\boldsymbol K\boldsymbol a-\boldsymbol f$ が残ります。
tree–cotree CG：全節点へつながるtree辺 $e_0,e_1,e_2$ を零に固定し、cotree辺 $e_3,e_4,e_5$ の正定値縮約系を解きます。$\eta>0$ ではtree辺に拘束反力を許した別問題です。
未処理CG：特異行列へ非適合右辺をそのまま渡します。勾配核方向の残差を消せず、破綻します。

非適合成分 η CG反復回数核の扱い

η=0、反復3回で射影CGとtree–cotree CGを比べると何が起きますか。

可解性、ゲージ、反力を一つの残差と呼ばない

$\eta=0$ では射影解とtree固定解の差は $\operatorname{range}\boldsymbol G_r$ にあり、$\boldsymbol C\boldsymbol G_r=0$ なので面curlとエネルギーは一致します。これは「二つのソルバーが同じ辺係数を返す」という意味ではありません。ゲージ条件が違えば $\boldsymbol a$ は違ってよく、物理量 $\boldsymbol C\boldsymbol a$ を比較します。

$\eta>0$ では話が変わります。射影法は右辺から非適合成分を取り除いた系を解きます。一方tree–cotree法は、固定したtree辺に反力が生じる拘束問題を解きます。どちらも計算は止められますが、元の $\boldsymbol K\boldsymbol a=\boldsymbol f$ を解いたとは主張できません。可解性指標 $\|\boldsymbol G_r^T\boldsymbol f\|$、元式残差 $\|\boldsymbol K\boldsymbol a-\boldsymbol f\|$、射影残差 $\|\boldsymbol P(\boldsymbol K\boldsymbol a-\boldsymbol f)\|$、拘束辺反力を分けて記録する必要があります。

この小模型と実務FEMの境界

ここで手計算しているのは、向き付き四面体の正確な接続行列と離散HodgeによるWhitney 1-formの最小模型です。高次HCurl要素、Piola写像、曲面幾何、要素向き、局所・大域自由度変換、積分と弱形式組立てはNGSolveなどの有限要素基盤へ任せます。実務でこの六辺行列をPythonから再構成する教材ではありません。

前処理では何を守る必要があるのか

前処理後のKrylov空間でも既知の勾配核を混ぜ戻さないことが重要です。射影法なら前処理の前後で核直交性を維持し、AMSなどの補助空間法なら勾配空間を前処理の構造として明示します。この可視化は前処理器そのものではなく、前処理付き反復法が報告すべき可解性と残差の基準を与えます。

穴のある領域でもtreeだけで十分か

十分ではありません。treeは局所的な勾配ゲージを固定しますが、多重連結領域にはcurlも発散も零で勾配ではない大域調和自由度が残り得ます。G-6-3の穴あきセル複体でtree後の基本閉路を面境界で商にし、G-6-4でcohomology周期とHodge最小代表を分けてください。

G-1-B4からG-6へ：単連結の核と穴の自由度を混ぜない

この四面体では $\ker\boldsymbol C=\operatorname{range}\boldsymbol G$ なので、curl-freeな辺場はすべて節点ポテンシャルの勾配です。穴のある複体では次の商空間が非零になり、tree固定後にも大域自由度が残ります。

H^1_h =\dfrac{\ker\boldsymbol C}{\operatorname{range}\boldsymbol G}, \qquad \dim H^1_h=b_1

ここで検査した量	穴あき複体で追加する検査	次に動かす図
$\boldsymbol C\boldsymbol G=0$、$\operatorname{nullity}\boldsymbol K=3$	$\ker\boldsymbol C$ のうち $\operatorname{range}\boldsymbol G$ でないclass	G-6-3 tree–cotreeと商空間
tree辺を零にするゲージ固定	cotree基本閉路から面境界を商にし、$b_1$個だけ残す	G-6-3 穴を回る $z_5$
curlエネルギーでゲージ違いを比較	周期を固定し、材料Hodgeで最小エネルギー代表を選ぶ	G-6-4 離散調和代表
非適合右辺を射影・反力で区別	独立循環を切断面またはcohomology基底の係数として与える	G-6-2 切断面

したがって、G-1-B4の「nullspaceを消す」という言い方を、そのまま穴あき領域へ延長してはいけません。消すのは局所gradゲージであり、物理的に指定すべき独立循環まで零にしないことが次の検査です。

実務ゲート G-1-B5 多要素HCurlソルバーは何を保存するか

一個の四面体で見た量を、多要素・複数連結成分・穴のある解析へ移すときは、ソルバーの「収束」だけを合否にしません。grad核の基底を $\boldsymbol Z_g$、穴に由来する調和基底を $\boldsymbol Z_h$、両者を並べた核基底を $\boldsymbol Z=[\boldsymbol Z_g\ \boldsymbol Z_h]$ とします。純粋なcurl–curl系なら

\boldsymbol K\boldsymbol Z=\boldsymbol0, \qquad \boldsymbol Z^{\mathsf T}\boldsymbol f=\boldsymbol0

が、拘束を加える前の可解性条件です。ただし、$\boldsymbol Z_g$ は同じ物理場を表す局所ゲージなので固定できますが、$\boldsymbol Z_h$ の係数は独立循環という物理条件になり得ます。調和成分を「nullspaceだから」と自動射影してはいけません。

段階	必ず保存する診断値	合格の読み方	失敗時に疑うもの
メッシュと離散複体	連結成分数 $b_0$、free節点・辺・面数、$\lVert\boldsymbol C\boldsymbol G\rVert$	essential境界条件を反映した自由度上で $CG=0$	辺・面向き、共有自由度、境界自由度の除去
局所grad核	$\operatorname{rank}\boldsymbol Z_g$、$\lVert\boldsymbol K\boldsymbol Z_g\rVert$	各連結成分の定数を除いたgrad自由度を再現	成分ごとのtree、離散gradient、拘束rank
大域位相	$b_1$、$\operatorname{rank}\boldsymbol Z_h$、周期行列 $\boldsymbol L^{\mathsf T}\boldsymbol Z_h$	独立閉路と基底の周期行列がfull rank	face boundaryとの商、切断面、cohomology基底
右辺	$\lVert\boldsymbol Z_g^{\mathsf T}\boldsymbol f\rVert$ と $\lVert\boldsymbol Z_h^{\mathsf T}\boldsymbol f\rVert$ を別記	grad非適合と大域循環条件の不足を分離	荷重組立て、境界条件、指定電流・循環
ゲージ・周期拘束	gauge残差、$\lVert\boldsymbol L^{\mathsf T}\boldsymbol a-\boldsymbol I\rVert$	gradだけを固定し、指定周期を保持	treeで大域cycleまで消した、周期符号
Krylov反復	元式残差、射影残差、自由DOF残差、反力、反復ごとの核直交	何を射影・拘束した系へ収束したか説明できる	非適合右辺、前処理後の核混入、誤った停止判定
物理量	$\boldsymbol C\boldsymbol a$、curlエネルギー、独立周期	ゲージを変えてもcurlと指定周期が一致	gauge依存量を物理解として比較した

周期を物理条件として与える場合は、調和基底を消す代わりに、例えば

\boldsymbol L^{\mathsf T}\boldsymbol a=\boldsymbol I

を追加した拘束系、またはcohomology基底係数を未知量へ加えた定式化を使います。ここで $\boldsymbol L$ は独立閉路の周期を読む作用素です。$\boldsymbol Z_h$ を含む全核を射影してから「収束した」と報告すると、指定すべき循環まで捨てる危険があります。

実装順序：rankを数えてから前処理を選ぶ

essential境界条件を反映した自由HCurl空間を確定し、連結成分ごとに接続とrankを検査する。
grad核 $\boldsymbol Z_g$ とcohomology基底 $\boldsymbol Z_h$ を別々に作り、$\boldsymbol C\boldsymbol Z_g=0$、$\boldsymbol C\boldsymbol Z_h=0$、周期行列のrankを確認する。
右辺を二つの核へ投影して可解性を診断し、ゲージ固定と物理的周期条件を選ぶ。
射影を使うなら前処理の前後で同じ核直交を保つ。補助空間法を使うならgrad空間を前処理構造へ明示する。
元式・射影系・拘束系のどれを解いたかと、反力・周期残差・curl不変量を同じrun記録へ残す。

有限要素基盤との境界：高次HCurlの要素向き、局所・大域DOF変換、Piola写像、曲面幾何、数値積分、弱形式組立てはNGSolveの空間・form・GridFunctionへ任せます。この表は手製基底実装ではなく、NGSolveが作った離散系を受け入れる診断契約です。

$\operatorname{rank}G=N_v-b_0$ をいつ使えるか

essential境界条件をまだ除いていない頂点–辺接続で使えるグラフrankです。自由DOFへ縮約した後、周期・対称・同一視境界がある場合は、元の節点数だけから決めず、実際の縮約作用素についてrankを取り直します。

質量項やゲージ罰則があれば核検査は不要か

不要にはなりません。質量項や罰則で行列が正則になっても、元のde Rham構造、低周波極限、条件数、物理的循環は残ります。正則化がgrad成分や大域周期をどれだけ変えたかを、$C a$、周期、エネルギーで別に検査します。

実測可視化 G-1-B6 NGSolve T–Ω解析で七段階ゲートを通す

G-1-B5の受入契約を、実際の多要素解析へ適用します。題材は、$z$方向へ流れる直線電流を囲む環状空気領域です。穴を一つ持つため $b_1=1$ であり、磁界をスカラーポテンシャルの勾配だけでは表せません。Radiaのcohomology engineがNGSolveメッシュから単位周期のHCurl 1-form $\boldsymbol h$ を作り、NGSolveのH1弱形式で

\vec H_h =-\nabla\phi_h+NI\,\boldsymbol h, \qquad \nabla\times\boldsymbol h=\vec0, \qquad \oint_{\gamma_{\rm hole}}\boldsymbol h\cdot d\vec\ell=1

を解きます。ここで $NI$ は電流ターン数、$\gamma_{\rm hole}$ は穴を一周する閉路です。比較する厳密解と物理的周期は

H_\varphi(r)=\dfrac{NI}{2\pi r}, \qquad \oint_{\gamma_r}\vec H\cdot d\vec\ell=NI

です。実行済みの正典はRadiaの result-bearing notebook「docs/cohomology/tomega_wire.ipynb」です。ブラウザはNGSolve 6.2.2606で保存した三メッシュの結果を再生し、有限要素基底や組立てをJavaScriptで再実装しません。電流スライダーだけは線形問題の比例則で保存値を拡大縮小します。

実測メッシュ電流ターン数 $NI$ 保存済み評価半径位相成分

この実測で、位相の正しさと点値精度をどう判定しますか。

実測値：位相量と点値誤差は同じ収束指標ではない

三つの実行はすべて一連結成分・一つの穴を持ち、$\boldsymbol C\boldsymbol G=\boldsymbol0$、$b_1=1$、$\operatorname{rank}(\boldsymbol L^{\mathsf T}\boldsymbol Z_h)=1$ を保ちました。選択した離散基本閉路上ではperiod matrixを単位行列へ正規化しています。下表の「穴周期」は、別に置いた幾何学的な円周上でGridFunctionを数値積分した値なので、補間・周回積分誤差を含みます。

maxh	四面体	HCurl(0) DOF	H1(2) DOF / free	$\lVert\operatorname{curl}h\rVert/\lVert h\rVert$	円周穴周期	Ampère相対誤差	$H_\varphi$平均点値誤差
12 mm	10,791	13,953	16,226 / 12,670	$8.56\times10^{-13}$	0.998430	$1.81\times10^{-4}$	2.002%
8 mm	14,370	18,807	21,916 / 16,604	$1.05\times10^{-12}$	0.999174	$2.84\times10^{-5}$	1.164%
6 mm	22,745	29,629	34,516 / 26,528	$1.06\times10^{-12}$	1.001830	$1.94\times10^{-4}$	0.750%

点値の平均誤差は2.002%から0.750%へ低下しました。一方、別円周上のAmpère誤差は単調ではありません。これは位相rankが壊れたのではなく、HCurl GridFunctionを幾何学的円へ補間して積分する観測誤差が混ざるためです。合格条件は「すべての数値が単調」ではありません。離散複体とperiod matrixのrankを厳密に検査し、円周積分と点値には別の許容差・細分収束を設けます。

G-1-B5の七段階をこのrunへ対応させる

メッシュ：$b_0=1$、頂点・辺・面・四面体数と $\lVert CG\rVert_\infty=0$ を保存する。
局所grad：$\operatorname{rank}G=N_v-1$ を確認し、H1の外側Dirichlet境界でスカラーポテンシャルの定数ゲージを固定する。
大域位相：$b_1=1$、cohomology基底数1、period matrix rank 1を確認する。
右辺：指定値 $NI$ とcohomology係数の個数を一致させる。勾配だけの反例は線形方程式を解けてもAmpère周期0の別問題である。
周期拘束：離散基本閉路では単位周期、別の穴円周では約1、可縮円周では $10^{-5}$ 程度であることを分けて保存する。
線形解法：この規模はNGSolveの直接逆行列を使用した。CGやAMSを使っていないrunへ反復収束を捏造しない。
物理量：Ampère周回積分、半径8点の $H_\varphi$、平均・最大点値誤差を保存し、位相検査と精度検査を閉じる。

実装所有境界：NGSolveがH1空間、自由DOF、弱形式、積分、直接解法、GridFunction評価を所有します。Radiaは最低次離散複体からcohomology classと単位周期を作り、T–Ωの物理的な大域自由度を与えます。高次HCurlの要素向きやPiola写像をブラウザやPython教材で再構成していません。

実測可視化 G-1-B7 場の誤差と観測演算子の誤差を分ける

G-1-B6では、離散period matrixが正しくても、幾何学的な円周積分と一点の $H_\varphi$ は別の誤差を持つことが分かりました。ここでは同じNGSolve T–Ω解を固定したまま、どこで、何点で、どの位相から測るかだけを変えます。厳密場を $\vec H$、有限要素場を $\vec H_h$、連続な観測量を $Q$、$M$点・開始位相 $\theta_0$ の観測演算子を $\mathcal O_{M,\theta_0}$ とすると、

Q(\vec H)-\mathcal O_{M,\theta_0}(\vec H_h) =\underbrace{Q(\vec H)-Q(\vec H_h)}_{\text{場の離散誤差}} +\underbrace{Q(\vec H_h)-\mathcal O_{M,\theta_0}(\vec H_h)}_{\text{観測演算子誤差}}

です。穴の離散基本閉路では $\boldsymbol L^{\mathsf T}\boldsymbol Z_h=\boldsymbol I$ を代数的に検査できます。一方、メッシュ辺と一致しない幾何学的な円 $\gamma_r$ では、

\mathcal O_{M,\theta_0}(\vec H_h) =\sum_{j=0}^{M-1}\vec H_h\!\left(\vec x(\theta_j)\right) \cdot\vec t(\theta_j)\,r\,\Delta\theta, \qquad \theta_j=\theta_0+j\Delta\theta

と評価するため、要素探索、補間、要素境界を横切る位置、曲線サンプリングが入ります。これはperiod rankの失敗ではありません。

有限要素場観測実験点評価半径円周サンプル数開始位相積分円半径

どの操作が有限要素場そのものを変えますか。

点値は境界へ近づくほど厳しく、局所的には単調収束しない

点評価は連続な $H(\operatorname{curl})$ エネルギー評価より強い要求です。$z=0$、$x=r$ 上で保存した実測の一部を示します。

$r$	近い境界までの距離	標準メッシュ誤差	細メッシュ誤差
5.5 mm	内壁0.5 mm	4.014%	2.905%
27.5 mm	両壁22.5 mm	0.381%	0.271%
40.0 mm	外壁10.0 mm	4.684%	3.034%
45.0 mm	外壁5.0 mm	11.753%	7.128%
49.5 mm	外壁0.5 mm	15.057%	13.876%

14点平均は4.890%から4.011%、最大は15.057%から13.876%へ下がりました。しかし $r=10,15,20,35$ mmでは細メッシュの一点誤差が標準より一時的に大きくなります。したがって「全体として改善」と「すべての点で単調」は別の主張です。

円周点数を増やしても、途中の値は単調とは限らない

半径27.5 mmの同じ円で、標準メッシュ・開始位相0のAmpère相対誤差は、24点で$3.19\times10^{-3}$、192点で$2.54\times10^{-5}$、384点で$1.73\times10^{-5}$、768点で$2.67\times10^{-5}$、6,144点で$2.96\times10^{-5}$でした。192点が後続すべてより正確に見えるのは偶然の相殺です。半刻みずらすと同じ$M$でも値が変わります。

細メッシュでは6,144点で二つの開始位相の相対誤差が$1.95\times10^{-6}$と$1.16\times10^{-6}$、周回積分差は$7.94\times10^{-5}$ Aまで縮みました。点数系列、開始位相系列、メッシュ系列を一つの収束線へ混ぜず、各軸を独立に振ります。

G-1-B7の受入順序

位相：まず離散基本閉路で$b_1$と$\boldsymbol L^{\mathsf T}\boldsymbol Z_h=\boldsymbol I$を検査する。
場：メッシュまたは次数を変え、エネルギー量・平均量・必要な局所量の収束を調べる。
観測：同じ$\vec H_h$を固定し、点位置、境界距離、円半径、分割数、開始位相だけを変える。
誤差床：観測を細かくしても残る差を場の離散誤差候補とし、別メッシュで再確認する。
報告：「period rank」「場の離散系列」「観測演算子系列」を別欄に保存する。

実装所有境界：点の要素探索とGridFunction評価はNGSolveへ任せます。幾何円のサンプリングは検証用の後処理であり、高次HCurl basisや要素向きの再実装ではありません。実務の線・面観測には、可能ならNGSolveの積分・trace・組立て作用素を優先します。

実測可視化 G-1-B8 一点値を有限幅probe観測へ置き換える

G-1-B7では、点値が境界距離と一個の要素に強く依存することを確認しました。そこで観測点 $\vec x_0=(r_0,0,0)$ の周囲へGaussian感度

w_{\sigma,\vec x_0}(\vec x) =\exp\!\left(-\dfrac{|\vec x-\vec x_0|^2}{2\sigma^2}\right)

を置き、固定した物理方向 $\vec e_y$ の磁界を

J_{\sigma,\vec x_0}(\vec H) =\dfrac{\displaystyle\int_\Omega w_{\sigma,\vec x_0}\,\vec H\cdot\vec e_y\,dV} {\displaystyle\int_\Omega w_{\sigma,\vec x_0}\,dV}

で読みます。これは有限体積のセンサ応答です。$\sigma$を固定すれば、Cauchy--Schwarz不等式から

|J_{\sigma}(\vec H)-J_{\sigma}(\vec H_h)| \leq \dfrac{\lVert w_\sigma\rVert_{L^2(\Omega)}} {\displaystyle\int_\Omega w_\sigma\,dV} \lVert H_y-H_{h,y}\rVert_{L^2(\Omega)} =V_{\mathrm{eff}}^{-1/2}\lVert H_y-H_{h,y}\rVert_{L^2(\Omega)}, \qquad V_{\mathrm{eff}} =\dfrac{\left(\int_\Omega w_\sigma\,dV\right)^2} {\int_\Omega w_\sigma^2\,dV}

となり、$L^2$場誤差で制御できる観測汎関数になります。ただし、$\sigma\to0$では$V_{\mathrm{eff}}\to0$となり、安定性係数は増大して一点評価へ戻ります。「平均すれば必ず正確」なのではなく、何を測るかを有限幅で定義すると誤差をノルムで制御できるという意味です。

有限要素場 probe中心観測定義 Gaussian幅 $\sigma$

有限要素場の細分効果だけを比較する操作はどれですか。

誤差を三つへ分ける

ブラウザの「厳密probe」は、解析解 $H_y=Ix/[2\pi(x^2+y^2)]$ を円筒座標へ移し、scaled modified Bessel関数を使った独立一次元積分で求めます。NGSolve側は同じGaussianを12次積分し、有限要素場と解析場を読みます。中心点を最終目標にした差は、符号付きで

H_y(\vec x_0)-J^{q}_{\sigma,h}(\vec H_h) =\underbrace{H_y(\vec x_0)-J_{\sigma}(\vec H)}_{\text{センサ幅による観測バイアス}} +\underbrace{J_{\sigma}(\vec H)-J^{q}_{\sigma,h}(\vec H)}_{\text{幾何・積分誤差}} +\underbrace{J^{q}_{\sigma,h}(\vec H)-J^{q}_{\sigma,h}(\vec H_h)}_{\text{有限要素場誤差}}

と分かれます。絶対値を別々に報告し、符号の相殺で総差だけが小さく見える状態を合格にしません。

実測で分かったこと

中心と観測	粗	標準	細	別に残る量
中央27.5 mm・一点	1.409%	0.381%	0.271%	一点評価
中央27.5 mm・$\sigma=3$ mm	1.733%	0.286%	0.134%	probe bias $4.3\times10^{-12}$%
外壁0.5 mm・一点	0.867%	15.057%	13.876%	粗の一点は偶然よく見える
外壁0.5 mm・$\sigma=12$ mm	0.504%	3.108%	1.292%	中心点に対するprobe bias 20.622%
中央27.5 mm・$\sigma=0.5$ mmの積分誤差	0.568%	0.941%	0.022%	$\sigma\ll h$では観測積分も未解像

中央で境界から十分離れた$\sigma\leq4$ mmでは、解析解のGaussian平均は中心点値と機械精度で一致しました。これは源のない近傍で$H_y$が調和関数となる平均値性の現れです。一方、内壁・外壁でGaussianが$\Omega$に切られると対称性がなくなり、有限要素誤差が下がっても中心点を表すセンサではなくなります。外壁近傍の$\sigma=12$ mmでは細メッシュのprobe汎関数誤差は1.292%まで下がりますが、中心点へのバイアスは20.622%あります。

G-1-B8の受入順序

観測定義：中心、物理成分、感度関数、$\sigma$、領域での切れ方を先に固定する。
独立真値：解析解または別経路で$J_\sigma(\vec H)$を求め、中心点値と混同しない。
積分：解析場を同じメッシュ・数値積分で読み、$\sigma/h$不足による観測積分誤差を測る。
場：同じ観測汎関数を固定したままメッシュまたは次数だけを変える。
センサ設計：場誤差と積分誤差が収束した後で、$\sigma$と空間分解能・中心点バイアスの交換条件を選ぶ。
報告：一点誤差、probe汎関数誤差、積分誤差、probe bias、総差を別欄へ保存する。

実装所有境界：NGSolveが有限要素場、曲面メッシュ、CoefficientFunction、12次数値積分を所有します。独立真値は解析場の円筒対称性とscaled Bessel積分で求めます。ブラウザは保存済み結果を再生し、要素内評価や積分をJavaScriptで模倣しません。

実測可視化 G-1-B9 随伴誤差推定と要素マーキングを分ける

G-1-B8の中央$\sigma=3$ mm Gaussian probeを、ここでは固定した目的汎関数にします。主問題を

\[ a(u,v)=F(v) \]

と書くと、磁界の$y$成分に対するpotential部分の観測は

J(v)=-\dfrac{\displaystyle\int_\Omega w_{\sigma,\vec x_0}\, \partial_y v\,dV}{\displaystyle\int_\Omega w_{\sigma,\vec x_0}\,dV}

です。同じ剛性作用素へ別の右辺を与える随伴問題

\[ a(v,z)=J(v) \]

を解けば、主問題残差$R(u_h)(v)=F(v)-a(u_h,v)$に対して

\[ J(u)-J(u_h)=R(u_h)(z-z_h) \]

となります。これがdual-weighted residual（DWR）の出発点です。しかし、誤差の符号付き総和を当てることと、各要素寄与の絶対値を大きい順に細分すれば必ず改善することは同じ定理ではありません。

表示する状態要素表示

DWR推定値のeffectivityが1に近いのに、naive DWR細分が悪化した理由はどれですか。

推定器は当たったが、naive markingは外れた

粗メッシュは10,791四面体、P1で2,273 DOFです。probe真値578.745 A/mに対し、P1は547.588 A/mで5.384%誤差、同じメッシュのP2豊富化は568.856 A/mで1.709%誤差でした。P2とP1の観測差は21.268 A/m、P2随伴から得たDWR推定値は22.176 A/mで、

\mathcal I_{\rm eff} =\dfrac{\eta_{\rm DWR}}{J(u_{P2})-J(u_{P1})} =1.043

です。これはP2豊富化補正に対するeffectivityであり、連続真値に対する完全な保証値とは呼びません。

一方、要素寄与$\eta_K=R_K(u_{P1})(z_{P2}-I_1z_{P2})$は

\sum_K\eta_K=22.176\ {\rm A/m},\qquad \sum_K|\eta_K|=591.342\ {\rm A/m}

解析	標識	最終DOF	四面体	probe誤差
P1基準	なし	2,273	10,791	5.384%
P2豊富化	同じメッシュ・次数だけ上昇	16,226	10,791	1.709%
energy余剰細分	上位1%・108要素	4,279	21,440	1.153%
naive $\|\eta_K\|$細分	上位6%・647要素	4,293	21,833	16.300%

細分後のDOF差は14、0.33%だけで、予算差では説明できません。両方とも$b_1=1$を再確認しています。energy標識とnaive DWR標識の重なりは31要素だけでした。naive DWR側はprobe近傍へ強く局在しましたが、再解析すると符号相殺と大域的な主問題誤差の釣合いが変わり、観測値は484.412 A/mへ悪化しました。随伴を否定する反例ではなく、cellwise signed residualを絶対値順位へ直結するnaive localizationへの反例です。

G-1-B9の受入順序

固定目標：中心、$\sigma$、物理成分、積分次数を変えずに$J$を定義する。
同じ作用素：主問題と随伴問題が同じ剛性行列、異なる右辺であることを保存する。
豊富化：P1主問題に対してP2主問題・P2随伴を解き、effectivityの分母を明記する。
符号：$\sum\eta_K$、$\sum|\eta_K|$、正負和、相殺比を別々に残す。
同一予算：標識率ではなく、閉包後の実DOFと要素数を比較する。
再解析：細分後にcohomology、$b_1$、period、目的汎関数を解き直し、推定値だけで合格させない。

次の弱点：安定な局在化には、face jumpを含む強残差、partition-of-unityまたはpatchwise DWR、符号を保つマーキング、複数回のsolve--estimate--mark--refine受入が必要です。次はこのnaive反例を基準に、局在化そのものを修復します。

実装所有境界：NGSolveがH1空間、cohomology基底、弱形式、P1/P2 GridFunction、同一剛性作用素、要素積分、局所細分を所有します。ブラウザは保存済みの主問題・随伴・再解析結果だけを再生します。

実測可視化 G-1-B10 PU-DWRでも自動採用せず、再解析でロールバックする

G-1-B9のcellwise DWR反例に対し、まず古典的なpartition of unity（PU）で要素寄与を節点patchへ集約します。P1節点基底を$\{\psi_i\}$、P2随伴余剰を$q=z_{P2}-I_1z_{P2}$とすると、

\sum_i\psi_i=1, \qquad \eta_i^{\rm PU}=R(u_{P1})(q\psi_i), \qquad \sum_i\eta_i^{\rm PU}=R(u_{P1})(q)

です。NGSolveでは積の勾配を展開し、

\eta_i^{\rm PU} =\int_\Omega\mu_0\vec H_{P1}\cdot \left[(\nabla z_{P2}-\nabla I_1z_{P2})\psi_i +(z_{P2}-I_1z_{P2})\nabla\psi_i\right]dV

をP1 LinearFormとして一度に組み立てました。PU総和は22.175859876 A/mでG-1-B9のDWR総和と一致します。絶対値総和はcellwiseの591.342 A/mから94.475 A/mへ下がり、相殺比は3.750%から23.473%へ改善しました。しかし、局在化が滑らかになったことと、その標識を採用してよいことはまだ別です。

再解析候補表示する証拠

PU相殺比が改善したのに、PU細分を採用しない理由はどれですか。

五候補を同じ固定probeで再解析する

PUの$|\eta_i|$上位66節点patchは、絶対値総和の82.039%を含み、閉包後4,274 DOFになりました。正寄与と負寄与を同じ捕捉率にする二側標識は、正15 patchで77.042%、負51 patchで77.774%を捕捉しました。hybridはenergy上位81要素とPU上位16 patchを合わせます。結果は次の通りです。

候補	標識	最終DOF	probe誤差	受入
P1基準	なし	2,273	5.384%	現在値
energy安全基準	108要素	4,279	1.153%	accept
cellwise $\|\eta_K\|$	647要素	4,293	16.300%	rollback
PU $\|\eta_i\|$	66 patch・797要素	4,274	15.448%	rollback
PU正負二側	15正＋51負 patch・890要素	4,568	15.445%	rollback
energy + PU hybrid	81要素＋16 patch・和集合291要素	4,199	9.615%	rollback

全候補で$b_1=1$、$CG=0$、period rank 1を再確認しました。したがって失敗は位相消失ではありません。PU恒等式はDWR総和の分割として正しく、相殺も減らしましたが、このP2豊富化差から作った一回のpatch順位は、細分後の非線形なメッシュ応答を予測できませんでした。よい局在化の必要条件を満たしても、採用の十分条件ではないという実測です。

受入判定を

{\rm accept}(m)= [b_1(m)=1]\land[\mathrm{rank}\,P(m)=1] \land[\mathrm{DOF}(m)\le B] \land[E_J(u_m)\lt E_J(u_{\rm current})]

とし、候補を別メッシュで解き直してから状態を進めます。この比較ではenergyだけをacceptし、他は現在の基準解へrollbackします。推定器の値を「次の解の保証」と読み替えないため、rollbackも適応解析アルゴリズムの一部です。

G-1-B10の受入順序

代数恒等式：$\sum_i\psi_i=1$と$\sum_i\eta_i^{\rm PU}=\eta_{\rm DWR}$を先に再計算する。
符号記録：正負和、絶対値和、捕捉率を別々に保存し、個数だけの正負均等化をしない。
候補生成：cellwise、PU、二側、hybrid、energyを別候補にし、標識集合を混同しない。
実予算：標識数ではなく、閉包後の実DOF・四面体数を記録する。
再解析：各候補でcohomology基底と主問題を作り直し、固定した$J_\sigma$を再積分する。
commit / rollback：位相、予算、目的誤差をすべて通った候補だけを次のcurrent meshへ進める。

次の弱点：PU化だけでは修復できませんでした。次はP3または独立mesh-enriched随伴、体積残差とinterior-face jumpの明示分解、equilibrated patch flux、複数cycleのtrust-region受入を比較し、失敗した候補がどの項を欠いたかを切り分けます。

実装所有境界：NGSolveがP1/P2空間、節点基底、LinearForm、積の勾配、要素向き、局所細分、cohomology、再解析を所有します。JavaScriptはPU基底やjumpを再実装せず、保存済み候補とaccept / rollbackだけを再生します。

実測可視化 G-1-B11 面残差と平衡化fluxをtrust-regionで受理する

G-1-B10では、PU-DWRの総和恒等式が成立しても、次のmeshで目的誤差が減るとは限りませんでした。そこで同じ直線電流T–Ω問題、中央27.5 mm・幅$sigma=3$ mmのGaussian probeを固定し、随伴をP2からP3へ上げます。P3随伴$z_3$とそのP1補間$I_1z_3$から$q=z_3-I_1z_3$を作ると、

\eta_{\rm DWR}=R(u_1)(q) =\int_\Omega\mu_0\vec H_1\cdot (\nabla z_3-\nabla I_1z_3)\,dV

です。基準meshでは$eta_{\rm DWR}=21.666214$ A/m、P3–P1のprobe差は21.082130 A/mで、その比は1.027705でした。P2の1.0427より1に近づきますが、これは豊富化差に対するeffectivityであって、連続真値に対する保証ではありません。

P1場を要素ごとに部分積分し、体積と面を明示的に分けます。

R(u_1)(q)=\sum_K\left[ -\int_K\nabla\!\cdot(\mu_0\vec H_1)q\,dV +\int_{\partial K}\mu_0\vec H_1\!\cdot\vec n_K\,q\,dS \right].

この最低次・直線四面体の実測では、強体積残差の絶対値和は$1.50\times10^{-11}$ A/mにすぎず、向き付き要素面和21.666214 A/mがdirect residualへ閉じました。NGSolveのskeleton積分とOther()で作るinterior-face jumpは21.642107 A/m、trace remainderは0.024107 A/mです。したがって「面jumpが支配的」は言えますが、traceを捨てて完全一致とはしません。

もう一つの候補は、NGSolveの$H(\mathrm{div})$法線連続性を使う平衡化fluxです。$\boldsymbol\sigma_h\in H(\mathrm{div})$、$p_h\in L^2$と平均値拘束$c$について、

(\boldsymbol\sigma_h,\boldsymbol\tau) +(\nabla\!\cdot\boldsymbol\tau,p_h) +(\nabla\!\cdot\boldsymbol\sigma_h,q) +c(q,1)+d(p_h,1) =(\mu_0\vec H_1,\boldsymbol\tau), \qquad \xi_K=\int_K\frac{|\mu_0\vec H_1-\boldsymbol\sigma_h|^2}{\mu_0}\,dV

を混合UMFPACK solveで解きます。基準meshの$\|\nabla\!\cdot\boldsymbol\sigma_h\|_{L^2}=1.88\times10^{-18}$であり、$\xi_K$は法線不連続を継ぎ合わせた候補生成器になります。ただし、主・随伴の保証付き目的量区間をまだ作っていないため、ここでは目的誤差の上下界とは呼びません。

trust-region cycle 表示する証拠

最終P1誤差0.068%だけを、未知解の問題で停止証明にできない理由はどれですか。

半径を縮め、改善しない候補は必ず戻す

cycle 1では基準meshの$\xi_K$上位108要素を細分し、P1 probe誤差を5.384%から0.876%へ下げてacceptしました。cycle 2では現在meshの上位1%（207要素）をenergy、jump、fluxで別々に試しましたが、すべて悪化しました。そこで半径を104要素へ縮め、flux候補だけが0.068%へ改善したためacceptします。cycle 3で半径を再び1%（337要素）へ広げると、三候補とも現在値より悪化したのでrollbackし、0.1%許容値を満たす現在meshを保持して停止します。

cycle・候補	標識半径	最終DOF	probe誤差	判定
基準P1	なし	2,273	5.384%	current
cycle 1・equilibrated flux	108要素	4,160	0.876%	accept
cycle 2 wide・energy / jump / flux	各207要素	8,291 / 5,141 / 7,928	1.190% / 9.842% / 1.292%	全rollback
cycle 2 shrink・energy / flux	各104要素	7,038 / 6,659	1.954% / 0.068%	energy rollback / flux accept
cycle 3 re-expand・energy / jump / flux	各337要素	10,546 / 8,109 / 10,168	0.209% / 10.009% / 1.933%	全rollback → stop

全11 solveで$b_1=1$、$CG=0$、period rank 1、実DOF 12,000以下を再確認しました。したがってrollbackは位相消失や予算超過ではなく、固定した目的量の悪化によるものです。

P3が教える「近すぎるP1」への注意

current mesh	P1 probe誤差	P3 probe誤差	P3-DWR / (P3−P1)
基準	5.384%	1.741%	1.0277
cycle 1 accept後	0.876%	0.614%	0.9885
cycle 2 shrink accept後	0.068%	0.991%	0.9871

最終meshではP1が解析真値578.745 A/mに偶然近い一方、P3は584.479 A/mへ行き過ぎます。さらにcycle 3のenergy候補もP1誤差0.209%へ悪化しました。したがって、単調なh細分やp上昇を仮定して「0.068%だから収束した」とは言えません。この教材でsafe stopできるのは、解析解を持つbenchmarkとして再解析後の真値誤差を測れたからです。

G-1-B11の受入順序

P3診断：P3–P1 correctionとDWRのeffectivityを測るが、連続真値の保証とは呼ばない。
残差分解：強体積、向き付き要素面、interior jump、trace remainderを別々に保存する。
平衡化：$H(\mathrm{div})\times L^2\times\mathbb R$混合問題を解き、法線連続性とdivergenceを検査する。
trust半径：wide候補が全滅したら半径を縮め、改善後は再拡大して局所最適らしさを試す。
再解析：各候補でcohomologyと主問題を作り直し、同じGaussian probe、位相、実DOFで判定する。
commit / rollback：改善した候補だけをcurrentへ進め、失敗候補は必ず破棄する。

次の弱点：未知解の実問題では、解析真値を使うpost-solve受入はできません。G-1-B12では平衡化した主・随伴場、独立mesh cross-check、目的量の保証付き区間を組み合わせ、真値なしでaccept / stopできる契約へ進めます。

実装所有境界：NGSolveがP1/P3空間、`hesse`・`grad`、skeletonの`Other()`、HDiv/L2/NumberSpace、混合組立、UMFPACK、局所細分を所有し、Radiaがcohomologyと単位周期を供給します。JavaScriptはbasis、面向き、Piola写像、平衡化を再実装せず、検証済み11 solveを再生します。

実測可視化 G-1-B12 解析真値なしで目的量をacceptし、安全にstopする

G-1-B11の最後では、P1 probeが解析解へ偶然近いことを使って停止しました。しかし実機形状では真値がありません。そこでB12は、穴を一つもつT–Ω断面を厚み方向へ押し出して観察し、角度依存の未知透磁率

\mu_r(x,y) =1+8\exp\!\left[ -\dfrac{(x-0.023)^2+(y-0.012)^2}{2(0.009)^2} \right]

を入れます。中心導体の電流は100 A、外周で$\phi=0$、内周は自然境界です。穴を一周する単位cohomology場

\vec h_c=\dfrac{1}{2\pi(x^2+y^2)} \begin{pmatrix}-y\\x\end{pmatrix}, \qquad \oint_C\vec h_c\cdot d\vec{\ell}=1

を固定し、$\vec H_h=I\vec h_c-\nabla\phi_h$をNGSolveで解きます。$\mu_r$が角度方向に変わるため$\phi_h$はゼロではなく、目的量の解析真値も使えません。

主問題と随伴問題の両方を平衡化する

Gaussian観測を

J(\vec H)=\dfrac{\displaystyle\int_\Omega w\,H_y\,d\Omega} {\displaystyle\int_\Omega w\,d\Omega}, \qquad w=\exp\!\left[-\dfrac{(x-0.0275)^2+y^2}{2(0.003)^2}\right]

と固定します。随伴$z$は、同じ双線形形式$a(u,v)=\int_\Omega\mu\nabla u\cdot\nabla v\,d\Omega$を使って

a(v,z)=-\dfrac{1}{W}\int_\Omega w\,\dfrac{\partial v}{\partial y}\,d\Omega

を満たします。主場には法線連続な$\vec{\sigma}_u\in H(\mathrm{div})$、随伴には

\vec g_z=-\dfrac{w}{W}\vec e_y,\qquad \vec\rho_z=\vec g_z+\vec s_z

を構成します。$\nabla\cdot\vec{\sigma}_u=0$、$\nabla\cdot\vec s_z=0$、内周の自然境界で両補正の法線条件を満たすよう、NGSolveの$H(\mathrm{div})\times L^2$混合問題を解きます。すると

\eta_u^2=\int_\Omega (\mu\vec H_h-\vec{\sigma}_u)\cdot\mu^{-1} (\mu\vec H_h-\vec{\sigma}_u)\,d\Omega,

\eta_z^2=\int_\Omega (\mu\nabla z_h-\vec\rho_z)\cdot\mu^{-1} (\mu\nabla z_h-\vec\rho_z)\,d\Omega

が主・随伴のエネルギー誤差を上から押さえます。残差補正

R_u(z_h)=\ell(z_h)-a(\phi_h,z_h)

を含めると、Cauchy–Schwarzから目的量は

\boxed{ J(\vec H)\in \left[ J(\vec H_h)+R_u(z_h)-\eta_u\eta_z,\qquad J(\vec H_h)+R_u(z_h)+\eta_u\eta_z \right] }

に入ります。B11の$\xi_K$一つを目的誤差と読み替えたのではなく、主と随伴の二つの許容場の積で初めて上下界を作っています。

保証段階表示する証拠

独立P3 meshの値が保証区間内に入ったとき、独立meshは何を証明しますか。

区間の包含でacceptし、半幅でstopする

解析真値は保存せず、各候補について次だけを判定します。

1. $b_1=1$、円周積分100 A、主・随伴のdivergence許容性、実DOF予算を満たす。

2. 新しい保証区間が現在区間に含まれる。

3. 別の境界分割則・別のP3 meshで得た339.767 A/mが区間内にある。

4. 相対半幅が0.4%以下ならsafe stopする。

段階	要素 / 主DOF	$J_h$ [A/m]	保証区間 [A/m]	相対半幅	判定
base P1 / RT0	352 / 204	348.821	[14.651, 682.991]	95.80%	current
p-up P2 / RT1	1,190 / 2,478	339.856	[328.299, 351.413]	3.40%	accept
h = 2.5 mm	2,962 / 6,076	339.689	[337.033, 342.345]	0.782%	accept
h = 2.0 mm	4,626 / 9,444	339.646	[338.417, 340.874]	0.362%	accept + safe stop

独立meshはmaxh 3.0 mm・P3・10,449 DOFで、主系列と異なる境界分割則とmesh fingerprintをもちます。その339.767 A/mが全候補区間に入ることを確認しますが、これを証明の根拠には数えません。証明半径は常に$\eta_u\eta_z$です。

G-1-B12の受入順序

未知解化：角度依存$\mu_r(x,y)$を入れ、解析真値を判定データから完全に外す。
主平衡化：$\vec{\sigma}_u$の法線連続性、divergence、内周自然境界を検査する。
随伴平衡化：Gaussian荷重$\vec g_z$を保持したまま、発散ゼロ補正$\vec s_z$を解く。
区間化：$J_h+R_u(z_h)\pm\eta_u\eta_z$を作り、推定値と保証半径を分けて保存する。
accept：位相・予算・許容性・独立mesh照合を通り、前区間へ包含された候補だけを採用する。
safe stop：解析真値誤差ではなく、保証区間の相対半幅0.4%で停止する。

保証範囲：この区間は、表示した計算領域と固定14次求積で定義した変分モデルに対する相補性保証です。混合場の許容性は主側$10^{-14}$、随伴補正側$10^{-9}$以下で検査します。浮動小数点区間演算による丸め誤差証明や、理想円と離散幾何のモデル化誤差まで含む保証とは呼びません。

次の弱点：G-1-B13では求積、曲線境界、線形ソルバー、丸めの誤差予算を別項に分け、保証区間を「離散変分モデル」から「連続幾何モデル」へ拡張します。

実装所有境界：NGSolveがH1/HDiv/L2空間、Piola写像、主・随伴組立、混合平衡化、UMFPACK、別mesh生成を所有します。JavaScriptは有限要素basisや平衡化を再実装せず、5 mesh・782検査で凍結した区間、標識点、判定だけを再生します。

次に進む電磁界数値解析

弱形式部分積分、境界項、Galerkin法から有限要素方程式を作ります。 離散Hodgeと材料則接続行列と材料行列を分け、辺・面自由度を組み立てます。 電磁力と節点力Maxwell応力、エネルギー微分、節点力を比較します。 位相・切断面・ゲージtreeが消すgrad核と、穴を回って残すべき大域自由度を分けます。 保証付き誤差評価相補場からエネルギー誤差の上下界を作ります。 高次基底と適応解析h細分・p次数上昇・hp適応を選びます。 開境界とDtN作用素無限外部を境界作用素へ縮約し、弱形式へ戻します。

主な参照資料

高橋則雄『三次元有限要素法―磁界解析技術の基礎―』：三次元要素、補間関数、Galerkin法、節点要素・辺要素、境界条件、電磁力。
高橋則雄『磁気工学の有限要素法』：二次元静磁界の一次三角形要素、四面体辺要素、非線形磁気特性、渦電流、磁性材料モデリング。
高橋則雄ほか『電気工学の有限要素法』：要素係数行列、全体行列、面積座標、高次要素、固定・自然・周期境界条件。
高橋則雄『有限要素法を用いた最適化』：設計変数、目的関数、自動分割、要素数と最適解の関係。
五十嵐一・亀有昭久・加川幸雄・西口磯春・A. Bossavit『新しい計算電磁気学』：辺要素、向き、双対メッシュ、位相的離散化。