PDEModelsの最良の方法

保存型対流と非保存型対流平滑化ステップ関数
スケールの違いで数値解が台無しになるメモリ負荷の高いPDEを解く
副次未知数量の不連続性単位
ディラックのデルタ関数

保存型対流と非保存型対流

一般的に，対流は保存型と非保存型でモデル化することができる．保存型はConservativeConvectionPDETermとして表され，をモデル化する．別の方法はConvectionPDETermを使うというものであり， をモデル化する． の位置が違うことに注目のこと．

スケールの違いで数値解が台無しになる

微分方程式モデルを数値的に解く場合，モデルのコンポーネントのスケールおよびその相対サイズ・絶対サイズに注意を払うことが大切である．例えば，形状と材料パラメータのスケールが大きく異なると解が台無しになることがある．また計算機イプシロンと比較した形状の絶対サイズが解の精度に悪影響を及ぼすこともある．一番の問題は，これが起こるときほぼ警告メッセージはなく，モデル作成者がこの影響を考慮しなければならないという点である．

この影響を調べるために，幅 ，閉込めポテンシャル の量子井戸を例に取る．ここではシュレーディンガー方程式を解き最初の数個の束縛固有状態を求める．

モデルのパラメータを定義する：

ここで は粒子の質量， は換算プランク定数である．パラメータは"SIBase"単位で指定されていることに注意する． したがって長さの単位は"Meters"になる．同時に井戸の長さはわずか "Angstroms"である．つまり領域の次元は "Meters"のスケールになる．また，換算プランク定数の"SIBase"単位は"Kilograms" "Meters"^2 /"Seconds"である．つまりこの換算プランク定数の長さの単位は "Meters"になる．

の大きさを調べる：

領域の長さの次元は，対応するモデルパラメータの長さの次元と数桁異なり，$MachineEpsilonよりも小さいことが分かる．

数値計算を続けよう．

FEMパッケージをロードする：

この問題の境界条件を構築する際に，物理学的な見地から気を付けなければならないのは，波動関数は井戸からの距離が無限大に近づくにつれてゼロに減衰する必要があるという点である．しかし，数値的に無限領域を定義することは不可能なので，十分大きい領域を発見的に選び，後で外側境界で となるようなディリクレ条件を定義する．ここの問題では，波動関数の減衰には井戸の長さの4倍の距離が適切と思われる．

問題の境界を定義する：
境界メッシュを定義する：
メッシュを定義し，そのワイヤフレームを表示する：
PDEモデルパラメータを定義する：
PDE演算子を定義する：

演算子の係数のスケールは10-20と10-39であり，計算機イプシロンよりずっと小さい．

計算機イプシロンを調べる：
ディリクレ境界条件を設定する：
固有値問題を解く：

固有値 vals を見ると，これらも計算機イプシロンより小さいことが分かる．これは解が正しくない可能性があることを示している．

ここで，得られた固有値と固有関数を見てみる．

固有値を"Millielectronvolts"に変換する：
確率密度をプロットする：

上のプロットでは，一見エネルギーの最初の状態の確率密度だけが量子井戸の基底状態で想定される関数形式に似ているように見える．それにもかかわらず，他の２つの確率密度は想定される形状にならない．

この時点で，シミュレーション結果を照合することは常に良い習慣だと言える．利用できる場合は解析結果を，または別の数値メソッドを使うかするとよい．ここでは解析解がある．解析解の詳細はここでは冗長になるので，以下のセクションで解析解を定義する．解析結果の代りに２つ目の数値結果が利用可能になってもプロセスは変わらない．

解析エネルギー固有値と数値エネルギー固有値の差分を取る：
解析値についての相対誤差を計算する：

解析固有値に対する誤差率に基づくと，得られた数値固有値のいずれも受け入れられないということは明らかである．特に2つ目と3つ目の固有値は，1つ目の状態と比較して，対応する解析値からは遠く離れている．これはプロットでより簡単に可視化できる．

解析エネルギー固有値と数値エネルギー固有値をプロットする：

