CAE-AI Lab 近畿大学 English

3Dで学ぶ電磁気学

G-2 微分形式による弱形式

F-6 Clebschポテンシャルまでは連続な場とポテンシャルを扱いました。G-2では、その場を有限要素法へ渡すため、静磁界の強形式をStokesの公式で弱形式へ変え、境界条件、ゲージ、Galerkin行列までを一続きに導きます。さらに F-5 Hodograph の連続作用素を受け取り、trial/test空間、Q1離散化、局所から全体への組立て、零空間診断まで実行します。次の G-3 Bossavitの離散Hodgeと材料則 で積分自由度と材料行列を深めます。

まず30秒でつかむ

式を読む順番

1. 強形式の左辺を残差 $R(\mathcal A)$ とみなし、任意の試験1-form $\mathcal W$ を掛けます。

2. 領域で積分し、Stokesの公式で未知場に掛かる微分を一段減らします。

3. その結果、領域内部の項と境界項が分かれ、Dirichlet条件と自然境界条件の役割が見えます。

4. 可視化で動かす矢印は解そのものではありません。試験関数の位置と向きを変え、どの局所残差が積分へ入るかを見ます。

この回の到達点

1. 強形式と弱形式が要求する滑らかさの違いを説明できる。

2. 外微分の積の公式とStokesの公式から境界項を導ける。

3. Dirichlet境界と自然境界を、試験1-formのトレースで区別できる。

4. ベクトルポテンシャルのゲージ核が離散行列へどう現れるか説明できる。

5. Whitney 1-formのGalerkin展開から $\boldsymbol K=\boldsymbol C^{\mathsf T}\boldsymbol M_\nu\boldsymbol C$ を組み立てられる。

6. 逆HodographのQ1最小二乗系について、境界互換性、局所scatter、低積分モード、平行移動核を別々に診断できる。

修士課程の電磁気学演習

弱形式の零空間は、物理的な零磁束と同じか

電磁気の問い: ベクトルポテンシャルへ勾配を加えると辺自由度は変わるが、面磁束と磁気エネルギーは変わらない。このゲージ核を残した行列と、tree-cotreeで除いた行列を区別できるか。

観察ポイント: ゲージ振幅、透磁率、拘束方法を変えたとき、辺未知数、$\boldsymbol C\boldsymbol a$、非零固有値、零空間次元のどれが変わるかを答える。

操作と観測: 可視化G-2-2で無拘束、1辺固定、tree-cotreeを切り替え、$\boldsymbol C\boldsymbol G=\boldsymbol0$、行列rank、零固有値、磁束、エネルギーを比較する。

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

離散系は $\boldsymbol K=\boldsymbol C^{\mathsf T}\boldsymbol M_\nu\boldsymbol C$ で、$\boldsymbol C\boldsymbol G=\boldsymbol0$ より $\boldsymbol a\mapsto\boldsymbol a+\boldsymbol G\phi$ は磁束を変えない。これは物理場が零という意味ではなく、同じ磁束を表すポテンシャルの冗長性である。

診断: ゲージ振幅で $\boldsymbol C\boldsymbol a$ やエネルギーが変わるなら、離散複体または境界条件が壊れている。正の材料Hodgeでもtree-cotree縮約後に零固有値が残るなら、位相自由度または拘束不足を調べる。

研究へ: ソルバーの特異性を数値的な失敗として隠さず、ゲージ核・cohomology・物理的零モードのどれかをrankと周期から分類する。

強形式：磁束密度からAmpereの法則へ

静磁界で磁気ベクトルポテンシャルを1-form $\mathcal A$ とすると、磁束密度2-formと磁界1-formは

\[ \mathcal B=d\mathcal A, \qquad \mathcal H=\star_{\nu}\mathcal B \]

です。ここで $\star_{\nu}$ は磁気抵抗率 $\nu=1/\mu$ と計量を含む材料Hodge starです。Ampereの法則 $d\mathcal H=\mathcal J$ へ代入すると、点ごとに満たす強形式は

\[ d\!\left(\star_{\nu}d\mathcal A\right)=\mathcal J \qquad\text{in }\Omega \]

となります。外側の $d$ まで直接評価するためには、$\mathcal A$ に強い滑らかさが必要です。辺要素を使う有限要素法では、試験1-form $\mathcal W$ を掛けて領域全体で積分し、外側の $d$ を試験関数へ移します。

外微分の積の公式とStokesの公式

1-form $\mathcal W$ と1-form $\mathcal H$ に対して、外微分の積の公式は

\[ d(\mathcal W\wedge\mathcal H) =d\mathcal W\wedge\mathcal H -\mathcal W\wedge d\mathcal H \]

です。これを領域 $\Omega$ で積分し、Stokesの公式

\[ \int_{\Omega}d(\mathcal W\wedge\mathcal H) =\int_{\partial\Omega} \operatorname{tr}\mathcal W\wedge\operatorname{tr}\mathcal H \]

