NDEigensystem::fembdcc NDEigenvalues::fembdcc NDSolve::fembdcc NDSolveValue::fembdcc ParametricNDSolve::fembdcc ParametricNDSolveValue::fembdcc InitializePDECoefficients::femcnmd
例題
例 (2)
k = Piecewise[{{{{10, 0}, {0, 10}}, x <= 1 / 2},
{{{50 * T[x, y] + 1, 0}, {0, 50}}, True}}];以下の偏微分方程式は,与えられた拡散係数に対してメッセージを発する:
NDSolveValue[{Inactive[Div][-k.Inactive[Grad][T[x, y], {x, y}], {x, y}] == 1, DirichletCondition[T[x, y] == 0, x == 0], DirichletCondition[T[x, y] == 1, x == 1]}, T, {x, y}∈Rectangle[]]
これは,自動的な線形化の際に,係数の導関数が形成できないことから起る.解決策として,係数をPiecewiseを返す行列としてではなく,Piecewiseのスカラー係数を含む行列として書くことができる:
k = {{Piecewise[{{10, x <= 1 / 2},
{50 * T[x, y] + 1, True}}], 0}, {0, Piecewise[{{10, x <= 1 / 2},
{50, True}}]}};;NDSolveValue[{Inactive[Div][-k.Inactive[Grad][T[x, y], {x, y}], {x, y}] == 1, DirichletCondition[T[x, y] == 0, x == 0], DirichletCondition[T[x, y] == 1, x == 1]}, T, {x, y}∈Rectangle[]]偏微分方程式係数として使われる,以下のコンパイル関数を考える:
cf = Compile[{rho, z}, rho + z];NDSolveValue[{1 == cf[rho, z] D[u[rho, z], z, z], DirichletCondition[u[rho, z] == 0, True]}, u, {rho, 0, 1}, {z, 0, 1}]最初の警告メッセージは,コンパイル関数が機械サイズの実数以外のもので評価されたことを示す.実際,このメッセージは,記号の引数でコンパイル関数を評価したことから発するものである:
cf[rho, z]このメッセージについては,ここにより詳しく説明されている.コンパイル関数は,コンパイルされたコードの本体を記号式として返し,その後NDSolveが偏微分方程式を解く.NDSolveが返す警告メッセージは,これが対流支配の偏微分方程式であることを示し,そのことについてはここに説明がある.メッセージを発しながらも,NDSolveが偏微分方程式の解を返すということに注意されたい.
記号入力について評価しないコンパイル関数を指定して,このような問題に対処せずに済むようにしたい場合もあるかもしれない:
cf = Compile[{rho, z}, rho + z, RuntimeOptions -> EvaluateSymbolically -> False];
cf[rho, z]cfをNDSolveで偏微分方程式の係数として使うと問題が起る.なぜ問題が起るのか,どのようにしてこれを避けることができるのかを理解するためには,少し脇道に逸れる必要がある.
偏微分方程式の設定は有限要素法をトリガする.何が有限要素法をトリガするのかについては,ここに説明がある.有限要素法は,ここに説明されているように,非常に特別の形式の偏微分方程式のみを扱うことができる:
この偏微分方程式を解くために起ることは,係数
が発散に引き込まれ,
項を調整することで以下のように補われるのである:
記号的に評価されないコンパイル関数は,この偏微分方程式の再構築の過程で問題を起す:
NDSolveValue[{1 == cf[rho, z] D[u[rho, z], z, z], DirichletCondition[u[rho, z] == 0, True]}, u, {rho, 0, 1}, {z, 0, 1}]
D[cf[rho, z], z]このためNDSolveがエラーを発する.これに対処する方法は2つある.1つ目は,可能であれば,記号係数を与える方法である.そうすると,NDSolveは解を求めることができる:
solution1 = NDSolveValue[{1 == (rho + z) D[u[rho, z], z, z], DirichletCondition[u[rho, z] == 0, True]}, u, {rho, 0, 1}, {z, 0, 1}]もう一つは,手動で導関数を構築して,係数を指定する方法である:
dcf = Div[(rho + z) * {{0, 0}, {0, 1}}, {rho, z}]solution2 = NDSolveValue[{1 == Inactive[Div][(cf[rho, z] * {{0, 0}, {0, 1}}).Inactive[Grad][u[rho, z], {rho, z}], {rho, z}] - dcf.Inactive[Grad][u[rho, z], {rho, z}] , DirichletCondition[u[rho, z] == 0, True]}, u, {rho, 0, 1}, {z, 0, 1}]