これらの問題を回避する方法として，よりよく整列するように形状および係数を適切にスケールするというものがある．これを行うために，井戸の長さを前に使った"Meters"の代りにスケール

に再定義する．

をスケール
に再定義する：
問題の境界を再定義する：
境界メッシュを再定義する：
メッシュを再定義し，そのワイヤフレームを表示する：

形状の次元が変わったので，PDE演算子と境界条件を再定義する必要がある．井戸の長さが再定義されたため，まずパラメータ"SchrodingerPotential"を更新しなければならない．さらにSchrodingerPDEComponentはデフォルトでは"SIBase"単位を使うため，長さの単位は"Meters"である．それに応じてPDE演算子を再定義するために，"ScaleUnits"パラメータを使う．この場合"ScaleUnits"->{"Meters"->"Angstroms"}を使う．これで物理的パラメータの長さの単位すべてが"Angstroms"に再スケールされ，必要なPDE演算子が提供される．

上と同じパラメータを使うが，"ScaleUnits"を加え"SchrodingerPotential"を更新する：
PDE演算子を再定義する：

ここで係数のスケールは1次と10次である．

ディリクレ境界条件を最定義する：

これを行ったら，通常通りに問題を解くことができる．

固有値問題を解く：

この場合，物理的パラメータの長さの単位を"Meters"から"Angstroms"に変えたので，得られた固有状態もそれに応じた単位を持つ．特にエネルギー固有値には，固有値を所望のエネルギー単位（この場合"Millielectronvolts"）に変換する前に"Meters"を"Angstroms"に変換することを除いて、対応する"SIBase"エネルギー単位を割り当てる．

固有値に適切な単位を与え，"Millielectronvolts"に変換する：
解析エネルギー固有値と数値エネルギー固有値の差分を取る：
解析値に対する相対誤差を計算する：
確率密度をプロットする：

形状の単位はAngstromsなので，波動関数の単位は通常のではなくになる．モデル作成者が波動関数で事後計算を行いたい場合は，このことを考慮に入れる必要がある．

また，確率密度と解析解を比較するためにも，単位の違いを考慮に入れる必要がある．解析波動関数の単位はである．そのため，解析確率密度出力および独立変数はだけスケールダウンされる．

数値確率密度の誤差をプロットする：

確率密度は量子井戸に囲まれた最初の3つの境界状態について，予想された関数形式を示している．さらに，エネルギーレベルの数値はその解析エネルギーレベルにほぼ一致する．

解析解

偶パリティ状態と奇パリティ状態の超越方程式を定義する：
超越方程式の根を計算する：

Rationalizeを使って，井戸の長さの10分の1の許容誤差を持つ方程式の厳密係数を得る．

解析エネルギー固有値をまとめて"Millielectronvolts"で示す：
固有値問題の解析解の関数を定義する：
最初のつの状態に対する解析確率密度をプロットする：

副次未知数量の不連続性

主要な未知変数から副次数量を計算しなければならない場合がある．このような副次数量は，主要な未知変数の第一導関数に比例することが多い．主要な未知変数とは異なり，副次変数は要素境界において不連続となり得る．

以下の例は，従属変数 の導関数を取る際に発生する不連続性を示している．また，副次数量 が要素数が増えるに従ってどのように正確になっていくかも示している．

解く方程式は係数 を持つ拡散項と空間で二次式的に変化するソース項である．方程式はから までの1D領域で解かれる．

ソース項は以下の方程式で表される：

ここで は定数， は独立変数である．

有限要素パッケージをロードする：
変数を定義する：
使用するパラメータを定義する：

3つの異なるメッシュを使う．最初のメッシュには2本の線の要素，2つのメッシュには4本線の要素，最後のメッシュには20本線の要素がある．

メッシュを生成する：
方程式を定義する：

行の両端で2つのDirichletConditionを指定する．左では値1，右では値0を指定する．

ディリクレ条件を定義する：
異なる領域に対するPDEを解く：

