WOLFRAM

Wolfram言語 ＆ システムドキュメントセンター

偏微分方程式の数値解法

Wolfram言語関数NDSolveには偏微分方程式を解くための広範な機能が備わっている．NDSolve独自の機能として，偏微分方程式と解の領域が記号形式で与えられると，NDSolveは問題の構造に最も適していると思われる数値解法を自動的に選ぶというものがある．通常，アルゴリズムの自動選択はうまくいくが，使用されるメソッドを理解しておくと便利である．そのことが，提供された解をよりよく理解し，よりよい解を得たり効率を向上させたりするのにメソッドオプションを使うのに役立つ．

NDSolve[{eqn1,eqn2,…},u,{x,xmin,xmax},{y,ymin,ymax}]矩形領域[xmin,xmax]×[ymin,ymax]上で関数 u[x,y]に対する{eqn1,eqn2,…}の数値解を求める
NDSolve[{eqn1,eqn2,…},u,{x,y}∈Ω]領域Ω上で関数 u[x,y]に対する{eqn1,eqn2,…}の数値解を求める
NDSolve[{eqn1,eqn2,…},u,{t,tmin,tmax},{x,y}∈Ω]時間 tmin と tmax の間で，領域Ω上で関数 u[t,x,y]に対する時間依存方程式{eqn1,eqn2,…}の数値解を求める

NDSolveで偏微分方程式の数値解を求める

NDSolveは，"TensorProductGrid"メソッドを介して表わされる有限要素法と有限差分法を使って，変微分方程式を離散化して解く．時間依存方程式の場合は，有限要素法あるいは有限差分法の空間離散化とともに，数値解法の線の方法が使われる． 線の方法については「線の方法」で詳述する．有限要素法はより汎用である．これに関してはそのチュートリアルで詳しく説明する．有限要素法は，定常偏微分方程式と時間依存偏微分方程式のどちらでも利用できるが，"TensorProductGrid"法は時間依存偏微分方程式でしか利用できない．

以下の例では，どのように方程式が入力され，どのようにメソッドが自動的に選ばれるかを簡単に説明する．

解を検証するのに便利なパッケージをロードする：
Wolfram Language code: Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
ゼロ境界条件の矩形領域上で，2変数のポワソン(Poisson)方程式を解く：
Wolfram Language code: uval = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == Exp[-(x ^ 2 + 10 y ^ 2)], u[x, -1] == u[x, 1] == u[-1, y] == u[1, y] == 0}, u, {x, -1, 1}, {y, -1, 1}]
Wolfram Language code: ContourPlot[uval[x, y], {x, -1, 1}, {y, -1, 1}]

上の入力は明示的に与えられた導関数の表記を使っている．問題によっては演算子形式を使った方がより表現しやすく正確である場合がある．デフォルトでは，一般的な演算子Grad，Div，Curl，Laplacianは自動的に評価されるが，変微分方程式を表現する場合は通常これらが評価されないようにした方がよい．この場合，Inactiveを使って演算子形式を維持することができる．一般に，Inactive[f]は，f の評価されない非アクティブな形式を表す．Inactive形式はActivateを使ってアクティブにすることができるが，NDSolveのようなコマンドも直接非アクティブな形式で使える．

微分演算子が評価されたInactiveな形式：
Wolfram Language code: Div[a[x, y] Grad[u[x, y], {x, y}], {x, y}]
Wolfram Language code: Inactive[Div][a[x, y] Inactive[Grad][u[x, y], {x, y}], {x, y}]

非アクティブな形式は，演算子形式の明示的な評価で失われた可能性のある意味を保持していることに注意する．例えば， が非連続ならば，導関数の式は不定であるが， が の不連続性に対応する不連続性を持つ場合は演算子式は正確に定義され得る．

演算子を非アクティブにするためには，第2引数のパターンで与えられた，非アクティブにしたい演算子にInactivateを使う．

Inactivateを使う：
Wolfram Language code: Inactivate[Div[a[x, y] Grad[u[x, y], {x, y}], {x, y}], Div | Grad]

方程式を入力する場合に最も便利で最も読みやすい方法は，タイプセット形式を使うことであろう．これは非アクティブな形式の出力をコピーペーストするか，それぞれの演算子のレファレンスページで与えられるエスケープ列を使うかすることで行える．

ラプラシアン(Laplacian)のタイプセット形式：
Wolfram Language code: Inactive[Laplacian][u[x, y], {x, y}]
Wolfram Language code: Activate[%]

上記の例では，明示的方程式としてゼロ境界条件が与えられた．他の方法としてDirichletCondition[eqn,pred]を使うものがある．これは pred がTrueである領域の境界上ならどこでも eqn が満足されなければならないことを指定する．一般化されたノイマン(Neumann)境界条件はNeumannValueを使って指定できる．これは有限要素法に特有のものなので，NeumannValueの説明は有限要素法のチュートリアルで説明されている．

通常，空間変数はある領域に限られており，NDSolveは表記を理解する．

