COMSOL加速空洞: 進行波型加速管(TMモード)偏微分方程式からモデルビルダー

進行波加速管では，軸対称の TM$_{01}$ モードを加速モードとして使います．このモードでは，電場は $E_\phi=0$ で，$E_r$ と $E_z$ が有限の値をとります．一方，磁場は $H_r=H_z=0$ で，$H_\phi$ のみが残ります．したがって，$(E_r,\,E_z)$ あるいは $H_\phi$ のいずれかを計算すれば，Maxwell 方程式から残りの場を求めることができます．未知数が一つで済む $H_\phi$ を計算する方が簡単なので，ここで示す COMSOL モデルでは $H_\phi$ を未知関数として計算します．

進行波加速管の加速セル

図1と表 1 に，計算モデルと寸法を示します．これは進行波型加速管の 1 周期分です．表 1 のパラメーターには，図1に示す寸法のほか，COMSOL の計算で使う値も含めています．cell_phase は左端と右端の RF 位相差 $\phi$，freq_search_min は固有周波数探索の開始周波数，omega_search_min はその角周波数です．N_cal_modes は求めるモード数です．

位相規約

このページでは，RF 場の時間依存性を $\exp(-i\omega t)$ とし，z 方向に進む波を $\exp(i k_z z)$ で表します．以降，$H_\phi(r,z)$ は時間因子を除いた複素振幅を表すものとします．時間依存性を含めた複素表示と 1 周期あたりの位相進みは，次のように書けます．

\begin{align} H_\phi(r,z,t)&=H_\phi(r,z)\exp(-i\omega t) \label{eq:phase_convention}\\[6pt] H_\phi(r,z+L)&=H_\phi(r,z)\exp(i k_z L) \label{eq:Hphi_phase_advance}\\[6pt] k_z&=\frac{\phi}{L}=\frac{\mathrm{cell\_phase}}{\mathrm{cell\_D}} \label{eq:kz_definition} \end{align}

ここで，$L$ は 1 周期の長さ，$\phi$ は 1 周期あたりの位相進みです．COMSOL の設定では，表 1 のように $L=\mathrm{cell\_D}$，$k_z=\mathrm{cell\_phase}/\mathrm{cell\_D}$ とします．この符号を使うと，CST Studio の phase advance と同じ位相定義になります．実数の場が必要な場合は，最後にこの複素表示の実部を取ります．この節では $H_\phi$ の位相だけを定義し，数値計算用の変数変換は後の節で導入します．

ここでは，軸対称の進行波型加速管 1 セルの固有角周波数と電磁場を COMSOL で計算するために，必要な式を整理します．COMSOL の Coefficient Form PDE では，

\begin{align} &\lambda^2e_au+\lambda d_a u+\acute{\nabla}\cdot\left(-c\acute{\nabla} u-\alpha u+\gamma\right)+\beta\cdot\acute{\nabla}u + au =f& \label{eq:COMSOL_PDE}\\[6pt] &\acute{\nabla}=\left[\frac{\partial}{\partial r},\,\frac{\partial}{\partial z}\right]& \label{eq:COMSOL_solve_PDE} \end{align}

を解きます．ユーザーは，式\eqref{eq:COMSOL_PDE}の係数 ($e_a,\,d_a,\,c,\,\alpha,\,\gamma,\,\beta,\,a$) とソース項 $f$ を設定します．さらに境界条件を指定すると，固有値 $\lambda$ と固有関数 $u$ が得られます．このページでは $\lambda=\omega$ として扱い，周波数は $\operatorname{Re}(\lambda)/(2\pi)$ で求めます．なお，COMSOL の係数 $c$ と光速は記号が重なります．以下では，光速を $c_0$ と書き，COMSOL では c_const を使います．

ここでは，図2の領域を計算します．これは，図1に示した軸対称加速セルを，二次元座標系 $(r,\,z)$ で表したものです．以降はこの座標系を使います．

物理的な偏微分方程式

軸対称構造の加速空洞内の電磁場を表す偏微分方程式は，

\begin{align} & \nabla^2\boldsymbol{E}+\left(\frac{\omega}{c_0}\right)^2\boldsymbol{E}=0 & & \nabla^2\boldsymbol{H}+\left(\frac{\omega}{c_0}\right)^2\boldsymbol{H}=0 & \end{align}