この問題には以下で与えられる解析解がある．解が4次になる点に注目する．つまり，補間関数で使える二次要素が完全には解を補間できないということである．

解析解を定義する：
要素境界を赤い点で示して異なる解を可視化する：

主要数量で見る故tができるように，解は要素数が多いほど正確になる．

解の補間関数は二次である．ここで解の導関数を計算すると，二次補間関数は実質的に線形になる．さらに線形要素が要素境界において不連続性を導入するというこてゃ有限要素の特徴である．しかしこれらの不連続性は次に示すようにメッシュを調整することで削減することができる．

解析式の導関数を計算する：
解析式の導関数を可視化する：
数値式の導関数を計算する：
2要素解の導関数を可視化する：
数値導関数と解析導関数の差を計算する：
不連続性の近くで数値導関数を評価する：
4要素解の導関数を可視化する：
数値導関数と解析導関数の差を計算する：
不連続性の近くで導関数を評価する：

各要素の交点で不連続性を見ることができるが，不連続性は最初のものよりも小さくなっている．

20要素解の導関数を可視化する：

20要素では不連続性は見られない．

数値導関数と解析導関数の差を計算する：
不連続性の近くで導関数を評価する：
解の導関数を近似するためにより多くの要素が使われるにつれ，どのように誤差が減衰していくかをプロットする：

ディラックのデルタ関数

正則化デルタ関数は，音響学等多数の物理分野で単極音源を実装するために使われる．

単極音源を使うためには，音源は になければならない．単極音源項 は以下のように書ける：

ここで は音源位置 におけるディラックのデルタ関数である．

しかし，ディラックのデルタ関数は離散化された空間領域では解けないので，数値シミュレーションでは問題が起こる．これは，ディラックのデルタ関数が音源位置 において特異だからである．従ってディラックのデルタ関数の近似が必要になる．ディラックのデルタ関数の近似の過程は正則化と呼ばれる．

利用できる正則化デルタ関数 にはいくつかある．例として以下がある：

ここで は正則化デルタ関数 のサポートを制御する正則化パラメータである．通常 はメッシュ間隔 に匹敵する大きさである．

を中心とする正則化デルタ関数 を設定する：
を中心とする正則化デルタ関数 を調べる．サポートの幅は に等しい：

を変更すると正則化デルタ関数の幅が変わるが，では空間積分は常に1である．

正則化デルタ関数上で積分する：

平滑化ステップ関数

いくつかの例では，熱流束 や表面温度 等の過渡パラメータに対する時間プロファイルを規定するために平滑化ステップ関数 を使う．平滑化ステップ関数以下のように定義される：

ここで関数 が到達できる最小値と最大値はそれぞれ と で表される．ステップの位置は で制御され，平滑化勾配は で制御される．

平滑化ステップ関数 を定義する：
平滑化ステップ関数 を可視化する：

メモリ負荷の高いPDEを解く

有限要素法のドキュメントには「メモリを多く必要とする偏微分方程式を解く」というセクションがあるので，PDEを解く上で知っておくと役に立つ．

単位

幾何学単位

幾何学単位は再スケールすることができる．これは「要素メッシュの生成」チュートリアルで説明してある．

材料単位

幾何学単位が材料単位と異なる場合，材料単位はスケールすることができる．

材料を設定し単位をスケールする：

内部的にはすべての材料データは"SIBase"単位に変換される．その結果長さのデフォルト単位は"Meters"になる．幾何学単位もメートルであれば，何も変更は必要ない．幾何学単位がメートルでなければ，PDEと材料特性を幾何学単位にスケールするか，機関学単位を"Meters"にスケールするかしなければならない．PDEと材料特性の単位をスケールするためには，パラメータ"ScaleUnits"を与えることができる．明示的に述べられていない場合，このノートブックの例ではデフォルトの"SIBase"単位を使う．

材料パラメータが数量として与えられていても，手順は同じである．

材料を設定して単位をスケールする：
パラメータから，スケールされた熱移動PDEを作成する：