を使います。$d\mathcal H=\mathcal J$ と $\mathcal H=\star_\nu d\mathcal A$ を代入して整理すると、

\[ \boxed{ \int_{\Omega}d\mathcal W\wedge\star_{\nu}d\mathcal A =\int_{\Omega}\mathcal W\wedge\mathcal J +\int_{\Gamma_N} \operatorname{tr}\mathcal W\wedge\operatorname{tr}\mathcal H_N } \]

を得ます。これが磁気ベクトルポテンシャルの弱形式です。外側の $d$ が試験1-formへ移ったため、$\mathcal A$ と $\mathcal W$ には1階の外微分だけを要求すればよくなります。

境界条件と試験空間

接線成分を与える境界を $\Gamma_D$、磁界を与える自然境界を $\Gamma_N$ とし、$\partial\Omega=\overline{\Gamma_D\cup\Gamma_N}$ と分けます。試験空間は

\[ \mathcal V =\left\{\mathcal W\in H\Lambda^1(\Omega) \;\middle|\; \operatorname{tr}\mathcal W=0\text{ on }\Gamma_D\right\} \]

です。既知の境界値をもつ持ち上げ $\mathcal A_D$ を使えば、未知量は

\[ \mathcal A\in\mathcal A_D+\mathcal V, \qquad \mathcal W\in\mathcal V \]

と書けます。$\Gamma_D$ では試験関数のトレースが消え、$\Gamma_N$ の境界積分だけが右辺に残ります。

ゲージ核と変分原理

$d^2=0$ なので、任意の0-form $\phi$ に対して

\[ d(\mathcal A+d\phi)=d\mathcal A, \qquad \mathcal A\sim\mathcal A+d\phi \]

です。この勾配成分は磁束密度を変えないため、Coulombゲージ、tree-cotreeゲージ、または商空間で核を除かなければ離散行列は特異になります。大域的な核と穴を回る循環はG-6で扱います。

以下の2次汎関数による議論では、$\star_\nu$ が線形・対称・正定値であることを仮定します。非線形材料では $\dfrac12\mathcal B\wedge\star_\nu\mathcal B$ をそのまま使わず、構成則を積分した磁気エネルギー密度へ置き換えます。同じ弱形式は、線形材料の磁気ポテンシャル汎関数

\[ \Pi(\mathcal A) =\dfrac12\int_{\Omega}d\mathcal A\wedge\star_{\nu}d\mathcal A -\int_{\Omega}\mathcal A\wedge\mathcal J -\int_{\Gamma_N} \operatorname{tr}\mathcal A\wedge\operatorname{tr}\mathcal H_N \]

の停留条件

\[ \delta\Pi(\mathcal A;\mathcal W)=0 \qquad\text{for every }\mathcal W\in\mathcal V \]

からも得られます。したがって、弱形式の左辺は磁気エネルギーの2次形式、右辺は電流源と自然境界が行う仕事として読めます。

Whitney-Galerkin行列

Whitney 1-form $w_e^1$ で未知量と試験関数を