です．軸対称 TM$_{0\,n}$ モードの場合，これらのベクトル方程式は，成分ごとの偏微分方程式として次のように書けます．

\begin{align} \def\pd#1#2{\frac{\partial #1}{\partial #2}} \def\pdd#1#2{\frac{\partial^2 #1}{\partial #2^2}} & \pdd{E_r}{r}+\frac{1}{r}\pd{E_r}{r}+\pdd{E_r}{z}-\frac{E_r}{r^2}+\left(\frac{\omega}{c_0}\right)^2E_r=0 \\[6pt] & \pdd{E_z}{r}+\frac{1}{r}\pd{E_z}{r}+\pdd{E_z}{z}+\left(\frac{\omega}{c_0}\right)^2E_z=0 \\[6pt] & \pdd{H_\phi}{r}+\frac{1}{r}\pd{H_\phi}{r}+\pdd{H_\phi}{z}-\frac{H_\phi}{r^2}+\left(\frac{\omega}{c_0}\right)^2H_\phi=0& \label{eq:solve_Hphi} \end{align}

これら三つの方程式は，式の形と境界条件が異なります．しかし，Maxwell 方程式と境界条件を一貫して扱えば，同じ固有周波数が得られます．このページでは，未知数が一つで済む式\eqref{eq:solve_Hphi}から出発します．

COMSOL 形式  式\eqref{eq:solve_Hphi}を式\eqref{eq:COMSOL_PDE}に合わせると，

\begin{align} \left(\frac{\omega}{c_0}\right)^2H_\phi+\acute{\nabla}\cdot\left[\acute{\nabla} H_\phi+\left(\frac{1}{r},\,0\right)H_\phi\right]=0 \end{align}

となります．式\eqref{eq:COMSOL_PDE}の係数は，

\begin{align} &e_a=\frac{1}{c_0^2}& &d_a=0& &c=-1& &\alpha=\left(-\frac{1}{r},\,0\right)& \nonumber\\[] &\gamma=\left(0,\,0\right)& &\beta=\left(0,\,0\right)& &a=0& &f=0& \label{ep:PDE:physics:coef} \end{align}

です．解は，固有値 $\lambda=\omega$，固有関数 $u=H_\phi$ です．

数値計算で安定な式

先に示した式\eqref{eq:solve_Hphi}は $1/r$ と $1/r^2$ を含む項があるため，軸上の $r=0$ では係数が形式的に発散します．連続的な解では，軸対称性から $H_\phi$ は $r$ に比例して 0 に近づくので，これらの特異項は互いに打ち消されます．しかし数値計算では，メッシュ点が $r=0$ に置かれると 0 で割る形になったり，$H_\phi$ やその微分のわずかな離散化誤差が $1/r^2$ の項によって大きく増幅されたりします．その結果，行列の条件数が悪くなる，収束が遅くなる，固有値や場分布に軸近傍の不自然な振動が混入する，といった不安定性が生じます．そこで，正則な関数 $v$ を用いて $H_\phi=v/r$ と変数変換を行います．変数変換の結果は，

\begin{align} & \pdiffB{2}{v}{r}-\frac{1}{r}\pdiffA{v}{r}+\pdiffB{2}{v}{z}+\left(\frac{\omega}{c_0}\right)^2 v=0 \label{eq:solve_v} \end{align}

です．

COMSOL 形式  式\eqref{eq:solve_v}を式\eqref{eq:COMSOL_PDE}に合わせると，

\begin{align} \left(\frac{\omega}{c_0}\right)^2v+\acute{\nabla}\cdot\acute{\nabla}v+\left(-\frac{1}{r},\,0\right)\cdot\acute{\nabla}v=0 \end{align}

となります．式\eqref{eq:COMSOL_PDE}の係数は，

\begin{align} &e_a=\frac{1}{c_0^2}& &d_a=0& &c=-1& &\alpha=\left(0,\,0\right)& \nonumber\\[] &\gamma=\left(0,\,0\right)& &\beta=\left(-\frac{1}{r},\,0\right)& &a=0& &f=0& \label{ep:PDE:v:coef} \end{align}

です．解は，固有値 $\lambda=\omega$，固有関数 $u=v$ です．

周期的境界条件の式

