有限要素法の使用上のヒント

はじめに

このチュートリアルでは,有限要素法をNDSolveおよび関連する関数(NDEigensystem等)と一緒に使用したときに起こり得る問題と,そのような問題を回避する方法について述べる.

有限要素法パッケージをロードする.

非線形の定常偏微分方程式の解の進捗を監視する

偏微分方程式の解を得るのに時間がかかることがある.これは,特に大規模な非線形偏微分方程式を扱う場合によくあることである.非線形偏微分方程式の解を得る過程を監視するために,EvaluationMonitorStepMonitorを使うことができる.次の例では,非線形偏微分方程式を解き,関数評価,ステップ,関数行列式評価の数を測定する.この例では特にPauseコマンドを加え,より詳しく調べられるように解を得る過程をゆっくりにする.偏微分方程式の解を得る過程それ自体がゆっくりである場合には,もちろんPauseを使用する必要はない.

非線形偏微分方程式と境界条件を設定する.
関数評価,ステップ,関数行列式評価の数を監視しながら,方程式を解く.

非定常の偏微分方程式の時間積分の進捗を監視する

偏微分方程式の時間積分には,時間がかかることがある.特に,境界条件あるいは係数が時間依存である,あるいは使う偏微分方程式の領域が三次元であるときに,これは真である.非定常の偏微分方程式の積分の進捗を監視する簡単な方法として,EvaluationMonitorを使う方法がある.

EvaluationMonitorは,他の場合と同じように,NDSolve内で使える.以下は,時間依存のディリクレ条件を持つ,波動方程式の例である.

1つの焦点に円形の穴を持つ,楕円形の領域を指定する.
領域のメッシュを作成し,可視化する.
時点 までアクティブな円形の穴上の非定常のディリクレ境界条件を指定する.
進捗を監視しながら,偏微分方程式を時間積分する.

非定常偏微分方程式のステップサイズプロット

非定常偏微分方程式のモデルは,時の経過とともに変化するシミュレーションを扱う.非定常偏微分方程式を解く際には,NDSolveや関連する関数は,与えられた初期条件で始め,それ以降の偏微分方程式の解を構築し始める.NDSolveは,適応時間積分法を使って,従属変数の進化を計算する.適応的ということは,取る時間ステップのサイズがある時点においてどのくらい動態が存在するかによって変わることを意味する.

他の有限要素ソフトウェアプログラムの中には,ユーザ自身が時間積分が問題なく行われていることをチェックしなければならないものもある.このようなことは,NDSolveを使う際には不要である. 理由が何であれ,確度や精度の目標(陰的に請求される場合もある)に従って,時間積分が行えない場合や,ステップサイズが小さすぎてうまくいかない場合には,NDSolveは,自動的にメッセージを発する.しかしながら,偏微分方程式の動態がある時点でどのようになっているかを理解するためには,NDSolveが取ったステップサイズを見ることが役に立つ.

シミュレーション過程中に,ステップサイズがどのように変化するかを可視化するためには,StepDataPlotを使う.これは,NDSolveの解で使われたステップサイズを対数スケールで返す.StepDataPlotを使うためには,Numerical Differential Equation Solvingパッケージをロードしなければならない.

ステップサイズプロットの可視化は,例を使うと分かりやすい.

矩形領域を定義する.

この例では,HeatTransferPDEComponentのデフォルトのパラメータ(.熱源 で与えられる.)を使って,時間依存の熱伝導方程式を設定する.

時間依存の熱伝導方程式を設定する.
時間と空間,および係数 に依存する関数 を指定する.
初期条件とシミュレーションの終了時間 を設定する.

方程式は,熱的に絶縁された境界条件を使って,すべての境界上で解く.これは,境界条件が与えられていないデフォルト設定である.

方程式を解く.
StepDataPlotを使ってステップサイズプロットを可視化する.

軸上でステップサイズの値が増えるステップサイズプロットは,時間依存ソルバがより大きな時間ステップを取ることを示す.他の有限要素法ソフトウェアの中には,ステップサイズの逆数をプロットすることによって,そのステップサイズプロットを表示するものもある.

メモリを多く必要とする偏微分方程式を解く

以下に挙げる例の中には,大量のRAMメモリを必要とするものがある.

有限要素法の特定の問題を解く際にNDSolveが必要とするメモリの量は,方程式の数と操作が行われている空間次元に比例する.3Dの結合系については,もとになっているハードウェアが提供するメモリの量が障壁となることがある.以下のセクションでは,この潜在的な障壁問題を軽減する,簡単に使えるオプションについて説明する.

まず$HistoryLengthを0に設定する.これで,Outに保存された前の結果が保存されないようにすることができる.

$HistoryLengthをゼロに設定する.

有限要素分析において,メモリに関して問題となる2つの事柄がある.離散化の作成と,方程式系の解(LinearSolve)である.NDSolveは,離散化と解法の際に必要となるメモリ量に影響するメッシュ生成とLinearSolveのステップにオプションを提供する.概算で,内部で生成されたメッシュの座標数に空間次数を掛けたものが,生成される離散化された系の行列の自由度(行と列の数)である.

例として,線形弾性から応力演算子を考えてみよう.これについては,固体力学のモノグラフに詳しい説明がある.応力演算子は,3D空間における3つの結合偏微分方程式の系である.このセクションで重要なのは,大規模な方程式系を作成することであるので,偏微分方程式自体についてこれ以上詳しく述べない.