\[ \mathcal A_h=\sum_e a_e w_e^1, \qquad \mathcal W=w_{e'}^1 \]

と展開すると、Galerkin行列と右辺は

\[ K_{e'e} =\int_{\Omega}dw_{e'}^1\wedge\star_{\nu}dw_e^1, \qquad f_{e'} =\int_{\Omega}w_{e'}^1\wedge\mathcal J +\int_{\Gamma_N}\operatorname{tr}w_{e'}^1\wedge\operatorname{tr}\mathcal H_N \]

となります。外微分を表す接続行列 $\boldsymbol C$ と、材料・計量を表す離散Hodge行列 $\boldsymbol M_\nu$ を分けると、

\[ \boldsymbol K\boldsymbol a=\boldsymbol f, \qquad \boldsymbol K=\boldsymbol C^{\mathsf T}\boldsymbol M_{\nu}\boldsymbol C \]

です。$\boldsymbol C$ はメッシュの接続だけ、$\boldsymbol M_\nu$ は形状・計量・材料だけを担います。

可視化 G-2-1 試験1-formが拾う局所寄与

下図は3次元辺要素の局所台を平面で切った教育用モデルです。橙色の向き付き辺 $e$ に対応するWhitney 1-form $w_e^1$ は、辺に接する2要素 $K_1,K_2$ の内部だけで非零になります。緑の矢印は幾何と接続から決まる $w_e^1$ なので、材料を変えても動きません。一方、面の濃さと右側の棒は局所積分密度 $\nu|dw_e^1|^2$ と要素寄与を表します。要素 $K_2$ の比透磁率を変え、基底は不変でも材料Hodgeによる重みが変わることを確認してください。

μr=6.0

各三角形の面積を $1/2$、$|dw_e^1|=2$ と規格化したモデルでは、自己剛性の局所寄与は

\[ K_{ee}^{(1)}=2\nu_1, \qquad K_{ee}^{(2)}=2\nu_2=\dfrac{2}{\mu_r}, \qquad K_{ee}=K_{ee}^{(1)}+K_{ee}^{(2)} \]

です。試験関数の台が局所的なので、全領域積分は要素積分の和として組み立てられます。辺や要素の向きが変われば $dw_e^1$ の符号は変わっても、自己エネルギー $dw_e^1\wedge\star_\nu dw_e^1$ は非負のままです。

可視化 G-2-2 修士演習：tree-cotreeでゲージ核を除く

最初のイメージ。 メッシュの各辺に「その辺に沿って測った磁気ベクトルポテンシャルの値」を書き込むと考えます。節点の高さを全体へ足して辺差分を変えても、一つの三角形を一周した合計は変わりません。この見えない自由度がゲージです。全域木の辺を基準として固定すると、閉路を作るcotree辺だけが物理的な周回を担い、重複した自由度を除けます。

可視化G-2-1では、一つの辺基底が二つの要素へ与える局所寄与を見ました。ここでは9節点、16辺、8三角形からなるメッシュ全体を組み立て、局所行列から大域行列へ進みます。節点値を辺差分へ写す接続行列を $\boldsymbol G$、辺周回を面へ写す接続行列を $\boldsymbol C$ とすると、

\[ \boldsymbol G\in\mathbb R^{16\times9}, \qquad \boldsymbol C\in\mathbb R^{8\times16}, \qquad \boldsymbol C\boldsymbol G=\boldsymbol0 \]

です。最後の式は「勾配を一周積分すると0」、すなわち $d^2=0$ の離散版です。磁気抵抗率を面ごとの正の重みとして離散Hodge行列 $\boldsymbol M_\nu$ に入れると、弱形式の大域行列は

\[ \boldsymbol K =\boldsymbol C^{\mathsf T}\boldsymbol M_\nu\boldsymbol C, \qquad \boldsymbol K\boldsymbol G\boldsymbol\phi =\boldsymbol0 \]

となります。したがって、節点ポテンシャル $\boldsymbol\phi$ から作る勾配 $\boldsymbol g=\boldsymbol G\boldsymbol\phi$ は、辺自由度を変えても面磁束を変えないゲージモードです。

無拘束では $\operatorname{rank}\boldsymbol C=8$ に対して辺未知数が16個あるため、

\[ \dim\ker\boldsymbol K =16-\operatorname{rank}\boldsymbol C =8 \]

です。1辺だけ固定しても未知数は15個なので零モードは7個残ります。赤破線で示す8本の全域木辺を固定し、残る8本のcotree辺だけを未知量にすれば、縮約行列は

\[ \boldsymbol K_c =\boldsymbol C_c^{\mathsf T}\boldsymbol M_\nu\boldsymbol C_c \in\mathbb R^{8\times8}, \qquad \ker\boldsymbol K_c=\{\boldsymbol0\} \]

となります。スライダーで $t$ を動かすと紫の辺自由度は変わりますが、面色と面内の $B_f=(\boldsymbol C\boldsymbol a)_f$ は変わりません。実際、

\[ \boldsymbol C(\boldsymbol a+t\boldsymbol g_j) =\boldsymbol C\boldsymbol a, \qquad \dfrac12(\boldsymbol a+t\boldsymbol g_j)^{\mathsf T} \boldsymbol K(\boldsymbol a+t\boldsymbol g_j) =\dfrac12\boldsymbol a^{\mathsf T}\boldsymbol K\boldsymbol a \]

だからです。一方、$\mu_r$ を変えると $\boldsymbol M_\nu$ が変わるため、青い非零スペクトルと磁気エネルギーは変化します。しかし $\boldsymbol C$ の接続は変わらないので、赤い零モードの本数は変わりません。位相が決める零空間と、材料・計量が決める非零固有値を同じ図で区別してください。

動かして確かめる順序。 まず「無拘束」でゲージ振幅 $t$ を動かし、紫の辺値は変わるのに面磁束とエネルギーが変わらないことを確認します。次に「1辺だけ固定」を選び、零モードが一つしか減らない理由を未知数とrankから数えます。「tree-cotree」へ切り替える前に零モード数を確認し、縮約行列が正定値になることをスペクトルで確かめます。最後に $\mu_r$ を変え、非零固有値は動く一方で $\boldsymbol C\boldsymbol G=0$ と零空間の次元が材料に依存しないことを説明してください。

G-2-H1 Chaplygin方程式の弱形式

F-5の連続Hodographで導いた強形式 $\mathcal L_hA=0$ を、ここからtrial/test空間、自然境界項、正値な双線形形式へ移します。材料則は既知係数 $a(q),b(q)$ に入り、有限要素法では境界条件と近似空間を明示して解きます。

境界を $\partial\Omega_h=\Gamma_D\cup\Gamma_N$ と分け、

\[ \begin{aligned} \boldsymbol K_h(q)&= \begin{pmatrix}a(q)&0\\0&b(q)\end{pmatrix},\\ \nabla_h A&= \begin{pmatrix}A_q\\A_\theta\end{pmatrix} \end{aligned} \]

を使って、

\[ \begin{aligned} A&=A_D &&\text{on }\Gamma_D,\\ \vec n_h\cdot\boldsymbol K_h\nabla_h A &=g_N &&\text{on }\Gamma_N \end{aligned} \]

を与えます。ここで $\vec n_h=(n_q,n_\theta)$ はHodograph平面の外向き単位法線です。試験空間と解空間を

\[ \begin{aligned} V_0&=\{v\in H^1(\Omega_h)\mid v=0\text{ on }\Gamma_D\},\\ V_{A_D}&=\{w\in H^1(\Omega_h)\mid w=A_D\text{ on }\Gamma_D\} \end{aligned} \]

とすると、弱形式は「$A\in V_{A_D}$ を求め、任意の $v\in V_0$ に対して」

\[ \boxed{ \begin{aligned} &\int_{\Omega_h}a(q)A_qv_q\,dq\,d\theta\\ &\quad+\int_{\Omega_h}b(q)A_\theta v_\theta\,dq\,d\theta\\ &=\int_{\Gamma_N}g_Nv\,ds_h \end{aligned} } \]

となります。$\mu(q)>0$ かつ $\{\mu(q)q\}'>0$ なら $a(q)>0$、$b(q)>0$ なので、この双線形形式は対称かつ正です。したがって、材料飽和は物理空間では非線形係数として現れる一方、$(q,\theta)$ 平面では既知の正係数を持つ線形楕円問題として扱えます。

G-2-H2 逆Hodographの一階Galerkin形と最小二乗弱形式

座標を与える境界を $\Gamma_D$ とし、

\[ \begin{aligned} \boldsymbol V_D &=\left\{ \vec r\in[H^1(\Omega_h)]^2 \mid \operatorname{tr}\vec r=\vec r_D \text{ on }\Gamma_D \right\},\\ \boldsymbol V_0 &=\left\{ \vec w\in[H^1(\Omega_h)]^2 \mid \operatorname{tr}\vec w=\vec0 \text{ on }\Gamma_D \right\} \end{aligned} \]

とします。最も直接的な一階Galerkin形は「$\vec r\in\boldsymbol V_D$ を求め、任意の $\vec v\in[L^2(\Omega_h)]^2$ に対して」

\[ \int_{\Omega_h} \vec v\cdot \mathcal D_{\mathrm C}\vec r \,dV_h =0 \]

です。Stokesの公式で微分を試験関数側へ移せば、

\[ \begin{aligned} &\int_{\Omega_h} \vec r\cdot \left[ -\vec v_\theta +\dfrac{\partial}{\partial q} \left(\boldsymbol C^{\mathsf T}\vec v\right) \right]dV_h\\ &\quad+ \int_{\partial\Omega_h} \left[ n_\theta\,\vec v\cdot\vec r -n_q\, \left(\boldsymbol C^{\mathsf T}\vec v\right)\cdot\vec r \right]ds_h =0 \end{aligned} \]

となります。これは座標を $L^2$ まで下げ、境界traceを別未知量として扱うultraweak形式へ進む入口です。

実装する最小二乗弱形式では、汎関数

\[ \mathcal J_{\mathrm{LS}}[\vec r] =\dfrac12\int_{\Omega_h} \left|\mathcal D_{\mathrm C}\vec r\right|^2dV_h \]

を $\boldsymbol V_D$ 上で最小化します。すなわち、$\vec r\in\boldsymbol V_D$ を求め、任意の $\vec w\in\boldsymbol V_0$ に対して

\[ \int_{\Omega_h} \mathcal D_{\mathrm C}\vec w\cdot \mathcal D_{\mathrm C}\vec r\,dV_h=0 \]

とします。これは残差を零にできない境界データにも最良近似を定義でき、次のQ1組立てへ直接移せます。

ここまでの式を「読んだ」で終わらせず、実際に離散化して解きます。次の演習では厳密解を描画結果として使いません。厳密解は境界値と誤差評価だけに使い、内部節点の $x,y$ はブラウザ内で組み立てた連立方程式から求めます。

可視化 G-2-H3 逆Hodograph最小二乗を離散化する

Hodograph長方形 $\Omega_h=[q_{\min},q_{\max}]\times[\theta_{\min},\theta_{\max}]$ を $N\times N$ 節点の双一次格子に分け、

\[ \vec r_h(q,\theta) =\sum_{a=1}^{N^2}N_a(q,\theta)\vec R_a, \qquad \vec R_a= \begin{pmatrix}x_a\\y_a\end{pmatrix} \]

と近似します。一つのセル $K$ の四隅を $\vec R_{00},\vec R_{10},\vec R_{01},\vec R_{11}$ とすれば、セル中心の微分は

\[ \begin{aligned} \left.\vec r_{q,h}\right|_K &=\dfrac{(\vec R_{10}+\vec R_{11})-(\vec R_{00}+\vec R_{01})}{2\Delta q},\\[3pt] \left.\vec r_{\theta,h}\right|_K &=\dfrac{(\vec R_{01}+\vec R_{11})-(\vec R_{00}+\vec R_{10})}{2\Delta\theta} \end{aligned} \]

です。各セルで

\[ \vec R_K =\left.\vec r_{\theta,h}\right|_K -\boldsymbol C(q_K,\theta_K) \left.\vec r_{q,h}\right|_K \]

を作り、全セルの残差を行列 $\boldsymbol A_h$ へ集めます。境界節点の既知寄与を $\vec c_h$、内部節点の未知座標を $\vec u_h$ と書くと、画面内で解いている式は

\[ \boxed{ \boldsymbol A_h^{\mathsf T}\boldsymbol A_h\vec u_h =-\boldsymbol A_h^{\mathsf T}\vec c_h } \]

です。これは上の $\mathcal J_{\mathrm{LS}}$ の停留条件を、セル中心積分した正規方程式です。材料指数 $k$ を動かすと $\beta=q/(k+1)$ と $\boldsymbol C$ が変わるため、同じHodograph長方形でも物理格子は別の形になります。

境界条件の互換性も動かして確かめます。$\delta_\Gamma=0$ では、べき乗材料 $\mu(q)=q^k$ の厳密逆写像

\[ \begin{aligned} x_{\mathrm{ex}}(q,\theta) &=-(k+1)\ln q+\dfrac{k}{2}\sin^2\theta,\\ y_{\mathrm{ex}}(q,\theta) &=\left(1+\dfrac{k}{2}\right)\theta -\dfrac{k}{4}\sin 2\theta \end{aligned} \]

を全周のDirichlet値に使います。$\delta_\Gamma>0$ では $q=q_{\max}$ の辺だけを、角を固定したまま

\[ \vec r_D^{\,\delta}(s) =\vec r_{\mathrm{ex}}(s) +\delta_\Gamma\sin(\pi s)\vec e_H, \qquad 0\le s\le1 \]

だけ押し出します。これは、第一階系と両立しない物理境界を全周へDirichlet指定したときに何が起こるかを見るための、意図的な過剰拘束です。

境界節点を指定できることと、内部で $\vec R_K=0$ を満たせることを分けて図を動かして確かめてください。

境界条件を満たすことと、内部方程式を満たすことを分ける

左の色は全状態で共通の $10^{-3}\le|\vec R_K|\le1$ 対数尺度です。操作後は灰色の互換境界解を固定し、不整合境界の現在解と同じ物理尺度で重ねます。赤セルは期待する $J\lt0$ から符号が変わったセルです。

読み方は三段階です。

1. まず $\delta_\Gamma=0$ のまま $m$ を上げます。$N=5,7,9,11$ と細分され、RMS残差とRMS座標誤差が下がることを確かめます。

2. 次に $k$ を動かします。方程式は $x,y$ に関して線形のままですが、材料係数 $\boldsymbol C$ と物理格子が変わります。

3. 最後に $\delta_\Gamma$ を上げます。細分しても残る最小二乗残差と、$\min|J|$ の低下を見ます。これは「線形化できた」ことと「任意の境界形状を自由に指定できる」ことが同義でない理由です。

画面の診断量は

\[ \varepsilon_{\mathrm{LS}} =\sqrt{\dfrac{1}{2|\mathcal T_h|} \sum_{K\in\mathcal T_h}|\vec R_K|^2}, \qquad J_{\min}=\min_{K\in\mathcal T_h} \left|\det(\vec r_{q,h},\vec r_{\theta,h})\right| \]

です。$J_{\min}$ が0へ近づけば、残差が小さくても逆写像が退化しかけています。実設計には厳密座標誤差を計算できないため、$\varepsilon_{\mathrm{LS}}$、$J_{\min}$、物理空間FEMとの場の差を組み合わせて判定します。

可視化 G-2-H4 一つのセルから全体行列へ

前のCanvasは連立方程式を解きましたが、ここでは一つのセル $K$ が何を全体行列へ足したかを分解します。セルの四節点座標を並べた局所ベクトルを $\vec u_K$、内部自由度だけを集めた全体ベクトルを $\vec u$ とすると、局所と全体の関係は

\[ \vec u_K =\boldsymbol S_K\vec u+\vec u_{D,K} \]

です。$\boldsymbol S_K$ はセルが参照する内部自由度を選ぶBoolean散布行列で、$\vec u_{D,K}$ は境界節点の既知座標です。セル中心の二成分残差を

\[ \vec R_K =\boldsymbol A_K\boldsymbol S_K\vec u+\vec c_K, \qquad \vec c_K=\boldsymbol A_K\vec u_{D,K} \]

と書けば、セルの最小二乗エネルギーと全体汎関数は

\[ \mathcal J_K =\dfrac{1}{2}|\vec R_K|^2\,\Delta q\Delta\theta, \qquad \mathcal J_h=\sum_{K\in\mathcal T_h}\mathcal J_K \]

です。微分して得る全体行列と右辺は

\[ \boxed{ \begin{aligned} \boldsymbol K_h &=\sum_K \boldsymbol S_K^{\mathsf T} \boldsymbol A_K^{\mathsf T}\boldsymbol A_K \boldsymbol S_K\,\Delta q\Delta\theta,\\ \vec f_h &=-\sum_K \boldsymbol S_K^{\mathsf T} \boldsymbol A_K^{\mathsf T}\vec c_K \,\Delta q\Delta\theta \end{aligned} } \]

となります。次のCanvasは $N=5$、内部18自由度の実行列を描きます。左では緑の $\vec r_{\theta,h}$ と青破線の $\boldsymbol C\vec r_{q,h}$ が一致すれば残差が0です。二つの矢印の先端を結ぶ赤線が、実際の $\vec R_K$ です。右は全体の $\boldsymbol K_h$ で、橙枠が選択セルから加算された成分です。

境界上の灰色×は既知のDirichlet節点です。自由な節点数から、座標未知数と橙枠の候補数を図を動かして確かめてください。

局所残差、散布、全体正規行列

左の物理格子をクリックしてもセルを選択できます。行列の青は正、茶は負、白は数値的に0です。非零パターンは格子の接続で決まり、値は材料指数と計量係数で変わります。

操作では次を確認してください。

1. セル番号を動かすと、橙枠だけが全体行列の別の場所へ移ります。これは $\boldsymbol S_K$ が変わるためです。

2. $k$ を動かしても行列の疎構造は変わりませんが、青・茶の濃さと矢印 $\boldsymbol C\vec r_q$ は変わります。接続と材料則が別の役割を持つことを示します。

3. $\delta_\Gamma$ を0へ戻すと、格子細分誤差だけが残ります。不整合を増やすと境界近傍の $|\vec R_K|$ と $\mathcal J_K$ が増え、その局所寄与が全体解へ伝わります。

この小系では正規方程式をCholesky分解して構造を見せています。ただし一般に

\[ \kappa_2\!\left(\boldsymbol A_h^{\mathsf T}\boldsymbol A_h\right) =\kappa_2(\boldsymbol A_h)^2 \]

なので、大規模・高次・強い材料異方性では条件数を二乗する正規方程式を無条件に選びません。QR分解、LSQR、適切な前処理、または一階混合形式を比較します。Canvasの「Cholesky対角比」は状態を眺める簡易指標であり、厳密なスペクトル条件数ではありません。

可視化 G-2-H5 境界条件と零空間を見分ける

全体行列を組めても、境界条件を与える前の逆写像は一意ではありません。$N\times N$ 節点の全座標を

\[ \vec u_0=(x_1,y_1,\ldots,x_{N^2},y_{N^2})^{\mathsf T} \]

とし、全セルの一階残差を積んだ設計行列を $\boldsymbol A_0$ と書きます。$N=5$、セル中心1点だけで積分すると

\[ \boldsymbol A_0\in\mathbb R^{32\times50}, \qquad \dim\ker\boldsymbol A_0 =50-\operatorname{rank}\boldsymbol A_0\ge18 \]

です。方程式数だけを数えても、境界上の写像を決めなければ多数の座標変形が残ることが分かります。拘束後に残す自由度を選ぶBoolean行列を $\boldsymbol P_F$ とすれば、実際に調べる行列は

\[ \boldsymbol A_F=\boldsymbol A_0\boldsymbol P_F, \qquad \ker\boldsymbol A_F=\{\vec0\} \Longleftrightarrow \sigma_{\min}(\boldsymbol A_F)>0 \]

です。境界条件は局所方程式を変えるのではなく、許される試行空間を狭めます。

もう一つ注意があります。双一次要素の残差をセル中心1点だけで積分すると、セル内で符号が入れ替わる変形を見落とすことがあります。それは低積分が作る零エネルギーモードです。そこで同じQ1要素を $2\times2$ Gauss点でも評価し、低積分モードと平行移動モードを分けるとともに、境界値を決めていないために残る自由度を比較します。零空間ベクトル $\vec v$ に対しては

\[ \boldsymbol A_F\vec v=\vec0, \qquad \left\|\boldsymbol A_F(\vec u+t\vec v)+\vec c\right\|_2 =\left\|\boldsymbol A_F\vec u+\vec c\right\|_2 \]

なので、$t$ を動かすと物理格子は変形しても離散残差は増えません。これが「行列は組めたが形状が決まらない」状態です。

中心1点・拘束なしでは零空間が18次元あります。積分不足と平行移動を、どの操作が取り除くか図を動かして確かめてください。

零空間モードと特異値スペクトル

左は灰色の基準格子へ実際の零空間モードを加えた格子です。赤点は固定自由度、緑矢印は節点変位です。右の赤い棒は数値的な零特異値を表します。

次の順に操作してください。

1. 「中心1点・4隅」のまま $j$ と $t$ を動かします。4隅を固定しても、セル中心だけでは見えない格子変形が残ります。

2. 積分を「$2\times2$ Gauss」へ変えます。Q1要素のセル内変化を拾うため、低積分由来のモードが減ります。

3. 拘束を「拘束なし」「左下1点」「4隅」「周囲境界」の順に変え、$\operatorname{rank}\boldsymbol A_F$ と $\dim\ker\boldsymbol A_F$ を比較します。周囲境界で零空間が消えたときだけ離散最小二乗解が一意になります。

4. $k$ を動かします。接続と拘束が同じなら行列寸法は変わりませんが、非零特異値の分布と条件の悪さは材料係数に応じて変化します。

初期値 $N=5$ では、「中心1点・4隅固定」が $32\times42$、rank 32、nullity 10です。同じ4隅固定を $2\times2$ Gaussへ変えると $128\times42$、rank 42、nullity 0になります。一方、Gauss積分でも無拘束なら $128\times50$、rank 48、nullity 2で、残る2本は物理格子全体の $x$・$y$ 平行移動です。周囲境界を固定すると $128\times18$、rank 18、nullity 0になります。これらはCanvas内で毎回行基本変形して得る値であり、表示用に決め打ちした数ではありません。

零空間が0でも、それだけで良い磁極形状とは限りません。最小非零特異値が極端に小さい場合は境界誤差が大きく増幅され、$J=0$ に近い写像を生む可能性があります。したがって、実装ではrank、特異値、最小Jacobian、物理空間FEMとの場の差を同時に監視します。

ただし、F-5の固定磁極可視化はこの弱形式を既知の $\Omega_h$ 上で直接解いた結果ではありません。物理空間非線形FEMの解を $(x,y)\mapsto(\theta,q)$ へ写しています。固定磁極が作る $\partial\Omega_h$ を再構成し、上の弱形式を直接解いて物理空間FEMと照合することが、ここから先のHodograph研究課題です。

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

接続行列は境界の境界が空であることを離散的に表し、

\[ \boldsymbol C\boldsymbol G=\boldsymbol0, \qquad \boldsymbol D\boldsymbol C=\boldsymbol0 \]

を満たします。材料が正なら離散Hodgeも正であり、ゲージ核を除いた空間で

\[ \boldsymbol a^{\mathsf T}\boldsymbol K\boldsymbol a =(\boldsymbol C\boldsymbol a)^{\mathsf T} \boldsymbol M_\nu(\boldsymbol C\boldsymbol a)\ge0 \]

です。位相的な恒等式と材料エネルギーの正値性を別々に検査できることが、この分解の利点です。次のG-3 Bossavitの離散Hodgeと材料則では、積分自由度、Whitney形式、離散Hodgeを詳しく扱い、リッツ線の異方性材料則へ進みます。

可視化 G-2-3 Stokes境界項の向き

1-form $\alpha=x\,dy$ なら

\[ d\alpha=dx\wedge dy \]

です。向き付き曲面 $S$ と、その向きから誘導される境界 $\partial S$ に対して

\[ \int_S d\alpha = \oint_{\partial S}\alpha \]

が成立します。三角形の頂点を $z$ 方向へ持ち上げると実面積 $|S|$ は増え、単位法線の $z$ 成分 $n_z$ は減りますが、その積 $n_z|S|$ は $xy$ 平面への向き付き投影面積なので変わりません。これは不具合ではなく、2-form $dx\wedge dy$ が測っている量そのものです。頂点の $y$ 座標も動かすと投影面積が変わり、面積分と境界線積分が同じ値で追従します。面の向きを反転すると両者の符号が同時に反転します。

可視化 G-2-4 試験関数の台を動かして要素寄与を追う

有限要素の双線形形式は、要素積分の和として組み立てます。

\[ a(\mathcal A,w_i) =\sum_{K\in\mathcal T_h} \int_K\nu\,d\mathcal A\wedge\star dw_i \]

局所台をもつ試験関数では、$w_i|_K=0$ の要素は行列成分へ寄与しません。

\[ \operatorname{supp}w_i =\bigcup\{K\in\mathcal T_h\mid w_i|_K\ne0\} \]

下の規格化実験では、中心から外周へ向く辺 $e_i$ を8三角形の扇形メッシュで共有させます。添字は周期的に $K_9=K_1$ と読み、

\[ \operatorname{supp}w_i=K_i\cup K_{i+1} \]

です。材料重みの山を

\[ m_k(c)=1+c\exp\!\left[-0.55(k-6)^2\right], \qquad k=1,\ldots,8 \]

とし、幾何係数を1に規格化すると、各要素の自己寄与と大域和は

\[ a_K(w_i,w_i) =\boldsymbol1_{K\subset\operatorname{supp}w_i}\,m_K(c), \qquad a_i=\sum_{K=1}^{8}a_K(w_i,w_i) \]

になります。したがって、辺を動かすと非零になる要素の場所は移りますが、その個数は内点辺なら2のままです。材料は台の形を変えず、拾った二要素の寄与の大きさだけを変えます。

橙色の台と、右図の非零棒の本数を先に見てください。

初期状態では $e_4$、$c=4$ に固定します。操作後は灰色破線の基準台を残し、橙色の現在台、8本の要素寄与、非零数、大域和を同じ尺度で比較します。辺と材料コントラストを別々に動かし、「接続が零・非零の場所を決め、材料が非零値の大きさを決める」ことを確かめてください。

可視化 G-2-5 体積項と境界項を分けて読む

部分積分後の弱形式を、抽象的な領域の絵ではなく実際のWhitney辺形状関数で確かめます。参照三角形

\[ K=\{(x,y)\mid x\ge0,\ y\ge0,\ x+y\le1\} \]

の重心座標を $\lambda_1=1-x-y$、$\lambda_2=x$、$\lambda_3=y$ とし、向き付き辺 $e_{ij}$ のWhitney 1-formを

\[ w_{ij}=\lambda_i\,d\lambda_j-\lambda_j\,d\lambda_i \]

とします。三辺を反時計回りに $e_{12},e_{23},e_{31}$ と向けると、どの基底も $dw_{ij}=2\,dx\wedge dy$ ですが、赤い自然境界 $\Gamma_N=e_{12}$ 上の接線トレースは

\[ I_e=\int_{\Gamma_N}\operatorname{tr}w_e =\begin{cases} 1,&e=e_{12},\\ 0,&e=e_{23},e_{31} \end{cases} \]

となります。つまり境界が赤く描かれているだけでは、すべての試験関数へ同じ境界荷重が入るわけではありません。その境界に自由度を持つ基底だけが自然境界項を拾います。

\[ R_e =2\nu a_e-f_{J,e}-h_N I_e \]

このCanvasでは $\nu=1$、電流源の局所荷重 $f_{J,e}=0.40$ とし、選んだ基底だけを使う一自由度モデルで $R_e=0$ となる $a_e$ を計算します。全要素行列では他の辺基底との相互項も加わりますが、ここでは「どの境界項がこの行へ入るか」に焦点を絞ります。

赤い Γₙ 上で接線積分が1になる基底を図を動かして確かめてください。 初期状態では、形状関数と局所行の答えを隠しています。

動かして確かめる順序。 まず赤い $\Gamma_N=e_{12}$ を見て、$w_{12},w_{23},w_{31}$ のどれが境界積分を持つか確かめます。操作後は、橙の選択辺、紫の形状関数、$I_e$、体積項 $2\nu a_e$、電流源項 $f_{J,e}$、自然境界項 $h_NI_e$ を同じ色で追ってください。次に $h_N$ を変え、$w_{12}$ の局所行だけが応答し、他の二行では境界項が0のままになる理由を接線トレースから説明してください。

理解の確認

1. 強形式から弱形式へ移ると、未知1-formに要求する滑らかさが弱くなるのはなぜでしょうか。

2. $\Gamma_D$ で消える試験関数と、$\Gamma_N$ に残る境界積分は、それぞれどの境界条件を表すでしょうか。

3. $\mathcal A\mapsto\mathcal A+d\phi$ の不定性を残したまま行列を解くと、何が起こるでしょうか。

4. $\boldsymbol C$ と $\boldsymbol M_\nu$ を分けると、材料変更とメッシュ接続変更をどう区別できるでしょうか。

5. 可視化で $\mu_r$ を大きくしたとき、$K_{ee}^{(2)}$ が小さくなる理由を磁気抵抗率から説明してください。

6. Chaplygin方程式の強形式をF-5に残し、G-2では境界項とtrial/testから始めると、連続理論と離散化のどこを独立に検査できるでしょうか。

7. 逆Hodograph最小二乗でDirichlet節点が指定境界に一致しても、内部セル残差が零とは限らない理由を説明してください。

8. Q1残差行列で、並進mode、境界条件不足、低減積分のzero-energy modeをrank・nullity・特異値からどう見分けますか。

参考: A. Bossavit, *Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements*, Academic Press (1998).