式\eqref{eq:solve_v}を直接使う場合，$v$ は波数 $k_z$ の Bloch 周期条件，または準周期条件を満たします．これは普通の意味での周期関数ではなく，1 周期進むと位相因子だけが掛かるという条件です．

\begin{align} v(r,z+L)=v(r,z)\exp(i k_z L) \end{align}

COMSOL の Periodic Condition では未知関数そのものを連続にしたいので，位相因子を方程式側へ移し，

\begin{align} v(r,z)=w(r,z)\exp(i k_z z) \label{eq:bloch_transform} \end{align}

と変数変換します．このように，Bloch 周期条件を満たす関数を位相因子とセル内の未知関数に分けることを，ここでは Bloch-Floquet 分解と呼びます．文献や分野によっては Bloch 変換，またはフロケ変換とも呼ばれます．

$L$ は周期長で，図1の場合は cell_D です．$k_z=\phi/L$ で，$\phi$ は 1 周期あたりの位相進みです．図3の赤色の $v$ は，波数 $k_z$ の Bloch 周期条件を満たす関数です．図では 1 セル分だけを描いていますが，これは普通の周期関数の一部ではなく，隣のセルへ行くと位相が $\exp(i k_z L)$ だけ進む関数です．

一方，緑色の $w$ は COMSOL の 1 セル領域 $0\le z\le L$ で解く未知関数です．計算では，$w$ が全空間で最初から周期関数として与えられていると考えるよりも，セル両端に

\begin{align} w(r,L)=w(r,0) \end{align}

という値の Continuity 条件を課す，と考える方が実装に近いです．ただし，COMSOL の Periodic Condition は値だけでなく，Coefficient Form PDE の flux も対応させます．この式での flux は，

\begin{align} \Gamma=-c\acute{\nabla}w-\alpha w+\gamma =\acute{\nabla}w+\left(0,\,ik_z\right)w \end{align}

です．したがって，$z=0$ と $z=L$ の対応境界では，値の連続に加えて，法線方向 flux の連続も課されます．周期境界が $z$ 方向の両端面であることを使うと，これは

\begin{align} \left(\pdiffA{w}{z}+ik_z w\right)_{z=L} = \left(\pdiffA{w}{z}+ik_z w\right)_{z=0} \end{align}

に対応します．さらに $w(r,L)=w(r,0)$ も課されるので，係数が両端で同じなら，結果として $\partial w/\partial z$ も両端で一致します．ただし，FEM の境界条件としては微分値を直接指定するというより，flux の連続を弱形式で課す，と理解するのが正確です．

式\eqref{eq:solve_v}に式\eqref{eq:bloch_transform}を代入します．

\begin{align} \pdiffA{v}{z}&=\left(\pdiffA{w}{z}+i k_z w\right)\exp(i k_z z)\\[6pt] \pdiffB{2}{v}{z}&=\left(\pdiffB{2}{w}{z}+2i k_z\pdiffA{w}{z}-k_z^2 w\right)\exp(i k_z z) \end{align}

したがって，$w$ についての方程式は，

\begin{align} \left(\frac{\omega}{c_0}\right)^2 w+\pdiffB{2}{w}{r}+\pdiffB{2}{w}{z}+2ik_z\pdiffA{w}{z}-k_z^2w-\frac{1}{r}\pdiffA{w}{r}=0 \label{eq:solve_w} \end{align}

となります．これが解くべき偏微分方程式です．

COMSOL 形式  式\eqref{eq:solve_w}を式\eqref{eq:COMSOL_PDE}に合わせると，

\begin{align} \left(\frac{\omega}{c_0}\right)^2w+\acute{\nabla}\cdot\left[\acute{\nabla}w+\left(0,\,ik_z\right)w\right]+\left(-\frac{1}{r},\,ik_z\right)\cdot\acute{\nabla}w-k_z^2w=0 \end{align}

となります．式\eqref{eq:COMSOL_PDE}の係数は，

\begin{align} &e_a=\frac{1}{c_0^2}& &d_a=0& &c=-1& &\alpha=\left(0,\,-ik_z\right)& \nonumber\\[] &\gamma=\left(0,\,0\right)& &\beta=\left(-\frac{1}{r},\,ik_z\right)& &a=-k_z^2& &f=0& \label{ep:PDE:w:coef} \end{align}