応力演算子,ヤング(Young)係数,ポアソン率を定義する.

一連の測定を通して,さまざまなオプションが解法時間と解を求めるのに必要なメモリ量に影響するかを示す.このために,ヘルパー関数を実装する.関数は,入力として任意数のNDSolveのオプションを取り,偏微分方程式を解くのに必要な時間とメモリ量(MB)で返す.これは,NDSolveが内部で設定する系の行列の大きさを推し量るためにも役に立つ.系の行列の大きさは,その自由度である.つまり,系の行列内の行と列の数である.

オプション付きのNDSolveを呼び出すヘルパー関数を作成する.測定された時間,必要なメモリ,解を求めるための自由度が表示される.

NDSolveのデフォルト設定は,正確な解がかなり効率的に求まるように設定される.さまざまな種類のコンピュータでこの例が適用できるように,四面体要素が使われる.これらの結果をより小さな行列で使う.

四面体要素がメッシュになるように強制する設定で,NDSolveを呼び出す.

デフォルト設定では,有限要素の離散化を行うNDSolveは,ソルバメソッドとしての"Pardiso" と一緒にLinearSolveを呼び出す.この方法は,非常に正確だが,使用できるメモリを多く必要とする直接のソルバを使う.

変形された構造を可視化する.

通常デフォルトの"Pardiso"ソルバ以外のソルバだとあまり効率的ではない.

デフォルトの"Pardiso"の代りに"Multifrontal"ソルバを使って,NDSolveを呼び出す.

"Pardiso"ソルバのオプションで,特定の問題を解くのに必要なメモリの量を減らすために使えるのは,アウトオブコア機能である.これは,LU成分が他のものと一緒にディスクに保存され,RAMメモリを消費しないことを意味する.このことによって,より大きな問題を解くことができる.1点留意しなければならないのは,解の過程の一部はディスクに保存されるが,システム行列がコンピュータのRAMに完全にフィットする必要があるということである.

NDSolveを細分化されたメッシュと"Pardiso"ソルバで呼び出し,データをアウトオブコア("OOC")に保存する.

使用されたメッシュがより細かいために,自由度の数もかなり増えたことに注意する.留意しなければならないのは,十分なRAMメモリが使用可能である場合には,アウトオブコアの機能の速度がゆっくりになるということである.このオプションは,NDSolveが使用できるメモリがなくなった場合にのみ使われるべきであり,勘や推測で使うものではない.

メモリ消費量を減らすもう1つの方法に,NDSolveでより少ない数の要素を持つメッシュを生成する方法がある.

NDSolveにより小さいメッシュを使うように指示する.

これで解を求めるのに必要な時間とメモリが軽減されたが,確度が下がる.正しい解を求めるのにメッシュにいくつの要素が必要であるかは,問題によって異なる.最初は粗いメッシュを使い,解が大きく異ならなくなるまでだんだん細かくしていくのがよい.

断面図で,計算された結果の違いを可視化する.
0から0.01までの範囲にスケールされ直されたImage3Dとして,ノード距離の差分を可視化する.

解を計算するのに必要なメモリ量をさらに軽減する方法に,一次で正確なメッシュだけを使う方法がある.

NDSolveに一次で正確なメッシュを使うように指示する.
断面図で,計算された結果の違いを可視化する.
0から0.01までの範囲にスケールされ直されたImage3Dとして,ノード距離の差分を可視化する.

これで計算時間とメモリ量が軽減されたが,確度が下がった.一次メッシュを使うと,より効率的に解けるより小さな系の行列が返される.一次メッシュを使うと,確度がかなり下がることがある.

大規模な問題を解くのに使えるもう一つの方法に,2つのステップを使って解く方法がある.まず粗く,一次の可能性があるメッシュで行う.その後,そこで求まった解を反復ソルバの最初のベクトルとして使う.

粗いメッシュを作成する.
粗いメッシュ上で問題を解く.
細かいメッシュを作成する.
粗い解を細かいメッシュ上で評価する.
最初のベクトルの長さを調べる.これは自由度に等しい.
反復ソルバと細かいメッシュ上の最初のベクトルを使って,偏微分方程式を解く.
計算された結果の違いを切断面で可視化する.
ノードの距離の違いを,0から0.01までの範囲でスケールし直したImage3Dとして可視化する.

ここで述べたオプションは組み合わせることもできる.NDSolveに与えることができるオプションの詳しい説明については,有限要素のチュートリアル「有限要素のためのNDSolveオプション」を参照されたい.

繰り返し評価する際のメモリ消費

偏微分方程式を繰り返して何度も解くと,連続する評価でメモリ消費量が増えていくことがある.これを示すために,パラメトリック関数を作成し,それを繰り返し評価する.

パラメトリック関数を作成する.

パラメトリック関数を複数回評価し,10回評価するごとにメモリ消費量をレポートにする.

パラメトリック関数を何度も評価し,メモリ消費量を表示する.

メモリ消費量が増加するのは,キャッシュデータの結果である.データは,評価速度を速めるためにキャッシュされる.メモリに制約がある場合には,キャッシュを随時消去することによって,保存されたデータを削除する.キャッシュの消去は,ClearSystemCacheを使って行う.

導関数を計算する際のキャッシュもメモリを消費する.