ポワソン方程式のきれいなタイプセット形式を解く：
Wolfram Language code: Ω = Rectangle[{-1, -1}, {1, 1}]; Subscript[Γ, D] = DirichletCondition[u[x, y] == 0, True];
Wolfram Language code: NDSolveValue[{Inactive[Laplacian][u[x, y], {x, y}] == E^-(x^2 + 10 y^2), Subscript[Γ, D]}, u, {x, y}∈Ω]

この解は上記で示したものと同じであり，どちらも有限要素法を使って計算されている．この例では と の両方に境界条件があるので，自動選択機能でこの方法が選ばれた．線の方法では，独立変数のうちの一つに対して初期値問題が必要とされる．

有限要素法は，解を計算するのにメッシュを使う．メッシュは解から抽出して見ることができる．

解からメッシュを抽出し，そのワイヤフレームを表示する：
Wolfram Language code: {mesh} = InterpolatingFunctionCoordinates[uval]
Wolfram Language code: mesh["Wireframe"]

有限要素法は非矩形の図形もうまく扱う．異なる領域上で問題を解くのに必要なのは，その領域を定義することだけである．

単位円板上でポワソン方程式を解く：
Wolfram Language code: Ω = Disk[]; Subscript[Γ, D] = DirichletCondition[u[x, y] == 0, True];
Wolfram Language code: duval = NDSolveValue[{Inactive[Laplacian][u[x, y], {x, y}] == E^-(x^2 + 10 y^2), Subscript[Γ, D]}, u, {x, y}∈Ω]

可視化コマンドの中には領域の表記を受け入れるものもあるので，それを直接使って解を可視化することができる．

解を表示する：
Wolfram Language code: ContourPlot[duval[x, y], {x, y}∈Ω]
解を見付けるために使用したメッシュを表示する：
Wolfram Language code: {mesh} = InterpolatingFunctionCoordinates[duval]; mesh["Wireframe"]

オプションMesh->Allを使って，プロットと一緒にメッシュを見ることができる場合もある．

解の曲面のプロットとメッシュを一緒に表示する：
Wolfram Language code: Plot3D[duval[x, y], {x, y}∈Ω, Mesh -> All]

NDSolveは独立変数を持つ方程式を扱う．その独立変数は，解がその変数内で，有限差分あるいは有限要素空間離散化で数値的線の方法を使って解かれる，指定された初期条件から進化するという点で，時間のようなものである．

例として，板の中心の温度が，一定の熱が外側から加えられ始めてから指定された温度に到達するのにかかる時間を求める問題を考えてみる．これは熱方程式でモデル化される．

矩形板での板加熱問題を解く：
Wolfram Language code: ph = NDSolveValue[{D[u[t, x, y], t] == D[u[t, x, y], x, x] + D[u[t, x, y], y, y], u[0, x, y] == If[-1 < x < 1 && -1 < y < 1, 0, 1], u[t, -1, y] == u[t, 1, y] == u[t, x, -1] == u[t, x, 1] == 1, WhenEvent[u[t, 0, 0] == .5, {T = t, "StopIntegration"}]}, u, {t, 0, 10}, {x, -1, 1}, {y, -1, 1}]; T
最終の時間で解をプロットする：
Wolfram Language code: Plot3D[ph[T, x, y], {x, -1, 1}, {y, -1, 1}]

領域が矩形の場合，NDSolveはデフォルトで有限差分法を使う．Methodオプションを使うと型を指定することや空間離散に対する特定の詳細を指定することも可能である．チュートリアルの「線の方法」および「有限要素法」には，オプションについての詳細とそれを使った多数の例題が記載されている．

有限要素を使って，矩形板の板加熱問題を解く：
Wolfram Language code: ph = NDSolveValue[{D[u[t, x, y], t] == D[u[t, x, y], x, x] + D[u[t, x, y], y, y], u[0, x, y] == If[-1 < x < 1 && -1 < y < 1, 0, 1], u[t, -1, y] == u[t, 1, y] == u[t, x, -1] == u[t, x, 1] == 1, WhenEvent[u[t, 0, 0] == .5, {T = t, "StopIntegration"}]}, u, {t, 0, 10}, {x, -1, 1}, {y, -1, 1}, MaxStepFraction -> 1, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> "FiniteElement"}}];T

空間geometryが不規則ならば，NDSolveは自動的に有限要素を使う．

円板について解く：
Wolfram Language code: Ω = Disk[]; Subscript[Γ, D] = DirichletCondition[u[t, x, y] == 1, True];
Wolfram Language code: cph = NDSolveValue[{Subscript[∂, t]u[t, x, y] == Inactive[Laplacian][u[t, x, y], {x, y}], Subscript[Γ, D], u[0, x, y] == If[x ^ 2 + y ^ 2 < 1, 0, 1], WhenEvent[u[t, 0, 0] == .5, {T = t, "StopIntegration"}]}, u, {t, 0, 1}, {x, y}∈Ω];T
最後に円板についての解をプロットする：
Wolfram Language code: Plot3D[cph[T, x, y], {x, y}∈Ω, Mesh -> All]