有限要素法の使用上のヒント
はじめに
このチュートリアルでは，有限要素法をNDSolveおよび関連する関数（NDEigensystem等）と一緒に使用したときに起こり得る問題と，そのような問題を回避する方法について述べる．
Needs["NDSolve`FEM`"]非線形の定常偏微分方程式の解の進捗を監視する
偏微分方程式の解を得るのに時間がかかることがある．これは，特に大規模な非線形偏微分方程式を扱う場合によくあることである．非線形偏微分方程式の解を得る過程を監視するために，EvaluationMonitorとStepMonitorを使うことができる．次の例では，非線形偏微分方程式を解き，関数評価，ステップ，関数行列式評価の数を測定する．この例では特にPauseコマンドを加え，より詳しく調べられるように解を得る過程をゆっくりにする．偏微分方程式の解を得る過程それ自体がゆっくりである場合には，もちろんPauseを使用する必要はない．
eqn = -Inactive[Div][(1/Sqrt[1 + Grad[u[x, y], {x, y}] . Grad[u[x, y], {x, y}]])*
Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0;
bc = DirichletCondition[u[x, y] == Sin[2π * (x + y)], True];Module[{pauseTime = 0.25, e = 0, s = 0, j = 0},
Monitor[NDSolveValue[{eqn, bc}, u, {x, y}∈Disk[], Method -> {"FiniteElement", "PDESolveOptions" -> {"FindRootOptions" -> {
Jacobian -> {Automatic, EvaluationMonitor :> (j++;Pause[pauseTime])},
EvaluationMonitor :> (e++;Pause[pauseTime]),
StepMonitor :> (s++;Pause[pauseTime])
}}}
],
StringJoin[
"
Function Evaluations: ", ToString[e],
"
Steps: ", ToString[s],
"
Jacobian evaluations: ", ToString[j]]]]非定常の偏微分方程式の時間積分の進捗を監視する
偏微分方程式の時間積分には，時間がかかることがある．特に，境界条件あるいは係数が時間依存である，あるいは使う偏微分方程式の領域が三次元であるときに，これは真である．非定常の偏微分方程式の積分の進捗を監視する簡単な方法として，EvaluationMonitorを使う方法がある．
EvaluationMonitorは，他の場合と同じように，NDSolve内で使える．以下は，時間依存のディリクレ条件を持つ，波動方程式の例である．
a = 3; b = 2;r = 3 / 10;
f = Sqrt[a ^ 2 - b ^ 2];
Ω = RegionDifference[Disk[{0, 0}, {a, b}], Disk[{f, 0}, r]];mesh = ToElementMesh[Ω, MaxCellMeasure -> 0.05];
mesh["Wireframe"]
t1 = 7 / 10;
Γ = DirichletCondition[u[t, x, y] == If[t < t1, Sin[t Pi / t1] ^ 2, 0], (x - f) ^ 2 + y ^ 2 ≤ (11 / 10r) ^ 2];Monitor[NDSolveValue[{D[u[t, x, y], {t, 2}] == Laplacian[u[t, x, y], {x, y}], u[0, x, y] == 0, Derivative[1, 0, 0][u][0, x, y] == 0, Γ}, u, {x, y}∈mesh, {t, 0, 2}
, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor]非定常偏微分方程式のステップサイズプロット
非定常偏微分方程式のモデルは，時の経過とともに変化するシミュレーションを扱う．非定常偏微分方程式を解く際には，NDSolveや関連する関数は，与えられた初期条件で始め，それ以降の偏微分方程式の解を構築し始める．NDSolveは，適応時間積分法を使って，従属変数の進化を計算する．適応的ということは，取る時間ステップのサイズがある時点においてどのくらい動態が存在するかによって変わることを意味する．
他の有限要素ソフトウェアプログラムの中には，ユーザ自身が時間積分が問題なく行われていることをチェックしなければならないものもある．このようなことは，NDSolveを使う際には不要である． 理由が何であれ，確度や精度の目標（陰的に請求される場合もある）に従って，時間積分が行えない場合や，ステップサイズが小さすぎてうまくいかない場合には，NDSolveは，自動的にメッセージを発する．しかしながら，偏微分方程式の動態がある時点でどのようになっているかを理解するためには，NDSolveが取ったステップサイズを見ることが役に立つ．
シミュレーション過程中に，ステップサイズがどのように変化するかを可視化するためには，StepDataPlotを使う．これは，NDSolveの解で使われたステップサイズを対数スケールで返す．StepDataPlotを使うためには，Numerical Differential Equation Solvingパッケージをロードしなければならない．
Needs["DifferentialEquations`NDSolveUtilities`"]
Ω = Rectangle[];この例では，HeatTransferPDEComponentのデフォルトのパラメータ（
，
，
．熱源 
は 
で与えられる．）を使って，時間依存の熱伝導方程式を設定する．
eqn = HeatTransferPDEComponent[{u[t, x, y], t, {x, y}}, <|"HeatSource" -> g[t, {x, y}, 2]|>]
g[t_, {x_, y_}, τ_] :=
(2 t / τ ^ 2)Exp[-t ^ 2 / τ ^ 2] Sin[Pi x] ^ 2 Sin[Pi y] ^ 2 - (2 Pi ^ 2(1 - Exp[-t ^ 2 / τ ^ 2])(Cos[2 Pi x]Sin[Pi y] ^ 2 + Sin[Pi x] ^ 2 Cos[2 Pi y]));
ic = {u[0, x, y] == 0};
tend = 10;方程式は，熱的に絶縁された境界条件を使って，すべての境界上で解く．これは，境界条件が与えられていないデフォルト設定である．
ufun = NDSolveValue[{eqn == 0, ic}, u, {t, 0, tend}, {x, y}∈Ω]StepDataPlot[ufun, FrameLabel -> {"time [s]", "step size [s]"}, PlotLabel -> "Convergence Plot"]
軸上でステップサイズの値が増えるステップサイズプロットは，時間依存ソルバがより大きな時間ステップを取ることを示す．他の有限要素法ソフトウェアの中には，ステップサイズの逆数をプロットすることによって，そのステップサイズプロットを表示するものもある．
メモリを多く必要とする偏微分方程式を解く
以下に挙げる例の中には，大量のRAMメモリを必要とするものがある．
有限要素法の特定の問題を解く際にNDSolveが必要とするメモリの量は，方程式の数と操作が行われている空間次元に比例する．3Dの結合系については，もとになっているハードウェアが提供するメモリの量が障壁となることがある．以下のセクションでは，この潜在的な障壁問題を軽減する，簡単に使えるオプションについて説明する．
まず$HistoryLengthを0に設定する．これで，Outに保存された前の結果が保存されないようにすることができる．
$HistoryLength = 0有限要素分析において，メモリに関して問題となる2つの事柄がある．離散化の作成と，方程式系の解(LinearSolve)である．NDSolveは，離散化と解法の際に必要となるメモリ量に影響するメッシュ生成とLinearSolveのステップにオプションを提供する．概算で，内部で生成されたメッシュの座標数に空間次数を掛けたものが，生成される離散化された系の行列の自由度（行と列の数）である．
例として，線形弾性から応力演算子を考えてみよう．これについては，固体力学のモノグラフに詳しい説明がある．応力演算子は，3D空間における3つの結合偏微分方程式の系である．このセクションで重要なのは，大規模な方程式系を作成することであるので，偏微分方程式自体についてこれ以上詳しく述べない．
op = SolidMechanicsPDEComponent[{{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}}, <|"YoungModulus" -> 10, "PoissonRatio" -> 1 / 3|>]一連の測定を通して，さまざまなオプションが解法時間と解を求めるのに必要なメモリ量に影響するかを示す．このために，ヘルパー関数を実装する．関数は，入力として任意数のNDSolveのオプションを取り，偏微分方程式を解くのに必要な時間とメモリ量（MB）で返す．これは，NDSolveが内部で設定する系の行列の大きさを推し量るためにも役に立つ．系の行列の大きさは，その自由度である．つまり，系の行列内の行と列の数である．
performanceTest[opts___] := Block[{measure, ndsol}, measure = AbsoluteTiming[MaxMemoryUsed[ndsol = NDSolveValue[{op == {0, 0, NeumannValue[-1, z == 1]}, DirichletCondition[{u[x, y, z] == 0, v[x, y, z] == 0, w[x, y, z] == 0}, z == 0]}, {u, v, w}, {x, y, z}∈Cuboid[{0, 0, 0}, {1, 1, 1}], opts]
] / (1024. ^ 2)];Print[
"Time -> ", measure[[1]], "
Memory -> ", measure[[2]],
"
Degrees of Freedom -> ", Length[ndsol[[1]]["ElementMesh"]["Coordinates"]] * 3];
ndsol
]NDSolveのデフォルト設定は，正確な解がかなり効率的に求まるように設定される．さまざまな種類のコンピュータでこの例が適用できるように，四面体要素が使われる．これらの結果をより小さな行列で使う．
result = performanceTest[Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshElementType" -> TetrahedronElement}}}];
デフォルト設定では，有限要素の離散化を行うNDSolveは，ソルバメソッドとしての"Pardiso" と一緒にLinearSolveを呼び出す．この方法は，非常に正確だが，使用できるメモリを多く必要とする直接のソルバを使う．
{ufun, vfun, wfun} = result;
mesh = ufun["ElementMesh"];
Show[MeshRegion[ElementMeshDeformation[mesh, {ufun, vfun, wfun}, "ScalingFactor" -> 1]], mesh["Wireframe"]]通常デフォルトの"Pardiso"ソルバ以外のソルバだとあまり効率的ではない．
performanceTest[Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshElementType" -> TetrahedronElement}, "PDESolveOptions" -> {"LinearSolver" -> {Automatic, Method -> "Multifrontal"}}}}];
"Pardiso"ソルバのオプションで，特定の問題を解くのに必要なメモリの量を減らすために使えるのは，アウトオブコア機能である．これは，LU成分が他のものと一緒にディスクに保存され，RAMメモリを消費しないことを意味する．このことによって，より大きな問題を解くことができる．1点留意しなければならないのは，解の過程の一部はディスクに保存されるが，システム行列がコンピュータのRAMに完全にフィットする必要があるということである．
performanceTest[Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshElementType" -> TetrahedronElement, "MaxCellMeasure" -> 0.00005}, "PDESolveOptions" -> {"LinearSolver" -> {Automatic, Method -> {"Pardiso", "OOC" -> True}}}}}];
使用されたメッシュがより細かいために，自由度の数もかなり増えたことに注意する．留意しなければならないのは，十分なRAMメモリが使用可能である場合には，アウトオブコアの機能の速度がゆっくりになるということである．このオプションは，NDSolveが使用できるメモリがなくなった場合にのみ使われるべきであり，勘や推測で使うものではない．
メモリ消費量を減らすもう1つの方法に，NDSolveでより少ない数の要素を持つメッシュを生成する方法がある．
result2 = performanceTest[Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshElementType" -> TetrahedronElement, "MaxCellMeasure" -> 0.00025}}}];
これで解を求めるのに必要な時間とメモリが軽減されたが，確度が下がる．正しい解を求めるのにメッシュにいくつの要素が必要であるかは，問題によって異なる．最初は粗いメッシュを使い，解が大きく異ならなくなるまでだんだん細かくしていくのがよい．
{ufun2, vfun2, wfun2} = result2;
Plot[{ufun2[1 / 2, y, 1] - ufun[1 / 2, y, 1], vfun2[1 / 2, y, 1] - vfun[1 / 2, y, 1], wfun2[1 / 2, y, 1] - wfun[1 / 2, y, 1]}, {y, 0, 1}, PlotRange -> All]s = 0.025;
data = Table[Sqrt[(ufun2[i, j, k] - ufun[i, j, k]) ^ 2 + (vfun2[i, j, k] - vfun[i, j, k]) ^ 2 + (wfun2[i, j, k] - wfun[i, j, k]) ^ 2], {k, 0, 1, s}, {i, 0, 1, s}, {j, 0, 1, s}];
Image3D[Reverse[Rescale[data, {0, 0.01}]]]//Rasterize解を計算するのに必要なメモリ量をさらに軽減する方法に，一次で正確なメッシュだけを使う方法がある．
result3 = performanceTest[Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshElementType" -> TetrahedronElement, "MeshOrder" -> 1}}}];
{ufun3, vfun3, wfun3} = result3;
Plot[{ufun3[1 / 2, y, 1] - ufun[1 / 2, y, 1], vfun3[1 / 2, y, 1] - vfun[1 / 2, y, 1], wfun3[1 / 2, y, 1] - wfun[1 / 2, y, 1]}, {y, 0, 1}, PlotRange -> All]s = 0.025;
data = Table[Sqrt[(ufun3[i, j, k] - ufun[i, j, k]) ^ 2 + (vfun3[i, j, k] - vfun[i, j, k]) ^ 2 + (wfun3[i, j, k] - wfun[i, j, k]) ^ 2], {k, 0, 1, s}, {i, 0, 1, s}, {j, 0, 1, s}];
Image3D[Reverse[Rescale[data, {0, 0.01}]]]//Rasterizeこれで計算時間とメモリ量が軽減されたが，確度が下がった．一次メッシュを使うと，より効率的に解けるより小さな系の行列が返される．一次メッシュを使うと，確度がかなり下がることがある．
大規模な問題を解くのに使えるもう一つの方法に，2つのステップを使って解く方法がある．まず粗く，一次の可能性があるメッシュで行う．その後，そこで求まった解を反復ソルバの最初のベクトルとして使う．
meshCoarse = ToElementMesh[Cuboid[{0, 0, 0}, {1, 1, 1}], MaxCellMeasure -> 0.0005, "MeshOrder" -> 1]AbsoluteTiming[MaxMemoryUsed[ndsol = NDSolveValue[{op == {0, 0, NeumannValue[-1, z == 1]}, DirichletCondition[{u[x, y, z] == 0, v[x, y, z] == 0, w[x, y, z] == 0}, z == 0]}, {u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}∈meshCoarse]
] / (1024. ^ 2)]meshFine = ToElementMesh[Cuboid[{0, 0, 0}, {1, 1, 1}], MaxCellMeasure -> 0.000019, "MeshElementType" -> TetrahedronElement]AbsoluteTiming[MaxMemoryUsed[
sv = Flatten[#["ValuesOnGrid"]& /@ EvaluateOnElementMesh[{x, y, z}, ndsol, meshFine]];
] / 1024. ^ 2]Length[sv]AbsoluteTiming[MaxMemoryUsed[result4 = NDSolveValue[{op == {0, 0, NeumannValue[-1, z == 1]}, DirichletCondition[{u[x, y, z] == 0, v[x, y, z] == 0, w[x, y, z] == 0}, z == 0]}, {u, v, w}, {x, y, z}∈meshFine, Method -> {"PDEDiscretization" -> {"FiniteElement", "PDESolveOptions" -> {"LinearSolver" -> {Automatic, Method -> {"Krylov", Method -> "BiCGSTAB", Tolerance -> 10 ^ -3, "Preconditioner" -> "ILU0", "StartingVector" -> sv}}}}}]
] / (1024. ^ 2)]{ufun4, vfun4, wfun4} = result4;
Plot[{ufun4[1 / 2, y, 1] - ufun[1 / 2, y, 1], vfun4[1 / 2, y, 1] - vfun[1 / 2, y, 1], wfun4[1 / 2, y, 1] - wfun[1 / 2, y, 1]}, {y, 0, 1}, PlotRange -> All]s = 0.025;
data = Table[Sqrt[(ufun4[i, j, k] - ufun[i, j, k]) ^ 2 + (vfun4[i, j, k] - vfun[i, j, k]) ^ 2 + (wfun4[i, j, k] - wfun[i, j, k]) ^ 2], {k, 0, 1, s}, {i, 0, 1, s}, {j, 0, 1, s}];
Image3D[Reverse[Rescale[data, {0, 0.01}]]]//Rasterizeここで述べたオプションは組み合わせることもできる．NDSolveに与えることができるオプションの詳しい説明については，有限要素のチュートリアル「有限要素のためのNDSolveオプション」を参照されたい．
繰り返し評価する際のメモリ消費
偏微分方程式を繰り返して何度も解くと，連続する評価でメモリ消費量が増えていくことがある．これを示すために，パラメトリック関数を作成し，それを繰り返し評価する．
pfun = ParametricNDSolveValue[{-Laplacian[u[x, y], {x, y}] == f, DirichletCondition[u[x, y] == 0, True]}, u, {x, y}∈Disk[], f]パラメトリック関数を複数回評価し，10回評価するごとにメモリ消費量をレポートにする．
Do[
m = MemoryInUse[];
pfun[i];
If[Mod[i, 10] == 0,
Print["Memory usage increase: ", MemoryInUse[] - m]]
, {i, 1, 100}]
メモリ消費量が増加するのは，キャッシュデータの結果である．データは，評価速度を速めるためにキャッシュされる．メモリに制約がある場合には，キャッシュを随時消去することによって，保存されたデータを削除する．キャッシュの消去は，ClearSystemCacheを使って行う．
SetSystemOptions["CacheOptions" -> {"Derivative" -> {"Cache" -> False}}];Do[
ClearSystemCache[];
m = MemoryInUse[];
pfun[i];
If[Mod[i, 10] == 0,
Print["Memory usage increase: ", MemoryInUse[] - m]]
, {i, 1, 100}]
解領域の外挿
このセクションは， 「有限要素のためのNDSolveオプション」に動かされた．
非連続係数のオーバーシュートまたはアンダーシュート
前のバージョンのWolfram言語では，領域内で不連続係数の補間を直接行う機能は限られていた． バージョン14ではこの問題が修正された．以下の例を考えてみよう．
Ω = RegionDifference[Rectangle[], Disk[{1 / 2, 1 / 2}, 1 / 6]];mesh = ToElementMesh[Ω, "RegionHoles" -> None, "RegionMarker" -> {{{1 / 5, 1 / 5}, 1}}, "MaxCellMeasure" -> 0.001];
mesh["Wireframe"["MeshElementStyle" -> {FaceForm[Gray], FaceForm[White]}]]それぞれの領域に2つの異なる値が割り当てられる（例えば，白色 0 と 灰色 1）とき，領域の境界部分で不連続が起る．
function = If[ElementMarker == 1, 0, 1];
difun = EvaluateOnElementMesh[{x, y}, function, mesh]EvaluateOnElementMeshはDiscontinuousInterpolatingFunctionを返すことに注意する．これはメッシュに複数の材料マーカーが与えられたためである．材料マーカーについての詳細は，「要素メッシュの生成」チュートリアルの「マーカー」セクションを参照されたい．
mesh["MeshElementMarkerUnion"]Plot3D[difun[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All, Mesh -> All]不連続補間関数は，境界領域で鮮明である．また最小値と最大値が与えられた関数の範囲内にある．
MinMax[difun["ValuesOnGrid"]]何らかの理由でそのほうが望ましい場合には，前のバージョンの動作に切り換えることもできる．"DiscontinuousInterpolation" Falseのオプションが指定されると， "EvaluateOnElementMesh"はInterpolatingFunction iを返す．
ifun = EvaluateOnElementMesh[{x, y}, function, mesh, "DiscontinuousInterpolation" -> False]MinMax[ifun["ValuesOnGrid"]]
Plot3D[ifun[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All]Plot3D[{ifun[x, y] - difun[x, y]}, {x, 0, 1}, {y, 0, 1}, PlotRange -> All]DiscontinuousInterpolatingFunctionがインターフェースで厳密にどのように動作するかは制御できる．
Plot[difun[x, 1 / 2], {x, 0, 1}]
difun[1 / 2 + 1 / 6, 1 / 2]インターフェースで値を選択することが望ましい場合もある．これは"MarkerPriority"を設定することによって可能である．DiscontinuousInterpolatingFunctionはメッシュの要素マーカーを見て，メッシュのマーカーによって不連続インターフェースでの評価の優先度を決める．
mesh["MeshElementMarkerUnion"]これは，インターフェースに点がある場合に，マーカー0に割り当てられた値がマーカー1に割り当てられた値よりも優先されることを意味する．
difun["MarkerPriority"]次に，DiscontinuousInterpolatingFunctionを再生成する．ただし，今回はマーカーの優先度が異なっている．
difun2 = EvaluateOnElementMesh[{x, y}, function, mesh, "MarkerPriority" -> {1, 0}]Plot[difun2[x, 1 / 2], {x, 0, 1}]difun2[1 / 2 + 1 / 6, 1 / 2]プロットは同じであるが，インターフェースにおいて，マーカー1に関連付けられたデフォルト値がマーカー0での値に優先するようになったことに注意する．
対流が卓越する方程式の安定化
係数についての詳細は，InitializePDECoefficientsを参照されたい．
対流成分 
が大きくなる状況では，偏微分方程式の解が下の例のように不安定になることがある．
β = {Cos[π / 3], Sin[π / 3]};c = 10 ^ -3;
op = Div[-c Grad[u[x, y], {x, y}], {x, y}] + β.Grad[u[x, y], {x, y}]Γ = {DirichletCondition[u[x, y] == 1, x == 0 || y == 0],
DirichletCondition[u[x, y] == 0, x == 1 || y == 1]};Ω = Rectangle[{0, 0}, {1, 1}];
mesh = ToElementMesh[Ω]ufun = NDSolveValue[{op == 1, Γ}, u, {x, y}∈mesh]
メッセージは，解が不安定である可能性を示唆する．ペクレ(Péclet)数は，無次元数であり，物理数量についての対流と拡散の間の比を表現する．ペクレ数は以下で計算される．
は対流成分，
は要素の最小直径，
は拡散係数である．ペクレ数がメッシュ次数よりも大きい場合には，解の安定性が問題になることがあり，メッセージが発せられる．
Plot3D[ufun[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All, ColorFunction -> "TemperatureMap"]結果はめちゃくちゃである．このような場合に行えることがいくつかある．しかしその前にまず，ペクレ数がどのように計算されたかを調べる．
hMin = Power[Min[Join@@mesh["MeshElementMeasure"]], 1 / mesh["EmbeddingDimension"]]Norm[β] * hMin / (2 * Min[Abs[c]])対流の項 
あるいは拡散の項 
が空間座標の関数である場合には，この項は単一の検定座標で評価される．この検定座標は，InitializePDECoefficientsで説明した"VerificationData"オプションを指定することによって変更することが可能である．
Γ
と 
における境界条件は，その方向への 
の流れを抑制する従属変数 
の値を指定する．これらは，「流出」を阻止しているのである．
VectorPlot[β, {x, y}∈Ω]
と 
におけるディリクレ条件を強制しないことによっても不安定性を避けることができる．
{DirichletCondition[u[x, y] == 1, x == 0 || y == 0]};境界条件は必ずしも選ばれないということもある．解の安定性を改善する1つの方法として，いわゆる人工的な拡散を加えることがある．この人工的な拡散は，偏微分方程式を安定化させるために拡散の項に加えられる．
d = 1 / mesh["EmbeddingDimension"];
mo = mesh["MeshOrder"];
Subscript[c, art] = d / mo * hMin * Norm[β]Norm[β] * hMin / (2 * Min[c + Subscript[c, art]])新たに計算されたペクレ数がメッシュ次数よりも小さいことに注意する．
ufunAD = NDSolveValue[{op + Div[-Subscript[c, art] Grad[u[x, y], {x, y}], {x, y}] == 1, Γ}, u, {x, y}∈mesh];Plot3D[ufunAD[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All, ColorFunction -> "TemperatureMap"]Plot[ufunAD[x, 0.8], {x, 0, 1}]これですでにかなり改善された．追加の拡散である程度解が不鮮明になっている．これをさらに改善するために，人工的な拡散が 
の流線に従うようにさせることができる．しかしこの方法は，流れ場 
に発散がない場合にしかうまくいかない．この方法は，流線風上ペトロフ・ガラーキン(Petrov–Galerkin)法(SUPG)と呼ばれる．
Div[β, {x, y}] == 0
τ = hMin;ufunSUPG = NDSolveValue[{op + Div[-(τ TensorProduct[β, β]).Grad[u[x, y], {x, y}], {x, y}] == 1, Γ}, u, {x, y}∈mesh];Plot3D[ufunSUPG[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All, ColorFunction -> "TemperatureMap"]Plot[{ufunAD[x, 0.8], ufunSUPG[x, 0.8]}, {x, 0, 1}, PlotLegends -> "Expressions"]SUPGを使った解は，純粋に人工的な拡散を使って計算した解よりも 
において拡散が小さい．
これらの方法の欠点は，偏微分方程式が変更されることである．したがって，これが常によい方法であるとは言えない．代りに，メッシュを変更することができる．メッシュ生成についての詳細は，「要素メッシュの生成」 を参照されたい．
mesh = ToElementMesh[Rectangle[{0, 0}, {1, 1}], "MaxBoundaryCellMeasure" -> 0.005, "MeshElementType" -> TriangleElement];
mesh["Wireframe"]ufun = NDSolveValue[{op == 1, Γ}, u, {x, y}∈mesh];Plot3D[ufun[x, y], {x, 0, 1}, {y, 0, 1}, ColorFunction -> "TemperatureMap"]Plot[ufun[x, 0.8], {x, 0, 1}]この解は，上の例に比べて改善されたが，計算時間は余計にかかった．実際問題として，メッシュの右上の境界を細分化するだけで十分である．この部分だけを細分化すると，計算時間が短縮される．
次のステップとして，ToGradedMeshを利用して，右側と上部に向かうにつれてぼかされている細かいメッシュを構築する．
m1 = ToGradedMesh[Line[{{0}, {1}}], <|"Alignment" -> "Right", "ElementCount" -> 150|>]mesh = ElementMeshRegionProduct[{m1, m1}]mesh["Wireframe"]Out[34]= [image]hMin = Power[Min[Join@@mesh["MeshElementMeasure"]], 1 / mesh["EmbeddingDimension"]]ufunGraded = NDSolveValue[{op == 1, Γ}, u, {x, y}∈mesh];Plot3D[ufunGraded[x, y], {x, 0, 1}, {y, 0, 1}, ColorFunction -> "TemperatureMap"]Plot[{ufun[x, 0.8], ufunGraded[x, 0.8]}, {x, 0, 1}, PlotRange -> {{0.95, 1}, All}, PlotLegends -> "Expressions"]対流が卓越する時間依存の方程式の安定化
以下のセクションでは，数値的な不安定が，対流が卓越する時間依存の方程式でどのように起るのかを示す．このセクションでは，拡散の欠如によって起る数値的な不安定を最小化するための人工的な拡散の使用を示す．
方程式は2つの従属変数 
および 
を持つ．系には拡散項ががないので，偏微分方程式の解が不安定になることがある．
系を解くために，係数 
および 
を指定してから，波動方程式の偏微分方程式演算子を設定する．
ρ = 1.2;c = 343;
SystemOp = {D[p[t, x], {t, 1}] + ρ c^2 * D[v[t, x], {x, 1}], D[v[t, x], {t, 1}] + (1/ρ) * D[p[t, x], {x, 1}]};境界条件として，角周波数 
の正弦波音源が左の境界に置かれる．
ω = 5025;
Γ = {DirichletCondition[p[t, x] == Sin[ω t], x == 0],
DirichletCondition[v[t, x] == Sin[ω t] / (ρ c), x == 0]};これは時間依存の偏微分方程式であるので，初期条件を指定しなければならない．
ics = {p[0, x] == 0, v[0, x] == 0};波動方程式を解くためには，正確な数値解を得るために，音波の波長 
が十分に細かいメッシュ 
によって解かれる必要がある．ここでは格子の大きさ 
を音波長の少なくとも60分の1より小さいものと設定する．
N[1 / 60 * λ /. {λ -> 2π c / ω}]
h = 0.005;
mesh = ToElementMesh[Line[{{0}, {1}}], "MaxCellMeasure" -> h]pfun = NDSolveValue[{SystemOp == {0, 0}, Γ, ics}, p, {t, 0, 0.002}, {x}∈mesh]
Plot[pfun[0.0006, x], {x, 0, 0.5}, PlotRange -> {-1.2, 1.2}]解のさざ波のような揺れは，数値的な不安定から来るものである．安定性を改善するために，偏微分方程式系に人工的な拡散を加えることができる．
偏微分方程式に従属変数が1つだけ含まれる場合にのみ，ペクレ数を使って必要な人工的な拡散を推定することができる．ここでは2つの未知の変数を持つ方程系であるので，手作業で拡散係数 
を選ばなければならない．
ϵ = 3;
Subscript[c, art] = ϵ * h / 2
は，人工的な拡散の影響を抑制する，手動で選ばれたパラメータである．偏微分方程式は人工的な拡散を加えた際に変わるので，ユーザは，忠実性を減らすことなく解を安定化する ϵ の値を1から10までの間で注意深く選ぶ必要がある．
pfun2 = NDSolveValue[{SystemOp == {Subscript[c, art] * D[p[t, x], {x, 2}], Subscript[c, art] * D[v[t, x], {x, 2}]}, Γ, ics}, p, {t, 0, 0.002}, {x}∈mesh]
Plot[pfun2[0.0006, x], {x, 0, 0.5}, PlotRange -> {-1.2, 1.2}]人工的な拡散によって安定性が改善されたことに注目する．次に，数値結果を解析結果と比べる．
pfunAna[t_, x_] := If[x ≤ c * t, Re[Exp[I * (-π / 2)] * Exp[-(I ω x / c)] * Exp[I ω t]], 0]tt = 0.002;
Plot[{pfun[tt, x] - pfunAna[tt, x], pfun2[tt, x] - pfunAna[tt, x]}, {x, 0, 1}, PlotRange -> {-0.015, 0.015}, AxesLabel -> {"x", "error"}, PlotLegends -> {"w/o artificial diffusion", "with artificial diffusion"}]偏微分方程式が変更されるので，人工的な拡散を加えることによって解の勾配の変化は滑らかになる．このため，波の前部でより大きなエラーが起るのである．
Plot[{pfun[tt, x], pfun2[tt, x], pfunAna[tt, x]}, {x, 0.67, 0.72}, PlotRange -> {-0.1, 0.1}, AxesLabel -> {"x", "p"}, PlotLegends -> {"w/o artificial diffusion", "with artificial diffusion", "Analytical"}]偏微分方程式の係数形式
この例では，有限要素法における偏微分方程式の係数形式の重要性が提示されている．単一の従属変数の偏微分方程式で係数形式は以下の方程式(1)で与えられる．
係数についての詳細は，InitializePDECoefficientsを参照されたい．
偏微分方程式は
で定義される．ここで 
は，求める解の従属変数である．係数 
，
，
，
はスカラー，
，
，
はベクトル，
は 
×
行列である．
この例の重要な点は，この方程式（およびその結合バージョン）のみが有限要素法で解ける方程式であるということである．この点を強調するために，より簡単な方程式(2)を考慮する．
方程式とその係数を設定する．偏微分方程式の空間定義域として
を使い，
で停止するように時間積分を設定する．拡散係数 
を1に設定し，これが内部で単位行列に変換される．減衰係数 
として，
の値には1，それ以外では2の値を持つ区分的連続関数を使う．
d[x_] = If[x < 1, 1, 2];
eq1 = d[x] D[u[t, x], t] - D[u[t, x], x, x] == 0;sol1 = NDSolveValue[{eq1, u[0, x] == Exp[x]}, u, {x, 0, 2}, {t, 0, 1}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}]
Plot[sol1[1, x], {x, 0, 2}]
eq2 = D[u[t, x], t] - 1 / d[x]D[u[t, x], x, x] == 0;sol2 = NDSolveValue[{eq2, u[0, x] == Exp[x]}, u, {x, 0, 2}, {t, 0, 1}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}]
Plot[sol2[1, x], {x, 0, 2}]Plot[sol1[1, x] - sol2[1, x], {x, 0, 2}]明らかに解1と解2は異なる．何が起こっているのかを調べるために，NDSolveのデータ構造から実際に解析された偏微分方程式の係数を抽出する小さな補助関数を書く．有限要素法の内部構造についての詳細は，有限要素プログラミングのチュートリアルを参照されたい．
getStateObject[eq_] := NDSolve`ProcessEquations[{eq, u[0, x] == Exp[x]}, u, {t, 0, 1}, {x, 0, 2}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}][[1]]getEquations[eq_] := Module[{temp},
temp = getStateObject[eq];
temp = temp["FiniteElementData"];
temp = temp["PDECoefficientData"];
(# -> temp[#])& /@ {"DampingCoefficients", "DiffusionCoefficients", "ConvectionCoefficients"}]getEquations[eq1]方程式(7)の係数は期待される通りのものである．係数を抽出する代りに，微分方程式がどのように解析されるかを調べる．
GetInactivePDE[getStateObject[eq1]]getEquations[eq2]//SimplifyGetInactivePDE[getStateObject[eq2]]方程式(9)の係数は予想と全く異なるものである．係数
が拡散係数 
に入っている．なぜだろうか．これを理解する鍵は，有限要素法がこのようなタイプの方程式しか解けないということにある．
項
の前には係数がないことに注意されたい．
のタイプの方程式を使うためには，
を 
に再設定し，
を調整して
によって引き起こされる導関数を取り除く必要がある．以下はその例である．
c = h[x];
β = -Div[{{h[x]}}, {x}];
Div[{{c}}.Grad[u[x], {x}], {x}] + β.Grad[u[x], {x}]Div[{{1 / d[x]}}.Grad[u[x], {x}], {x}] - Div[{{1 / d[x]}}, {x}]//Simplifyこれは事実上明示的に 
を指定することと同じである．数学的には，不連続性のため，係数 
はディラクデルタ関数のようなものであり，時間積分ではどのみち削除されてしまう．
eq3 = D[u[t, x], t] + Inactive[Div][{{-1 / d[x]}}.Inactive[Grad][u[t, x], {x}], {x}] == 0;getEquations[eq3]不連続の拡散係数を持つ，時間依存の偏微分方程式については，係数が拡散項から削除されるような形で方程式を並べる方法が最も望ましい．
形式的な偏微分方程式
形式的な偏微分方程式とは，Inactiveコンポーネントを持つ方程式である．Inactiveの偏微分方程式がタイプセットされる場合，Inactiveの偏微分方程式の成分は灰色で表示される．このセクションでは，いつ偏微分方程式が非アクティブでなければならないかについて説明する．
Inactiveの偏微分方程式を使う主な理由は，偏微分方程式係数時期尚早に評価されることを防ぐことにある．例えば，勾配の評価された発散はNeumannValueが何を意味するかに影響を与える．
Div[{{x}}.Grad[u[x], {x}], {x}] == 0Inactive[Div][{{x}}.Inactive[Grad][u[x], {x}], {x}] == 0非アクティブなバージョンでは，方程式に関連して，可能なNeumannValueが何を意味するかについて厳密に制御することができる．アクティブな形式ではこれはできない．NDSolveは引数を未評価のままにすることができないので，NDSolveおよび関連の関数は，自動的にこの操作を行うことはできない．つまり，NDSolveはAttributes，HoldFirst，HoldAllのいずれも持たない．
Inactive形式の方程式が必要となるもう一つの例は，偏微分方程式係数が非対称の場合である．例は，平面応力演算子である．ヤング率とポアソン比については，鉄鋼の値が使われる．
YoungsModulus = 200 10 ^ 9;
pst = (| |
| -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| Inactive[Div][{{-(Y/(1 - ν^2)), 0}, {0, -((Y*(1 - ν))/(2*(1 - ν^2)))}} . Inactive[Grad][u[t, x, ... /(1 - ν^2))}, {-((Y*(1 - ν))/(2*(1 - ν^2))), 0}} . Inactive[Grad][v[t, x, y], {x, y}], {x, y}] |
| Inactive[Div][{{0, -((Y*(1 - ν))/(2*(1 - ν^2)))}, {-((Y*ν)/(1 - ν^2)), 0}} . Inactive[Grad][u[t ... 1 - ν))/(2*(1 - ν^2))), 0}, {0, -(Y/(1 - ν^2))}} . Inactive[Grad][v[t, x, y], {x, y}], {x, y}] |) /. {Y -> YoungsModulus, ν -> 26 / 100};長さ10メートル，高さ0.7メートルで左側で固定された梁の最初の固有値と固有関数を計算する．
ρ = 7800;length = 10;height = 0.7;
{vals, funs} = NDEigensystem[{Thread[ρ * D[{u[t, x, y], v[t, x, y]}, {t, 1}] == pst], DirichletCondition[{u[t, x, y] == 0., v[t, x, y] == 0.}, x == 0]}, {u, v}, t, {x, y}∈Rectangle[{0, 0}, {length, height}], 1];nfreq = Sqrt[vals[[1]]] / (2 π)比較のために，オイラー・ベルヌーイの梁理論を利用した最初の共振周波数を解析的に計算する．
intertia = (width * (height ^ 3)) / 12.;
freq = 0.56 / (length ^ 2) * Sqrt[(intertia * YoungsModulus) / (height * width * ρ)]結果における差分を比べる．差分の約半分は，異なる理論的なアプローチによるものである．
freq - nfreq次に，Activeバージョンの偏微分方程式を調べ，結果の共振周波数を計算する．
ρ = 7800;length = 10;height = 0.7;
{valsWrong, funsWrong} = NDEigensystem[{Thread[ρ * D[{u[t, x, y], v[t, x, y]}, {t, 1}] == Activate[pst]], DirichletCondition[{u[t, x, y] == 0., v[t, x, y] == 0.}, x == 0]}, {u, v}, t, {x, y}∈Rectangle[{0, 0}, {length, height}], 1];
nfreq = Sqrt[valsWrong[[1]]] / (2 π)上の結果は間違っている．何が起こったのだろうか．簡単に言えば，Inactiveではないバージョンは非対称な係数には使えないのである．4つの係数行列のもとの非アクティブな偏微分方程式を見ると，対角ではない2つは対称ではない．
pst言い換えれば，NDEigensystem（およびその他の NDSolve等の関連関数）が意図された形式で偏微分方程式を構文解析する前に，評価されるのである．
pde = Inactive[Div][{{0, 1}, {2, 0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}]Activate[pde]Div[{{0, 1}, {2, 0}}.Grad[u[x, y], {x, y}], {x, y}]評価された偏微分方程式から非対称な係数を回復する方法はない．このため，Inactiveが導入されたのである．非対称な係数に対しては，必ずInactiveバージョンの偏微分方程式を使わなければならない．
構文解析された係数を調べるために，NDEigensystemの時間依存方程式の入力にはNDSolve`ProcessEquationsが使える．NDEigensystemへの入力が演算子形式である場合には，NDSolve`ProcessEquationsの時間依存方程式に変換しなければ，入力に使えない．任意の初期条件
と任意の時間積分区間
を指定しなければならないことに注意する必要がある．
{ndstate} = NDSolve`ProcessEquations[{Thread[ρ * D[{u[t, x, y], v[t, x, y]}, {t, 1}] == pst], DirichletCondition[{u[t, x, y] == 0., v[t, x, y] == 0.}, x == 0], u[0, x, y] == 0, v[0, x, y] == 0}, {u, v}, {t, 0, 1}, {x, y}∈Rectangle[{0, 0}, {length, height}], Method -> {"PDEDiscretization" -> {"MethodOfLines", {"SpatialDiscretization" -> "FiniteElement"}}}];
ndstate["FiniteElementData"]["PDECoefficientData"]["DiffusionCoefficients"]非アクティブな形式の方程式が必要である2つ目の例は，以下の偏微分方程式のように，拡散係数内に従属変数の微分を持つ非線形の偏微分方程式である．
op = Inactive[Div][{{-Derivative[1][u][x]}}. Inactive[Grad][u[x], {x}], {x}];NDSolveValue[{op == 1, DirichletCondition[u[x] == 1., x == 0], DirichletCondition[u[x] == 0., x == 1]}, u, {x}∈Line[{{0}, {1}}], InitialSeeding -> {u[x] == 1 - x}]ここで方程式が非アクティブな形式ではない場合には，NDSolveはこの偏微分方程式を解くことができない．
NDSolveValue[{Activate[op] == 1, DirichletCondition[u[x] == 1., x == 0], DirichletCondition[u[x] == 0., x == 1]}, u, {x}∈Line[{{0}, {1}}], InitialSeeding -> {u[x] == 1 - x}];
現在NDSolveは，係数のいずれかに一次よりも高次の微分を持つ従属変数の定常偏微分方程式を解くことはできない．これは，以下のように説明される．アクティブな演算子は係数形式
に合致する必要がある．まず，拡散係数を 
として設定する．
Activate[op]
Div[{{-2 Derivative[1][u][x]}}.Grad[u[x], {x}], {x}]残りの項
は，係数形式で 
も設定することによって相殺される．
Div[{{-2 Derivative[1][u][x]}}.Grad[u[x], {x}], {x}] + D[2 Derivative[1][u][x], {x}]Derivative[1][u][x] == Activate[op]しかし，
を評価すると，
係数に一次より高次の導関数微分が返される．これは，現在NDSolveの有限要素法では解くことができない．
D[2 Derivative[1][u][x], {x}]有限要素法におけるInactiveの目的は，方程式(11)で与えられる偏微分方程式の係数形式を形式的な偏微分方程式として設定することである．方程式(12)の数学演算子のすべてが非アクティブであってもよい．言い換えれば，この機能はDivとGradを使って表現できるので，Inactiveの偏微分方程式にはInactive[Div]とInactive[Grad]の成分を含むことができるが，Inactive[D]等の非アクティブなDの成分あるいはInactive[Derivative]等の非アクティブなDerivativeの成分は含むことはできない．
ときには，偏微分方程式がNDSolveによってどのように解析されるかを調べると便利なことがある．これは，GetInactivePDEを使って行える．
{state} = NDSolve`ProcessEquations[{D[u[t, x], t] - Sin[x] D[u[t, x], x, x] == Sin[x], u[0, x] == Exp[x]}, u, {t, 0, 1}, {x, 0, 2}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}];
GetInactivePDE[state]より詳しい情報は，偏微分方程式の係数形式のセクションを参照されたい．
ノイマン値と形式的な偏微分方程式
このセクションでは，形式的な偏微分方程式とNeumannValueとの相互関係について説明する．定常的な場合に偏微分方程式は以下のように表わすことができる．
係数についての詳細は，InitializePDECoefficientsを参照されたい．
それに対して，一般化されたノイマン境界値は，値 
を指定する．値 
は，境界の部分上の外向きの法線での流束を指示する．
ノイマン境界値の微分については，「偏微分方程式と境界条件」で説明されている．方程式(13)における係数 
，
，
は，NeumannValueの使用によって直接指定されているわけではないことに注意することが重要である．係数 
，
，
は，方程式(14)で指定されている． NeumannValueを使用すると，
が境界上の 
の値に置き換えられる．
DivとGradの使用のみでは，方程式(15)を完全に指定するには不十分である．以下の指定を考える．
が
として与えられ，
が単位行列，
がゼロベクトルである場合を考える．発散の評価の結果，以下が返されることに注意する．
ここで 
は 
，
は 
の発散（この例では
）である． 方程式(16)と(17)は等しい．
c = IdentityMatrix[2];
α = {x, y};
Div[-c.Grad[u[x, y], {x, y}] - α u[x, y], {x, y}] == Div[-c.Grad[u[x, y], {x, y}], {x, y}] - α.Grad[u[x, y], {x, y}] - Div[α, {x, y}] u[x, y]しかし，評価で問題が起こる．方程式は等しいが，ノイマン指定は等しくない．方程式(18)のノイマン値は，陰的に以下に変更されている．
このことについてさらに説明するために，「有限要素プログラミング」のチュートリアルで説明されるように，偏微分方程式を解析し，有限要素データから偏微分方程式データを抽出するヘルパー関数を作成する．
getPDECoefficients[eqn_] := First[NDSolve`ProcessEquations[
{eqn, DirichletCondition[u[x, y] == 0, True]}, u, {x, y}∈Disk[]]]["FiniteElementData"]["PDECoefficientData"]pdec = getPDECoefficients[Div[-c.Grad[u[x, y], {x, y}] - α * u[x, y], {x, y}] == 0]PDECoefficientDataは，解析された偏微分方程式係数を含むデータ構造である．
pdec["DiffusionCoefficients"]
pdec["ConvectionCoefficients"]
pdec["ReactionCoefficients"]pdec["ConservativeConvectionCoefficients"]評価中に，係数 
は，
と 
の項に分割される．微分方程式を入力する正しい構造は，Inactivateによって，あるいは直接Inactiveを使用することによって評価される係数である．
ipde = Inactive[Div][-c.Inactive[Grad][u[x, y], {x, y}] + Inactive[Times][-α, u[x, y]], {x, y}] == 0ipdec = getPDECoefficients[ipde]
ipdec["DiffusionCoefficients"]
ipdec["ConservativeConvectionCoefficients"]ipdec["ConvectionCoefficients"]ipdec["ReactionCoefficients"]γ の項がある場合には，以下の非アクティブ偏微分方程式を使うこともある．
γ = {y, x ^ 2};
ipde = Inactive[Div][Inactive[Plus][
-c.Inactive[Grad][u[x, y], {x, y}],
Inactive[Times][-α, u[x, y]],
γ
], {x, y}] == 0ipdec = getPDECoefficients[ipde]
ipdec["LoadDerivativeCoefficients"]以下の例では，実用的な面から見た，ノイマン値についてアクティブと非アクティブの偏微分方程式の意味の違いを説明する．この例では，バイアスポテンシャル 
について，領域 
に限定される拡散を示す．
Ω = ImplicitRegion[-5 ≤ x ≤ 5 && -5 ≤ y ≤ 5 && x ^ 2 + y ^ 2 ≥ 1 ^ 2, {x, y}];
bp = Function[{x, y}, -E ^ (-x ^ 2 - y ^ 2)];
α = Grad[bp[x, y], {x, y}];ディリクレ境界条件は，領域の上下(
)をそれぞれポテンシャル1と0で設定する．
bc = {DirichletCondition[u[x, y] == 1, y == -5], DirichletCondition[u[x, y] == 0, y == 5]};残りの2つの辺(
)には，自然な境界条件（ノイマンゼロ）が適用される．つまり，この2つの辺では，流束はゼロである．
op = Inactive[Div][Inactive[Plus][-Inactive[Grad][u[x, y], {x, y}], Inactive[Times][-α, u[x, y]]], {x, y}]ufunI = NDSolve[{op == 0, bc}, u, {x, y}∈Ω];ipdec = getPDECoefficients[op == 0];
{ipdec["DiffusionCoefficients"], ipdec["ConservativeConvectionCoefficients"]}
Plot[Evaluate[{
(Derivative[1, 0][u][5, y] + u[5, y] Derivative[1, 0][bp][5, y]), (Derivative[1, 0][u][-5, y] + u[-5, y] Derivative[1, 0][bp][-5, y])
} /. ufunI], {y, -5, 5}]
におけるゼロ流束の境界条件（2つの辺上でゼロ流束）で，もとの項がない場合には，
におけるボックスへの流束は，
おけるボックスからの流束と同じはずである．
Plot[Evaluate[Subtract[
(Derivative[0, 1][u][x, 5] + u[x, 5] Derivative[0, 1][bp][x, 5]), (Derivative[0, 1][u][x, -5] + u[x, -5] Derivative[0, 1][bp][x, -5])
] /. ufunI], {x, -5, 5}]ボックスへの流束は，ボックスからの流束と等しい．比較のために，アクティブ偏微分方程式を解く．
ufunA = NDSolve[{Activate[op] == 0, bc}, u, {x, y}∈Ω];apdec = getPDECoefficients[Activate[op] == 0];
{apdec["DiffusionCoefficients"], apdec["ConservativeConvectionCoefficients"]}
{apdec["ConvectionCoefficients"], apdec["ReactionCoefficients"]}
Plot[Evaluate[Function[y,
(Derivative[0, 1][u][x, y]) /. ufunA] /@ {5, -5}], {x, -5, 5}, PlotRange -> {-0.2, 0}]ときには，偏微分方程式がNDSolveによってどのように解析されるかを調べると便利なことがある．これは，GetInactivePDEを使って行える．
{state} = NDSolve`ProcessEquations[{D[u[t, x], t] - Sin[x] D[u[t, x], x, x] == Sin[x], u[0, x] == Exp[x]}, u, {t, 0, 1}, {x, 0, 2}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}];
GetInactivePDE[state]より詳しい情報は，偏微分方程式の係数形式のセクションを参照されたい．
偏微分方程式係数の効率的な評価
偏微分方程式の係数は定数一定でなくてもよい．この例では，拡散係数 
はシミュレーションの領域における位置によって変わる．
指定された係数を評価する順序を理解することが重要である．順序がシミュレーションのパフォーマンスに影響を及ぼすことがあるからである．まず，内部境界のあるメッシュを生成する．
mesh = ToElementMesh[Annulus[], "RegionHoles" -> None, MaxCellMeasure -> 0.001];mesh["Wireframe"["MeshElement" -> "BoundaryElements"]]メッシュを設定したら，非特定の拡散係数 
を持つ偏微分方程式を設定する．
pde = Inactive[Div][c*Inactive[Grad][u[x, y], {x, y}], {x, y}] == 1;bc = DirichletCondition[u[x, y] == 1, x ^ 2 + y ^ 2 == 1];拡散係数 
は，内部領域では拡散係数の値が 
となり，その領域の外側では拡散係数が 
となるように設定される．
radius = x ^ 2 + y ^ 2 ≤ 1 / 2 ^ 2;
cInner = 1;
cOuter = 10;
c = If[radius, cInner, cOuter]time1 = AbsoluteTiming[NDSolveValue[{pde, bc}, u, {x, y}∈mesh];][[1]]シミュレーションは予期されるように行われた．拡散係数 
をより詳しく見てみると，If文の分岐が評価されなかったことが分かる．
cIf文がその分岐引数を評価しないのは，If文の属性が評価を妨げる設定になっているからである．
Attributes[If]拡散係数 
におけるIf文の分岐引数が評価されないということは，呼出しのたびに変数が検索されなければならないことを意味する．その結果，If文をコンパイルすることができず，標準評価器が使われる．この場合，この評価器を使った方が時間がかかる．
係数関数のすべての引数が評価されるように関数を作成するには，Withを利用することができる．
c = With[{cInner = cInner, cOuter = cOuter}, If[radius, cInner, cOuter]]これでIf文のすべての引数が評価されていることが見て取れる．
pde = Inactive[Div][c*Inactive[Grad][u[x, y], {x, y}], {x, y}] == 1;
time2 = AbsoluteTiming[NDSolveValue[{pde, bc}, u, {x, y}∈mesh];][[1]]time1 / time2ここで使ったような小さい例ではこれはあまり大きな問題ではないように見えるが，大きなシミュレーションにおいては，完全に評価された係数は考慮に入れるべきことである．
If文に内容を入れるためにWithを使うよりも簡単に使えるもう一つの方法に，Piecewiseを使う方法がある． Piecewiseはそれに関連付けられた属性を持たないので，その引数を評価する．
c = Piecewise[{{cInner, radius}}, cOuter]Piecewise文ではすべての引数が評価されることに注意する．
pde = Inactive[Div][c*Inactive[Grad][u[x, y], {x, y}], {x, y}] == 1;
time3 = AbsoluteTiming[NDSolveValue[{pde, bc}, u, {x, y}∈mesh];][[1]]time1 / time3スピード面では，Piecewise構造はWith/If構造と同じくらい効率的である．
HeatTransferPDEComponentのような偏微分方程式成分を使って偏微分方程式が設定されている場合には，効率的な係数の生成が自動的に行われることに注意する．
pars = <||>;
pars["MassDensity"] = Piecewise[{{Subscript[ρ, tungsten], y <= 0}}, Subscript[ρ, steel]];
pars["SpecificHeatCapacity"] = Piecewise[{{Subscript[Cp, tungsten], y <= 0}}, Subscript[Cp, steel]];
pars["ThermalConductivity"] = Piecewise[{{Subscript[k, tungsten] * IdentityMatrix[2], y <= 0}}, Subscript[k, steel] * IdentityMatrix[2]];pars[Subscript[ρ, tungsten]] = 19000;
pars[Subscript[Cp, tungsten]] = 134;
pars[Subscript[k, tungsten]] = 163;
pars[Subscript[ρ, steel]] = 7500;
pars[Subscript[Cp, steel]] = 470;
pars[Subscript[k, steel]] = 44;HeatTransferPDEComponent[{T[x, y], {x, y}}, pars]記号のパラメータ値が実際の数値に置き換えられたことに注意する．
係数がSetDelayedの後ろに隠れている場合にもほとんど同じ問題が起る．
regionMember = RegionMember[Line[{{-1}, {1}}]];
ifCoefficient[x_] := If[regionMember[{x}], 1, 2];
pwCoefficient[x_] := Piecewise[{{1, regionMember[{x}]}}, 2];? ifCoefficient? pwCoefficientifCoefficient2[x_] := Evaluate[If[regionMember[{x}], 1, 2]];
pwCoefficient2[x_] := Evaluate[Piecewise[{{1, regionMember[{x}]}}, 2]];? ifCoefficient2? pwCoefficient22つ目の場合には，偏微分方程式モデル内の係数として使われるときに，検定は記号的に評価され，より効率的である．
偏微分方程式の係数内の補間関数の効率的な評価
偏微分方程式の係数も補間関数になり得る．偏微分方程式が解かれるメッシュと同じものを補間関数が利用する場合には，補間関数を非常に効率的に評価することが可能である．
mesh = ToElementMesh[Rectangle[], MaxCellMeasure -> 0.0005]data = Function[{x, y}, {x, y, Sqrt[x ^ 2 + y ^ 2]}]@@@mesh["Coordinates"];
fun1 = Interpolation[data, "InterpolationOrder" -> 1]RepeatedTiming[NDSolveValue[{-Laplacian[u[x, y], {x, y}] + fun1[x, y] * u[x, y] == 1, DirichletCondition[u[x, y] == 0, x == 0]}, u, Element[{x, y}, mesh]];]fun1["ElementMesh"] === meshfun2 = ElementMeshInterpolation[mesh, Sqrt[Total[mesh["Coordinates"] ^ 2, {2}]]]RepeatedTiming[NDSolveValue[{-Laplacian[u[x, y], {x, y}] + fun2[x, y] * u[x, y] == 1, DirichletCondition[u[x, y] == 0, x == 0]}, u, Element[{x, y}, mesh]];]fun2["ElementMesh"] === mesh補間関数のメッシュが，偏微分方程式に使われたメッシュと同じである場合には，計算に必要な時間は短くなる．
最も効率的な方法は，可能であれば，式を直接与えることである．
RepeatedTiming[sol = NDSolveValue[{-Laplacian[u[x, y], {x, y}] + Sqrt[x ^ 2 + y ^ 2] * u[x, y] == 1, DirichletCondition[u[x, y] == 0, x == 0]}, u, Element[{x, y}, mesh]];]NeumannValueと境界導関数の関係
以下の例では，NeumannValueと境界導関数の関係を説明し，パラメータが変更しやすいように設定する．まず，領域
を指定する．
pos1 = 1;pos2 = 2;
Ω = ImplicitRegion[True, {{x, pos1, pos2}}]
偏微分方程式は，従属変数 
を1つ持つ．これは一次元系であるので，拡散係数 
は1×1行列であり，これは数として指定できる．
c = -2;
op = D[c u[x], {x, 2}]
r = 1;
Subscript[Γ, D] = DirichletCondition[u[x] == r, x == pos1];g = -2;q = -1;
Subscript[Γ, N] = NeumannValue[g - q u[x], x == pos2];nufun = NDSolveValue[{op == Subscript[Γ, N], Subscript[Γ, D]}, u, {x}∈Ω]aufun = DSolveValue[{op == 0, u[pos1] == r, -c u'[pos2] == g - q u[pos2]}, u, x]方程式のDSolveValueの形式化において，一般化されたノイマン境界方程式がどのように明示的に与えられなければならないか，そしてそれらの係数がどのように偏微分方程式の係数に一致しているかに注意しなければならない．
解析結果を数値結果と比べるために，この2つの差分のプロットを作成する．
Plot[{aufun[x] - nufun[x]}, {x, pos1, pos2}]従属変数名の順序
偏微分方程式の結合系を分析する場合に，NDSolveで実装された有限要素法が予期しない結果を返すということがある．このセクションでは，どのような場合にこれが起こるのかを説明し，この問題を避ける方法を提示する． 最後に，さらに詳しく知りたいという読者のために，より詳細な説明を提供する．
Ω = Rectangle[{0, 0}, {1, 1}];eqn = {
Laplacian[u[x, y], {x, y}] - D[w[x, y], x] == 0,
Laplacian[v[x, y], {x, y}] - D[w[x, y], y] == 0,
D[u[x, y], x] + D[v[x, y], y] == 0};bcs = {u[x, 0] == 0, u[x, 1] == 0, v[x, 0] == 0, v[x, 1] == 0, w[1, y] == 0, w[0, y] == 1};pde = {eqn, bcs};この例における従属変数は，u，v，w である．系を解き，従属変数 w を可視化する．
{xVel, yVel, pressure} = NDSolveValue[pde, {u, v, w}, {x, y}∈Ω];
DensityPlot[pressure[x, y], {x, y}∈Ω]次の例では，従属変数の名前を w から p に変更し，方程式系をもう一度解く．
{xVel, yVel, pressure} = NDSolveValue[pde /. w -> p, {u, v, p}, {x, y}∈Ω];
DensityPlot[pressure[x, y], {x, y}∈Ω]その結果は予期しないものとなった．何が起こったかを詳しく説明する前に，予期される動作が必ず起こるようにする2つの方法について説明する．
1つ目は，NDSolveの"DependentVariable"メソッドオプションを通して従属変数を指定する方法である．
{xVel, yVel, pressure} = NDSolveValue[pde /. w -> p, {u, v, p}, {x, y}∈Ω, DependentVariables -> {u, v, p}];
DensityPlot[pressure[x, y], {x, y}∈Ω]"DependentVariable"メソッドオプションを指定することは，従属変数が指定の次数で離散化されるように強制する．このことについては，下で詳しく説明する．
この例の場合には，NDSolveは，偏微分方程式の離散化の際に従属変数が使うべき補間次数を指定するための有限要素法のオプションを与える．補間次数の指定は，陰的に方程式と従属変数を同じ次数に強制する．
{xVel, yVel, pressure} = NDSolveValue[pde /. w -> p, {u, v, p}, {x, y}∈Ω, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];
DensityPlot[pressure[x, y], {x, y}∈Ω]通常，それぞれの従属変数に同じ補間次数を指定しても問題ない．補間次数の低減も問題の根幹にあるものではない．この場合には，方程式は流体方程式（ストークスの流れ）であり，安定化は，圧力変数のInterpolation次数が速度 u および v の補間次数よりも1小さい場合には強制することができる．
根底にある問題が何であるかを知るために，両方の場合について構築された離散系の行列を表示する．
NDSolveValue[pde /. w -> p, {u, v, p}, {x, y}∈Ω, Method -> {"FiniteElement", "PDESolveOptions" -> {"LinearSolver" -> {(
Print[MatrixPlot[#]];
LinearSolve[##])&}}}];
NDSolveValue[pde /. w -> p, {u, v, p}, {x, y}∈Ω,
DependentVariables -> {u, v, p}, Method -> {"FiniteElement", "PDESolveOptions" -> {"LinearSolver" -> {(
Print[MatrixPlot[#]];
LinearSolve[##])&}}}];
この動作をさらに調べるために，NDSolveの状態オブジェクトを作成する．
{state1} = NDSolve`ProcessEquations[pde /. w -> p, {u, v, p}, {x, y}∈Ω, Method -> {"FiniteElement"}];FiniteElementDataのFEMMethodDataの従属変数の順序は{p,u,v}であることに注意する．
state1["FiniteElementData"]["FEMMethodData"]["VariableData"]PDECoefficientDataの拡散と対流の係数を調べると，結合偏微分方程式の係数の構造が分かる．
state1["FiniteElementData"]["PDECoefficientData"]["DiffusionCoefficients"]//MatrixFormstate1["FiniteElementData"]["PDECoefficientData"]["ConvectionCoefficients"]//MatrixForm有限要素の補間次数法のオプションが設定されると，係数構造は，従属変数の正規順序に調整される．
{state2} = NDSolve`ProcessEquations[pde /. w -> p, {u, v, p}, {x, y}∈Ω, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];state2["FiniteElementData"]["PDECoefficientData"]["DiffusionCoefficients"]//MatrixFormstate2["FiniteElementData"]["PDECoefficientData"]["ConvectionCoefficients"]//MatrixForm実際，方程式の入力の順序を変更して，従属変数の正規順序にマッチさせると，予期される通りの結果が得られる．
{pressure, xVel, yVel} = NDSolveValue[{eqn[[{3, 1, 2}]], bcs} /. w -> p, {p, u, v}, {x, y}∈Ω];
DensityPlot[pressure[x, y], {x, y}∈Ω]それでは，なぜこの順序変更が自動的に行われないのであろうか．
先天的な問題として，指定された境界条件は，係数の対角成分にしか適用できないということがある．これは，境界条件が入力されたときに，特定の境界条件がどのように偏微分方程式の特定の方程式に関連付けられるかということの情報が失われるためである．
bc1 = Flatten[state1["FiniteElementData"]["BoundaryConditionData"]["Stationary"]];
Transpose[{bc1[[All, 1]][[All, 2]], bc1[[All, -1]]}]ここでは，それぞれのサブリストは2つの部分に分かれている．最初の部分は，離散化の際に境界条件の2つ目の部分が適用される行と列を示す．
もう一つの方法は，ディリクレ境界条件をそれぞれが関係している方程式に関連付ける方法である．
pde2 = Thread[Equal[eqn[[All, 1]], {DirichletCondition[u[x, y] == 0, y == 0] + DirichletCondition[u[x, y] == 0, y == 1],
DirichletCondition[v[x, y] == 0, y == 0] + DirichletCondition[v[x, y] == 0, y == 1],
DirichletCondition[w[x, y] == 1, x == 0] + DirichletCondition[w[x, y] == 0, x == 1]
}]];{xVel, yVel, pressure} = NDSolveValue[pde2 /. w -> p, {u, v, p}, {x, y}∈Ω];
DensityPlot[pressure[x, y], {x, y}∈Ω]{xVel, yVel, pressure} = NDSolveValue[pde /. w -> p, {u, v, p}, {x, y}∈Ω, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];
DensityPlot[pressure[x, y], {x, y}∈Ω]解の検証
NDSolveや関連関数から得た微分方程式の解を検証する方法として，手元の微分方程式に解を再挿入してみるということがある．
bc = DirichletCondition[u[x] == 0, True];op = -Laplacian[u[x], {x}] - 1;
sol = NDSolveValue[{op == 0, bc}, u, {x, 0, π}]verify = op /. u -> solPlot[verify, {x, 0, π}]上の例を見ると，微分方程式に解を再挿入する方法が，微分方程式の解を検証するのに好ましいと思われるかも知れないが，必ずしもそうとは言えない．
op = -Laplacian[u[x], {x}] + u[x] - 1;
sol = NDSolveValue[{op == 0, bc}, u, {x, 0, π}];
Plot[Evaluate[op /. u -> sol], {x, 0, π}]この結果は，解が前の例ほど正確ではないことを示唆する．実際にそうだろうか．
前の例では，解が微分方程式に再挿入された場合に，InterpolatingFunctionのより高次な微分が計算されなければならないという陰的な仮定がされた．しかし，微分回復の過程に誤差が伴わないとは限らず，
の尺度までうまく計算できない．二次よりも高次のメッシュを使うと状況が改善され，一次のメッシュを使うとより悪化する．試してみよう．
Method -> {"FiniteElement", "MeshOptions" -> {"MeshOrder" -> 1}}理想的には，微分方程式の解は，解の検証中にさらなる数値誤差をもたらさずにチェックすることが望ましい．そのためには，作られた解のメソッドが便利である．微分方程式の解析解が存在していたなら，計算された改はこの解析解に対して確かめることができる．作られた解のアイディアは，任意解を選ぶためのものであり，選ばれた解が実際の解であるように，もとの微分方程式を修正するためのものである．
arbitrarySol = u -> Sing = op /. arbitrarySol選ばれた解析解をもとの微分方程式に挿入すると，新しい右辺の項 
ができる．もとの方程式がこの新しい右辺 
に等しくなるように設定されると，選択された解であるSinは，新しい微分方程式の解に構築される．
sol = NDSolveValue[{op == g, bc}, u, {x, 0, π}];Plot[Evaluate[sol[x] - u[x] /. arbitrarySol], {x, 0, π}]この選ばれた解析解の場合の誤差は，解の形式NDSolveValueをもとの方程式に再挿入した場合のものよりもずっと小さくなる．この方法の利点は，解の検証中に新たな数値誤差がもたらされないということである．