導関数を計算する際にキャッシュをオフにする.
システムキャッシュを消去してからパラメトリック関数を何度も評価し,メモリ使用量を表示する.

解領域の外挿

このセクションは, 「有限要素のためのNDSolveオプション」に動かされた.

非連続係数のオーバーシュートまたはアンダーシュート

前のバージョンのWolfram言語では,領域内で不連続係数の補間を直接行う機能は限られていた. バージョン14ではこの問題が修正された.以下の例を考えてみよう.

2つの領域を持つ要素メッシュを作成して可視化する.
メッシュを離散化して可視化する.

それぞれの領域に2つの異なる値が割り当てられる(例えば,白色 0 と 灰色 1)とき,領域の境界部分で不連続が起る.

メッシュ上で不連続関数を評価する.

EvaluateOnElementMeshDiscontinuousInterpolatingFunctionを返すことに注意する.これはメッシュに複数の材料マーカーが与えられたためである.材料マーカーについての詳細は,「要素メッシュの生成」チュートリアルの「マーカー」セクションを参照されたい.

メッシュのメッシュ要素マーカーを調べる.
不連続補間関数を可視化する.

不連続補間関数は,境界領域で鮮明である.また最小値と最大値が与えられた関数の範囲内にある.

補間関数の極値を調べる.

何らかの理由でそのほうが望ましい場合には,前のバージョンの動作に切り換えることもできる."DiscontinuousInterpolation" Falseのオプションが指定されると, "EvaluateOnElementMesh"InterpolatingFunction iを返す.

要素メッシュ上での評価からInterpolatingFunctionを作成する.
補間関数の極値を調べる.

補間関数の極値は,指定された範囲であるを超える値である.

補間関数を可視化する.
InterpolatingFunctionDiscontinuousInterpolatingFunctionの差を可視化する.

DiscontinuousInterpolatingFunctionがインターフェースで厳密にどのように動作するかは制御できる.

で不連続値をプロットする:
のインターフェースでの不連続補間関数を評価する:

インターフェースで値を選択することが望ましい場合もある.これは"MarkerPriority"を設定することによって可能である.DiscontinuousInterpolatingFunctionはメッシュの要素マーカーを見て,メッシュのマーカーによって不連続インターフェースでの評価の優先度を決める.

メッシュのメッシュ要素マーカーを調べる:

これは,インターフェースに点がある場合に,マーカー0に割り当てられた値がマーカー1に割り当てられた値よりも優先されることを意味する.

DiscontinuousInterpolatingFunctionのマーカー優先度を調べる:

次に,DiscontinuousInterpolatingFunctionを再生成する.ただし,今回はマーカーの優先度が異なっている.

DiscontinuousInterpolatingFunctionにカスタムのマーカー優先度を使う:
不連続関数をプロットする:
インターフェースを評価する:

プロットは同じであるが,インターフェースにおいて,マーカー1に関連付けられたデフォルト値がマーカー0での値に優先するようになったことに注意する.

対流が卓越する方程式の安定化

対流と拡散についての方程式を考える.

係数についての詳細は,InitializePDECoefficientsを参照されたい.

対流成分 が大きくなる状況では,偏微分方程式の解が下の例のように不安定になることがある.

対流・拡散の方程式の演算子を指定する.
境界条件を指定する.
領域とメッシュを指定する.
大きな対流成分を持つ対流・拡散の方程式を解く.

メッセージは,解が不安定である可能性を示唆する.ペクレ(Péclet)数は,無次元数であり,物理数量についての対流と拡散の間の比を表現する.ペクレ数は以下で計算される.

は対流成分, は要素の最小直径, は拡散係数である.ペクレ数がメッシュ次数よりも大きい場合には,解の安定性が問題になることがあり,メッセージが発せられる.

大きな対流成分を持つ対流・拡散の方程式の解を可視化する.

結果はめちゃくちゃである.このような場合に行えることがいくつかある.しかしその前にまず,ペクレ数がどのように計算されたかを調べる.

メッシュ要素の最小直径を計算する.
対流と拡散の係数,および最小要素直径からペクレ数を計算する.

対流の項 あるいは拡散の項 が空間座標の関数である場合には,この項は単一の検定座標で評価される.この検定座標は,InitializePDECoefficientsで説明した"VerificationData"オプションを指定することによって変更することが可能である.

指定した境界条件によって不安定性が生じることもある.

指定された境界条件.

における境界条件は,その方向への の流れを抑制する従属変数 の値を指定する.これらは,「流出」を阻止しているのである.

のベクトル場を調べる:

におけるディリクレ条件を強制しないことによっても不安定性を避けることができる.

安定した境界条件.

境界条件は必ずしも選ばれないということもある.解の安定性を改善する1つの方法として,いわゆる人工的な拡散を加えることがある.この人工的な拡散は,偏微分方程式を安定化させるために拡散の項に加えられる.

必要となる人工的な拡散の量を推定する.
人工的な拡散を加えることによって,安定した結果が返されることを確かめる.

新たに計算されたペクレ数がメッシュ次数よりも小さいことに注意する.

人工的な拡散(AD)を使って,対流成分が大きい偏微分方程式を解く.
大きな対流成分を持つ対流・拡散の方程式の解を可視化する.
解の断面図を可視化する.