です．解は，固有値 $\lambda=\omega$，固有関数 $u=w$ です．ここで $\alpha_z$ と $\beta_z$ を分けて入れることが重要です．この分け方により，強形式の $+2ik_z\,\partial w/\partial z$ が再現されます．

解くべき方程式が定まったので，次に境界条件を考えます．最初に COMSOL で境界条件を設定する式を示し，その後，各変数での境界条件を確認します．

完全導体

ここでは，加速管の完全導体の境界条件を示します．図1の赤色の部分が完全導体の境界条件を設定する部分です．この完全導体の境界条件は，電磁場の物理的な境界条件から導かれます．

COMSOL 形式

この完全導体の境界条件は COMSOL では，Flux/Source を使います．これは，

\begin{align} -\vm{n}\cdot\left(-c\acute{\nabla}u-\alpha u+\gamma\right)=g-qu \label{eq:BC:COMSOL} \end{align}

です．$u$ は固有関数，$\vm{n}$ は境界の法線単位ベクトルでドメインの外向きです．係数 ($c,\,\alpha,\,\gamma$) は，偏微分方程式を決めた時の係数です．残りの変数 $(g,\,q)$ はユーザーが設定します．

物理的な偏微分方程式

境界条件のうち，完全導体 (金属) の表面部分は，

\begin{align} &\vm{H}\cdot\vm{n}=0& &\vm{E}\parallel\vm{n}& \label{eq:EM_BC} \end{align}

です．$\vm{n}$ は境界表面の法線ベクトルです．ここで問題としている軸対称 TM$_{0\,n}$ モードでは，$\vm{H}\cdot\vm{n}=0$ は自動的に満たされます．問題は，$\vm{E}\parallel\vm{n}$ です．ここでは $H_\phi$ を未知数にしているので，

\begin{align} \nabla\times\vm{H}=\varepsilon\pdiffA{\vm{E}}{t} \end{align}

から計算します．これから，電場は

\begin{align} \vm{E}=\frac{i}{\omega\varepsilon}\left[-\pdiffA{H_\phi}{z}\hat{r}+\frac{1}{r}\pdiffA{(rH_\phi)}{r}\hat{z}\right] \end{align}

です．これが金属境界 $\vm{n}=(n_r,\,n_z)$ と平行になるためには，$\vm{n}_\perp=(-n_z,\,n_r)$ と電場の内積がゼロになる必要があります．これは，

\begin{align} 0=\vm{E}\cdot\vm{n}_\perp&=\pdiffA{H_\phi}{z}n_z+\frac{1}{r}\pdiffA{(rH_\phi)}{r}n_r \nonumber \\[6pt] &=\left[\pdiffA{H_\phi}{r}+\frac{H_\phi}{r},\,\pdiffA{H_\phi}{z}\right]\cdot\vm{n} \label{eq:BC:physics} \end{align}

です．この最後の式が完全導体の境界条件になります．

COMSOL 形式  式\eqref{eq:BC:physics}を式\eqref{eq:BC:COMSOL}に合わせます．係数 ($c,\,\alpha,\,\gamma$) は，式\eqref{ep:PDE:physics:coef}の通りです．すると，式\eqref{eq:BC:COMSOL}は

\begin{align} -\vm{n}\cdot\left(\pdiffA{H_\phi}{r}+\frac{H_\phi}{r},\,\pdiffA{H_\phi}{z}\right)=g-qu \end{align}

です．これを式\eqref{eq:BC:physics}と一致させるためには，COMSOL では $g=0$，$q=0$ と設定します．

数値計算で安定な式

ここでは，前節「数値計算で安定な式」の場合の境界条件を考えます．これは，変数変換 $H_\phi=v/r$ により得られる偏微分方程式です．完全導体の境界条件は，式\eqref{eq:BC:physics}に変換 $H_\phi=v/r$ を適用すると得られます．この結果は

\begin{align} \left[\pdiffA{v}{r},\,\pdiffA{v}{z}\right]\cdot\vm{n}=0 \label{eq:BC:v} \end{align}

です．

COMSOL 形式  式\eqref{eq:BC:v}を式\eqref{eq:BC:COMSOL}に合わせます．係数 ($c,\,\alpha,\,\gamma$) は，式\eqref{ep:PDE:v:coef}の通りです．すると，式\eqref{eq:BC:COMSOL}は