これですでにかなり改善された.追加の拡散である程度解が不鮮明になっている.これをさらに改善するために,人工的な拡散が の流線に従うようにさせることができる.しかしこの方法は,流れ場 に発散がない場合にしかうまくいかない.この方法は,流線風上ペトロフ・ガラーキン(PetrovGalerkin)法(SUPG)と呼ばれる.

流れ場 に発散がないかどうかをテストする.
調整パラメータ を最小のメッシュ直径に設定する.
流線風上の人工的な拡散(SUPG)を使って,対流成分が大きい偏微分方程式を解く.
SUPGを使って計算した解を可視化する.
純粋に人工的な拡散と流線風上化された拡散を比べる.

SUPGを使った解は,純粋に人工的な拡散を使って計算した解よりも において拡散が小さい.

これらの方法の欠点は,偏微分方程式が変更されることである.したがって,これが常によい方法であるとは言えない.代りに,メッシュを変更することができる.メッシュ生成についての詳細は,「要素メッシュの生成」 を参照されたい.

細分化された境界でメッシュを生成,可視化する.
事前定義されたメッシュ上で,大きな対流成分を持つ対流・拡散の方程式を解く.
大きな対流成分を持つ対流・拡散の方程式の解を可視化する.
解の断面図を可視化する.

この解は,上の例に比べて改善されたが,計算時間は余計にかかった.実際問題として,メッシュの右上の境界を細分化するだけで十分である.この部分だけを細分化すると,計算時間が短縮される.

次のステップとして,ToGradedMeshを利用して,右側と上部に向かうにつれてぼかされている細かいメッシュを構築する.

1Dグレーデッドメッシュを構築する.
グレーデッドメッシュから領域の積を作成する.
グレーデッドメッシュを可視化する.
メッシュ要素の最小半径を計算する.
大きな対流成分を持つ対流・拡散の方程式を事前定義されたメッシュ上で解く.
大きな対流成分を持つ対流・拡散の方程式の解を可視化する.
臨界終点での解の断面図を可視化する.

対流が卓越する時間依存の方程式の安定化

以下のセクションでは,数値的な不安定が,対流が卓越する時間依存の方程式でどのように起るのかを示す.このセクションでは,拡散の欠如によって起る数値的な不安定を最小化するための人工的な拡散の使用を示す.

例として一次の結合音波方程式を考える.

方程式は2つの従属変数 および を持つ.系には拡散項ががないので,偏微分方程式の解が不安定になることがある.

系を解くために,係数 および を指定してから,波動方程式の偏微分方程式演算子を設定する.

波動方程式を設定する.

境界条件として,角周波数 の正弦波音源が左の境界に置かれる.

境界条件を設定する.

これは時間依存の偏微分方程式であるので,初期条件を指定しなければならない.

初期条件を指定する.

波動方程式を解くためには,正確な数値解を得るために,音波の波長 が十分に細かいメッシュ によって解かれる必要がある.ここでは格子の大きさ を音波長の少なくとも60分の1より小さいものと設定する.

メッシュセルの大きさの上限を計算する.
メッシュの必要条件を満たすために,格子の大きさを と設定する.
方程式系を解く.
従属変数 の解を時点 で可視化する.

解のさざ波のような揺れは,数値的な不安定から来るものである.安定性を改善するために,偏微分方程式系に人工的な拡散を加えることができる.

偏微分方程式に従属変数が1つだけ含まれる場合にのみ,ペクレ数を使って必要な人工的な拡散を推定することができる.ここでは2つの未知の変数を持つ方程系であるので,手作業で拡散係数 を選ばなければならない.

拡散係数 を設定する.

は,人工的な拡散の影響を抑制する,手動で選ばれたパラメータである.偏微分方程式は人工的な拡散を加えた際に変わるので,ユーザは,忠実性を減らすことなく解を安定化する ϵ の値を1から10までの間で注意深く選ぶ必要がある.

各方程式に人工的な拡散の項を加えた方程式系を解く.
時点 における 従属変数 の解を可視化する.

人工的な拡散によって安定性が改善されたことに注目する.次に,数値結果を解析結果と比べる.

解析解を設定する.
不安定な解および安定化した解を解析解と比べる.

偏微分方程式が変更されるので,人工的な拡散を加えることによって解の勾配の変化は滑らかになる.このため,波の前部でより大きなエラーが起るのである.

波の前部を調べる.

偏微分方程式の係数形式

この例では,有限要素法における偏微分方程式の係数形式の重要性が提示されている.単一の従属変数の偏微分方程式で係数形式は以下の方程式(1)で与えられる.

係数についての詳細は,InitializePDECoefficientsを参照されたい.

偏微分方程式はで定義される.ここで は,求める解の従属変数である.係数 はスカラー, はベクトル,× 行列である.

この例の重要な点は,この方程式(およびその結合バージョン)のみが有限要素法で解ける方程式であるということである.この点を強調するために,より簡単な方程式(2)を考慮する.

方程式とその係数を設定する.偏微分方程式の空間定義域としてを使い,で停止するように時間積分を設定する.拡散係数 を1に設定し,これが内部で単位行列に変換される.減衰係数 として, の値には1,それ以外では2の値を持つ区分的連続関数を使う.

係数 を指定し,上の方程式(3)を設定する.
有限要素法をソルバとして強制して,方程式を解く.
で解をプロットする.

方程式(4)を以下のように再公式化することもできる.

方程式(5)を設定する.
方程式を解く.
で解をプロットする.
解1と解2の差分をプロットする.

明らかに解1と解2は異なる.何が起こっているのかを調べるために,NDSolveのデータ構造から実際に解析された偏微分方程式の係数を抽出する小さな補助関数を書く.有限要素法の内部構造についての詳細は,有限要素プログラミングのチュートリアルを参照されたい.

NDSolveから状態データを抽出するための補助関数を指定する.
偏微分方程式の係数を抽出するための補助関数を指定する.
方程式(6)について偏微分方程式の係数を調べる.

方程式(7)の係数は期待される通りのものである.係数を抽出する代りに,微分方程式がどのように解析されるかを調べる.

解析された微分法的式を再構築する.

ここまでは何も変な部分はない.

方程式(8)について偏微分方程式の係数を調べる.
解析された微分方程式を再構築する.

方程式(9)の係数は予想と全く異なるものである.係数 が拡散係数 に入っている.なぜだろうか.これを理解する鍵は,有限要素法がこのようなタイプの方程式しか解けないということにある.

の前には係数がないことに注意されたい.のタイプの方程式を使うためには, に再設定し, を調整してによって引き起こされる導関数を取り除く必要がある.以下はその例である.

方程式(10)について同じ方法を取ると,以下が得られる.

これは事実上明示的に を指定することと同じである.数学的には,不連続性のため,係数 はディラクデルタ関数のようなものであり,時間積分ではどのみち削除されてしまう.

偏微分方程式の係数を調べる.

不連続の拡散係数を持つ,時間依存の偏微分方程式については,係数が拡散項から削除されるような形で方程式を並べる方法が最も望ましい.

形式的な偏微分方程式

形式的な偏微分方程式とは,Inactiveコンポーネントを持つ方程式である.Inactiveの偏微分方程式がタイプセットされる場合,Inactiveの偏微分方程式の成分は灰色で表示される.このセクションでは,いつ偏微分方程式が非アクティブでなければならないかについて説明する.

Inactiveの偏微分方程式を使う主な理由は,偏微分方程式係数時期尚早に評価されることを防ぐことにある.例えば,勾配の評価された発散はNeumannValueが何を意味するかに影響を与える.

勾配の発散を評価する.
勾配のInactive発散を評価する.

非アクティブなバージョンでは,方程式に関連して,可能なNeumannValueが何を意味するかについて厳密に制御することができる.アクティブな形式ではこれはできない.NDSolveは引数を未評価のままにすることができないので,NDSolveおよび関連の関数は,自動的にこの操作を行うことはできない.つまり,NDSolveAttributesHoldFirstHoldAllのいずれも持たない.

Inactive形式の方程式が必要となるもう一つの例は,偏微分方程式係数が非対称の場合である.例は,平面応力演算子である.ヤング率とポアソン比については,鉄鋼の値が使われる.

平面応力演算子を定義する.

長さ10メートル,高さ0.7メートルで左側で固定された梁の最初の固有値と固有関数を計算する.

梁の固有関数を計算する.
固有値から共振周波数を計算する.

比較のために,オイラー・ベルヌーイの梁理論を利用した最初の共振周波数を解析的に計算する.

結果における差分を比べる.差分の約半分は,異なる理論的なアプローチによるものである.

これは予測される通りである.

次に,Activeバージョンの偏微分方程式を調べ,結果の共振周波数を計算する.

上の結果は間違っている.何が起こったのだろうか.簡単に言えば,Inactiveではないバージョンは非対称な係数には使えないのである.4つの係数行列のもとの非アクティブな偏微分方程式を見ると,対角ではない2つは対称ではない.

言い換えれば,NDEigensystem(およびその他の NDSolve等の関連関数)が意図された形式で偏微分方程式を構文解析する前に,評価されるのである.

評価された偏微分方程式から非対称な係数を回復する方法はない.このため,Inactiveが導入されたのである.非対称な係数に対しては,必ずInactiveバージョンの偏微分方程式を使わなければならない.

構文解析された係数を調べるために,NDEigensystemの時間依存方程式の入力にはNDSolve`ProcessEquationsが使える.NDEigensystemへの入力が演算子形式である場合には,NDSolve`ProcessEquationsの時間依存方程式に変換しなければ,入力に使えない.任意の初期条件と任意の時間積分区間を指定しなければならないことに注意する必要がある.

非アクティブな形式の方程式が必要である2つ目の例は,以下の偏微分方程式のように,拡散係数内に従属変数の微分を持つ非線形の偏微分方程式である.

.

拡散係数 は,従属変数 を含むので非線形である.

方程式を設定する.
非線形方程式を解く.

ここで方程式が非アクティブな形式ではない場合には,NDSolveはこの偏微分方程式を解くことができない.

アクティブな形式の非線形方程式を解くと,メッセージが発せられる.

現在NDSolveは,係数のいずれかに一次よりも高次の微分を持つ従属変数の定常偏微分方程式を解くことはできない.これは,以下のように説明される.アクティブな演算子は係数形式 に合致する必要がある.まず,拡散係数を として設定する.

拡散係数 を係数形式に挿入すると,以下の結果が返される.

残りの項 は,係数形式で も設定することによって相殺される.

しかし, を評価すると, 係数に一次より高次の導関数微分が返される.これは,現在NDSolveの有限要素法では解くことができない.