\begin{align} -\vm{n}\cdot\left(\pdiffA{v}{r},\,\pdiffA{v}{z}\right)=g-qu \end{align}

です．これを式\eqref{eq:BC:v}と一致させるためには，COMSOL では $g=0$，$q=0$ と設定します．

周期的境界条件の式

最後に，前節「周期的境界条件の式」の場合の境界条件を考えます．これは，変数変換 $v=w\exp(i k_z z)$ により得られる偏微分方程式です．完全導体の境界条件は，式\eqref{eq:BC:v}に変換 $v=w\exp(i k_z z)$ を適用すると得られます．この結果は

\begin{align} \left[\pdiffA{w}{r},\,\pdiffA{w}{z}+ik_z w\right]\cdot\vm{n}=0 \label{eq:BC:w} \end{align}

です．

COMSOL 形式  式\eqref{eq:BC:w}を式\eqref{eq:BC:COMSOL}に合わせます．係数 ($c,\,\alpha,\,\gamma$) は，式\eqref{ep:PDE:w:coef}の通りです．すると，式\eqref{eq:BC:COMSOL}は

\begin{align} -\vm{n}\cdot\left(\pdiffA{w}{r},\,\pdiffA{w}{z}+ik_z w\right)=g-qw \end{align}

です．これを式\eqref{eq:BC:w}と一致させるためには，COMSOL では $g=0$，$q=0$ と設定します．$ik_z w$ の項は，式\eqref{ep:PDE:w:coef}の $\alpha_z=-ik_z$ によって Flux/Source の左辺へすでに含まれています．

対称軸 ($r=0$)

軸対称構造の対称軸 ($r=0$) 付近では，電場の $r$ 成分はゼロでなければなりません．また，$\nabla\times\vm{H}=\varepsilon\pdiffA{\vm{E}}{t}$ と時間依存性 $\exp(-i\omega t)$ から，軸近傍では，

\begin{align} H_\phi=-\frac{i\omega\varepsilon r}{2}E_z \end{align}

となります．これから $H_\phi$ は $r$ に比例して 0 に近づくことがわかります．したがって，対称軸での $H_\phi$ の境界条件は，Dirichlet 条件となり，$H_\phi=0$ です．

このため，$v=rH_\phi$ と $w=v\exp(-i k_z z)$ は $r=0$ で 0 となり，かつ滑らかである必要があります．COMSOL では，軸上に Dirichlet 条件 $u=0$ を設定します．滑らかな解では，軸近傍で $\pdiffA{v}{r}=0$ も満たされます．

周期的境界条件

Bloch 位相は式\eqref{eq:solve_w}に含めたので，COMSOL の Periodic Condition では Type of periodicity を Continuity にします．このとき，周期境界で COMSOL が課す条件は，

\begin{align} u_{dst}&=u_{src} \\ -\vm{n}_{dst}\cdot\left(-c\acute{\nabla}u-\alpha u+\gamma\right)_{dst}&=\vm{n}_{src}\cdot\left(-c\acute{\nabla}u-\alpha u+\gamma\right)_{src} \end{align}

です．ここでの $u$ は COMSOL の未知関数で，本文中の $w$ です．式\eqref{ep:PDE:w:coef}の係数を使うと，flux は

\begin{align} \Gamma=\acute{\nabla}w+\left(0,\,ik_z\right)w \end{align}

となります．したがって，周期境界では $w$ の値の一致に加えて，対応する法線方向 flux も一致します．$z=0$ と $z=L$ の両端面では，実質的に

\begin{align} w(r,L)&=w(r,0) \\[6pt] \left(\pdiffA{w}{z}+ik_z w\right)_{z=L} &= \left(\pdiffA{w}{z}+ik_z w\right)_{z=0} \end{align}

に対応します．物理場の位相差 $\exp(i k_z L)$ は，プロット式や後処理で $\exp(i k_z z)$ を掛けることで復元します．

180 度モードだけを検算する場合は，Bloch 変換を使わず，物理場 $v$ に Antiperiodicity を直接課す方法も使えます．ただし，この方法は診断用です．Bloch 変換版の式\eqref{eq:solve_w}と同時に Antiperiodicity を使うと，位相を二重に入れることになります．

ここでは，関数 $w$ から電場と磁場を計算する式を整理します．式\eqref{eq:bloch_transform}と $H_\phi=v/r$ から，時間因子を除いた複素振幅の磁場は，