有限要素法におけるInactiveの目的は,方程式(11)で与えられる偏微分方程式の係数形式を形式的な偏微分方程式として設定することである.方程式(12)の数学演算子のすべてが非アクティブであってもよい.言い換えれば,この機能はDivGradを使って表現できるので,Inactiveの偏微分方程式にはInactive[Div]Inactive[Grad]の成分を含むことができるが,Inactive[D]等の非アクティブなDの成分あるいはInactive[Derivative]等の非アクティブなDerivativeの成分は含むことはできない.

ときには,偏微分方程式がNDSolveによってどのように解析されるかを調べると便利なことがある.これは,GetInactivePDEを使って行える.

すでに解析された非アクティブな偏微分方程式を調べる.

より詳しい情報は,偏微分方程式の係数形式のセクションを参照されたい.

ノイマン値と形式的な偏微分方程式

このセクションでは,形式的な偏微分方程式NeumannValueとの相互関係について説明する.定常的な場合に偏微分方程式は以下のように表わすことができる.

係数についての詳細は,InitializePDECoefficientsを参照されたい.

それに対して,一般化されたノイマン境界値は,値 を指定する.値 は,境界の部分上の外向きの法線での流束を指示する.

ノイマン境界値の微分については,「偏微分方程式と境界条件」で説明されている.方程式(13)における係数 は,NeumannValueの使用によって直接指定されているわけではないことに注意することが重要である.係数 は,方程式(14)で指定されている. NeumannValueを使用すると,が境界上の の値に置き換えられる.

DivGradの使用のみでは,方程式(15)を完全に指定するには不十分である.以下の指定を考える.

これは対応するノイマン値である.

として与えられ, が単位行列, がゼロベクトルである場合を考える.発散の評価の結果,以下が返されることに注意する.

ここで の発散(この例では)である. 方程式(16)と(17)は等しい.

上の方程式が等しいことを確かめる.

しかし,評価で問題が起こる.方程式は等しいが,ノイマン指定は等しくない.方程式(18)のノイマン値は,陰的に以下に変更されている.

このことについてさらに説明するために,「有限要素プログラミング」のチュートリアルで説明されるように,偏微分方程式を解析し,有限要素データから偏微分方程式データを抽出するヘルパー関数を作成する.

NDSolveの有限要素データを生成し,偏微分方程式の係数データを抽出するヘルパー関数を作成する.
偏微分方程式を設定し,偏微分方程式の係数データを抽出する.

PDECoefficientDataは,解析された偏微分方程式係数を含むデータ構造である.

拡散係数は, を持つ行列である.
対流係数 は, である.
反応係数は, の値を持つ.
保存型の対流係数はない.

評価中に,係数 は, の項に分割される.微分方程式を入力する正しい構造は,Inactivateによって,あるいは直接Inactiveを使用することによって評価される係数である.

偏微分方程式を非アクティブにする.
偏微分方程式の係数データを抽出する.
拡散係数は, である.
対流係数は, である.
対流係数に存在しない.
反応係数は存在しない.

γ の項がある場合には,以下の非アクティブ偏微分方程式を使うこともある.

項を持つ偏微分方程式を非アクティブ化する.
偏微分方程式の係数データを抽出する.
負荷導関数係数は, である.

以下の例では,実用的な面から見た,ノイマン値についてアクティブと非アクティブの偏微分方程式の意味の違いを説明する.この例では,バイアスポテンシャル について,領域 に限定される拡散を示す.

穴がある長方形領域を指定する.
バイアスポテンシャルとその勾配 を指定する.

ディリクレ境界条件は,領域の上下()をそれぞれポテンシャル1と0で設定する.

でディリクレ境界条件を指定する.

残りの2つの辺()には,自然な境界条件(ノイマンゼロ)が適用される.つまり,この2つの辺では,流束はゼロである.

モデル化する偏微分方程式は,以下で与えられる.

におけるその対応するノイマン値は以下の通りである.

非アクティブ偏微分方程式を指定する.
非アクティブ偏微分方程式を解く.
偏微分方程式係数を調べる.

まず,におけるノイマン値の充足を調べる.

非アクティブ演算子の正規導関数 を,陰的に指定されたノイマンゼロ値の辺上でプロットする.正規導関数の方向は である.

におけるゼロ流束の境界条件(2つの辺上でゼロ流束)で,もとの項がない場合には, におけるボックスへの流束は,おけるボックスからの流束と同じはずである.

非アクティブ演算子の正規導関数 をディリクレ条件を持つ辺上でプロットする.正規導関数の方向は である.

ボックスへの流束は,ボックスからの流束と等しい.比較のために,アクティブ偏微分方程式を解く.

アクティブ偏微分方程式を解く.
偏微分方程式係数を調べる.

係数は,対流と反応の係数として解析される.

対流と反応の係数を調べる.

アクティブの場合,計算される境界条件は である.

アクティブ演算子の正規導関数 をディリクレ条件を持つ辺上でプロットする.正規導関数の方向は である.

ときには,偏微分方程式がNDSolveによってどのように解析されるかを調べると便利なことがある.これは,GetInactivePDEを使って行える.

すでに解析された非アクティブな偏微分方程式を調べる.

より詳しい情報は,偏微分方程式の係数形式のセクションを参照されたい.

偏微分方程式係数の効率的な評価

偏微分方程式の係数は定数一定でなくてもよい.この例では,拡散係数 はシミュレーションの領域における位置によって変わる.

指定された係数を評価する順序を理解することが重要である.順序がシミュレーションのパフォーマンスに影響を及ぼすことがあるからである.まず,内部境界のあるメッシュを生成する.

内部境界と細かい解像度を持つメッシュを作成する.
内部境界を可視化する:

メッシュを設定したら,非特定の拡散係数 を持つ偏微分方程式を設定する.

ポワソン偏微分方程式を設定する.
境界条件を設定する.

拡散係数 は,内部領域では拡散係数の値が となり,その領域の外側では拡散係数が となるように設定される.

領域の部分によって値が異なる 係数を設定する.
境界条件を持つ方程式を解き,必要な時間を測る.

シミュレーションは予期されるように行われた.拡散係数 をより詳しく見てみると,If文の分岐が評価されなかったことが分かる.

係数を調べる.

If文がその分岐引数を評価しないのは,If文の属性が評価を妨げる設定になっているからである.

Ifの属性を調べる.

拡散係数 におけるIf文の分岐引数が評価されないということは,呼出しのたびに変数が検索されなければならないことを意味する.その結果,If文をコンパイルすることができず,標準評価器が使われる.この場合,この評価器を使った方が時間がかかる.

係数関数のすべての引数が評価されるように関数を作成するには,Withを利用することができる.

値が評価された 係数を作成する.

これでIf文のすべての引数が評価されていることが見て取れる.

新しい 係数でポワソン偏微分方程式を設定する.
スピードが上がったことを確かめる.

ここで使ったような小さい例ではこれはあまり大きな問題ではないように見えるが,大きなシミュレーションにおいては,完全に評価された係数は考慮に入れるべきことである.

If文に内容を入れるためにWithを使うよりも簡単に使えるもう一つの方法に,Piecewiseを使う方法がある. Piecewiseはそれに関連付けられた属性を持たないので,その引数を評価する.

評価された値で 係数を作成する.

Piecewise文ではすべての引数が評価されることに注意する.

新しい 係数でポワソン偏微分方程式を設定する.
スピードが上がったことを確かめる.

スピード面では,Piecewise構造はWith/If構造と同じくらい効率的である.

HeatTransferPDEComponentのような偏微分方程式成分を使って偏微分方程式が設定されている場合には,効率的な係数の生成が自動的に行われることに注意する.

材料のパラメータ parsにおける架空の材料境界で設定する.
同じパラメータ pars 内の材料パラメータ値を設定する.
偏微分方程式モデルを設定する.

記号のパラメータ値が実際の数値に置き換えられたことに注意する.

係数がSetDelayedの後ろに隠れている場合にもほとんど同じ問題が起る.

分岐の関数を使う係数のSetDelayed定義.
係数の定義を調べる.
同じ係数が評価されるように定義する.
その係数の定義を調べる.

2つ目の場合には,偏微分方程式モデル内の係数として使われるときに,検定は記号的に評価され,より効率的である.

偏微分方程式の係数内の補間関数の効率的な評価

偏微分方程式の係数も補間関数になり得る.偏微分方程式が解かれるメッシュと同じものを補間関数が利用する場合には,補間関数を非常に効率的に評価することが可能である.

細分化されたメッシュの領域を作成する.
補間関数を作成する.
補間関数を係数として,偏微分方程式を解き,時間を測る.
補間関数のメッシュが計算メッシュと同じものではないことを確認する.
計算メッシュと同じメッシュ上に補間関数を作成する.
補間関数を係数として,偏微分方程式を解き,時間を測る.
補間関数のメッシュが計算メッシュと同じであることを確かめる.

補間関数のメッシュが,偏微分方程式に使われたメッシュと同じである場合には,計算に必要な時間は短くなる.

最も効率的な方法は,可能であれば,式を直接与えることである.

式を偏微分方程式の係数として使う.

NeumannValueと境界導関数の関係

以下の例では,NeumannValueと境界導関数の関係を説明し,パラメータが変更しやすいように設定する.まず,領域を指定する.

領域を指定する.

これが解くべき偏微分方程式である.

偏微分方程式は,従属変数 を1つ持つ.これは一次元系であるので,拡散係数 は1×1行列であり,これは数として指定できる.

偏微分方程式の演算子を設定する.

ディリクレ境界条件は,従属変数 について値 を指定する.

ノイマン値は以下の形式で指定される.

次に,境界条件を設定する.

方程式の境界条件を指定する.
方程式系を解く.
DSolveValueを使って解析解を計算する.

方程式のDSolveValueの形式化において,一般化されたノイマン境界方程式がどのように明示的に与えられなければならないか,そしてそれらの係数がどのように偏微分方程式の係数に一致しているかに注意しなければならない.

解析結果を数値結果と比べるために,この2つの差分のプロットを作成する.

解析解と数値解の差分をプロットする.

従属変数名の順序

偏微分方程式の結合系を分析する場合に,NDSolveで実装された有限要素法が予期しない結果を返すということがある.このセクションでは,どのような場合にこれが起こるのかを説明し,この問題を避ける方法を提示する. 最後に,さらに詳しく知りたいという読者のために,より詳細な説明を提供する.

まず,適切な境界条件を持つ方程式系と領域を設定する.

領域を指定する.
3つの方程式の系を指定する.
その境界条件を指定する.
偏微分方程式系を指定する.

この例における従属変数は,uvw である.系を解き,従属変数 w を可視化する.