\begin{align} \vm{H}(r,z)&=H_\phi(r,z)\hat{\phi} \\[6pt] H_\phi(r,z)&=\frac{w(r,z)}{r}\exp(i k_z z) \label{eq:EM:Hphi_from_w}\\[6pt] H_r(r,z)&=0,\qquad H_z(r,z)=0 \end{align}

です．時間依存性を含める場合は，式\eqref{eq:phase_convention}と同じように $\exp(-i\omega t)$ を掛けます．実数の場を描く場合は，最後に実部を取ります．磁束密度が必要な場合は $\vm{B}=\mu\vm{H}$ です．

電場は Maxwell 方程式

\begin{align} \nabla\times\vm{H}=\varepsilon\pdiffA{\vm{E}}{t} \end{align}

から求めます．時間依存性を $\exp(-i\omega t)$ としているので，複素振幅では

\begin{align} \nabla\times\vm{H}=-i\omega\varepsilon\vm{E}, \qquad \vm{E}=\frac{i}{\omega\varepsilon}\nabla\times\vm{H} \end{align}

です．軸対称 TM モードでは $E_\phi=0$ であり，

\begin{align} \vm{E} &=\frac{i}{\omega\varepsilon} \left[ -\pdiffA{H_\phi}{z}\hat{r} +\frac{1}{r}\pdiffA{(rH_\phi)}{r}\hat{z} \right] \label{eq:EM:E_from_Hphi} \end{align}

となります．式\eqref{eq:EM:Hphi_from_w}を代入すると，$w$ から電場成分は次のように計算できます．

\begin{align} E_r(r,z) &= \frac{k_z w-i\pdiffA{w}{z}}{\omega\varepsilon r} \exp(i k_z z) \label{eq:EM:Er_from_w}\\[6pt] E_z(r,z) &= \frac{i}{\omega\varepsilon r} \pdiffA{w}{r} \exp(i k_z z) \label{eq:EM:Ez_from_w}\\[6pt] E_\phi(r,z)&=0 \end{align}

これらが，$w$ から求める電磁場の複素振幅です．表示用の実場が必要な場合は，これらの複素振幅の実部を取ります．

最終的に COMSOL に入力する主な設定を表 2 にまとめます．Dependent Variable 名を w とした場合です (本文中の $w$ に対応)．

COMSOL Coefficient Form PDE の設定

項目	設定
位相規約	$\exp(i k_z z)\exp(-i\omega t)$．CST Studio の phase advance と同じ向きにする．
kz	cell_phase/cell_D
ea	1/(c_const)^2
da, c, f	da=0，c=-1，f=0
alpha	alpha_r=0，alpha_z=-i*kz
beta	beta_r=-1/r，beta_z=i*kz
a	-kz^2
PEC 壁	Flux/Source で g=0，q=0．
周期境界	Periodic Condition の Type of periodicity は Continuity．値と Coefficient Form PDE の flux が対応する．
対称軸	Dirichlet 条件で u=0．
電磁場	磁場   $H_\phi$: (w/r*exp(i*kz*z))， 電場   $E_r$: (kz*w-i*d(w, z))/(real(lambda)*epsilon_0*r)*exp(i*kz*z)     $Ez$: i/(real(lambda)*epsilon_0*r)*d(w, r)*exp(i*kz*z)
周波数	real(lambda)/(2*pi)

上の設定で，CST Studio の phase advance と COMSOL の計算結果がほぼ一致することを確認しました．表 3 は，同じセル形状で得られた比較結果です．

CST Studio と COMSOL の固有周波数比較

Phase Shift [deg/cell]	CST Studio [GHz]	COMSOL [GHz]	$|df/f|$
0	2.835561430	2.835587571	9.219E-06
30	2.837391814	2.837417856	9.178E-06
60	2.842387559	2.842413328	9.066E-06
90	2.849200159	2.849225554	8.913E-06
120	2.855999278	2.856024292	8.758E-06
150	2.860968053	2.860992786	8.645E-06
180	2.862784953	2.862809583	8.603E-06

180 度では，Bloch 変換を使わない Antiperiodicity の直接計算でも同じ周波数が得られました．これは，セル形状，メッシュ，基準 PDE，位相規約が大きく外れていないことの確認になります．

参考資料

更新履歴