系を解き,従属変数 w を可視化する.

次の例では,従属変数の名前を w から p に変更し,方程式系をもう一度解く.

従属変数 wp で置き換えた系を解き,従属変数 p を可視化する.

その結果は予期しないものとなった.何が起こったかを詳しく説明する前に,予期される動作が必ず起こるようにする2つの方法について説明する.

1つ目は,NDSolve"DependentVariable"メソッドオプションを通して従属変数を指定する方法である.

従属変数を指定し,p の解をプロットする.

"DependentVariable"メソッドオプションを指定することは,従属変数が指定の次数で離散化されるように強制する.このことについては,下で詳しく説明する.

この例の場合には,NDSolveは,偏微分方程式の離散化の際に従属変数が使うべき補間次数を指定するための有限要素法のオプションを与える.補間次数の指定は,陰的に方程式と従属変数を同じ次数に強制する.

従属変数の補間次数を指定し,p の解をプロットする.

通常,それぞれの従属変数に同じ補間次数を指定しても問題ない.補間次数の低減も問題の根幹にあるものではない.この場合には,方程式は流体方程式(ストークスの流れ)であり,安定化は,圧力変数のInterpolation次数が速度 u および v の補間次数よりも1小さい場合には強制することができる.

根底にある問題が何であるかを知るために,両方の場合について構築された離散系の行列を表示する.

オプションが与えられていない場合の,構築された系の行列を表示する.
順序のオプションが与えられた場合の,構築された系の行列を表示する.

この動作をさらに調べるために,NDSolveの状態オブジェクトを作成する.

NDSolveの状態オブジェクトを作成する.

FiniteElementDataFEMMethodDataの従属変数の順序は{p,u,v}であることに注意する.

状態オブジェクトの従属変数を調べる.

PDECoefficientDataの拡散と対流の係数を調べると,結合偏微分方程式の係数の構造が分かる.

結合偏微分方程式の係数を表示する.
結合偏微分方程式の対流係数を表示する.

有限要素の補間次数法のオプションが設定されると,係数構造は,従属変数の正規順序に調整される.

NDSolveの状態オブジェクトを作成する.
結合偏微分方程式の拡散係数を表示する.
結合偏微分方程式の対流係数を表示する.

実際,方程式の入力の順序を変更して,従属変数の正規順序にマッチさせると,予期される通りの結果が得られる.

方程式の順序を変えて,結合偏微分方程式を解く.

それでは,なぜこの順序変更が自動的に行われないのであろうか.

先天的な問題として,指定された境界条件は,係数の対角成分にしか適用できないということがある.これは,境界条件が入力されたときに,特定の境界条件がどのように偏微分方程式の特定の方程式に関連付けられるかということの情報が失われるためである.

解析された境界条件の内部データを調べる.

ここでは,それぞれのサブリストは2つの部分に分かれている.最初の部分は,離散化の際に境界条件の2つ目の部分が適用される行と列を示す.

もう一つの方法は,ディリクレ境界条件をそれぞれが関係している方程式に関連付ける方法である.

ディリクレ境界条件が関連付けられている結合偏微分方程式系を作成する.
結合偏微分方程式を解き,解の一部を可視化する.

最後に,最適の方法として,補間次数を指定する方法がある.

従属変数の補間次数を指定し,p の解をプロットする.

解の検証

NDSolveや関連関数から得た微分方程式の解を検証する方法として,手元の微分方程式に解を再挿入してみるということがある.

境界条件を設定する.
微分方程式を設定し,その解を求める.
解を微分方程式に再挿入して,解を検証する.
検証を可視化する.

上の例を見ると,微分方程式に解を再挿入する方法が,微分方程式の解を検証するのに好ましいと思われるかも知れないが,必ずしもそうとは言えない.

微分方程式の再挿入した解を可視化する.

この結果は,解が前の例ほど正確ではないことを示唆する.実際にそうだろうか.

前の例では,解が微分方程式に再挿入された場合に,InterpolatingFunctionのより高次な微分が計算されなければならないという陰的な仮定がされた.しかし,微分回復の過程に誤差が伴わないとは限らず, の尺度までうまく計算できない.二次よりも高次のメッシュを使うと状況が改善され,一次のメッシュを使うとより悪化する.試してみよう.

NDSolveValueに対して1のメッシュ次数を指定するメッシュオプション.

理想的には,微分方程式の解は,解の検証中にさらなる数値誤差をもたらさずにチェックすることが望ましい.そのためには,作られた解のメソッドが便利である.微分方程式の解析解が存在していたなら,計算された改はこの解析解に対して確かめることができる.作られた解のアイディアは,任意解を選ぶためのものであり,選ばれた解が実際の解であるように,もとの微分方程式を修正するためのものである.

Sinを任意解として選ぶ.
選ばれた解をもとの微分方程式に挿入する.

選ばれた解析解をもとの微分方程式に挿入すると,新しい右辺の項 ができる.もとの方程式がこの新しい右辺 に等しくなるように設定されると,選択された解であるSinは,新しい微分方程式の解に構築される.

新しい方程式を作った解で解く.
数値解と解析解との間の誤差を可視化する.

この選ばれた解析解の場合の誤差は,解の形式NDSolveValueをもとの方程式に再挿入した場合のものよりもずっと小さくなる.この方法の利点は,解の検証中に新たな数値誤差がもたらされないということである.