有限要素プログラミング
はじめに
NDSolveは，有限要素法を使って偏微分方程式を解くための，高レベルでワンステップのインターフェースを提供する．しかし，解くプロセスのステップを詳しく制御したいという場合もある．NDSolve`FEM`パッケージは，解くプロセスの各ステップのおいて広範な制御を与える，低レベルインターフェースを提供する．
Needs["NDSolve`FEM`"]NDSolve`FEM`パッケージの低レベル関数は，さまざまな用途に使用できる．
- NDSolveを方程式のプリプロセッサとして使う
NDSolve内の有限要素データ
NDSolveフレームワークには，NDSolve`ProcessEquations，NDSolve`Iterate，NDSolve`ProcessSolutionsという低レベル関数と，すべてのタイプの微分方程式に対して解法プロセスへの何らかのアクセスを提供するデータオブジェクトNDSolve`StateDataが含まれる．これらの関数の一般的な用法は，「コンポーネントとデータ構造」に説明されている．
残りのNDSolveと一貫して，有限要素法は，そのデータにNDSolve`StateDataを通してアクセスすることができる．次のセクションでは，NDSolve`StateDataオブジェクトから有限要素データにアクセスする方法について見てゆく．有限要素データはNDSolve`FEM`FiniteElementDataオブジェクトに保存され，"FiniteElementData"の特性を通してアクセスされる．
{state} = NDSolve`ProcessEquations[{Laplacian[u[x, y], {x, y}] == 1, DirichletCondition[u[x, y] == 0, True]}, u, {x, 0, 1}, {y, 0, 1}]この問題のように，問題に時間依存がない場合には，NDSolve`StateDataオブジェクトの表示形式は，このことを"SteadyState"を表示することによって示す．非定常方程式については，積分が行われる時間区分が示される．
femdata = state["FiniteElementData"]解法プロセス中に生成されたデータオブジェクトはすべて，例えばFiniteElementDataのように，そのコンストラクタ関数として文書化されている．このチュートリアルでは，これらの関数の使用方法について説明する．
FiniteElementDataオブジェクトの表示形式は，離散化された偏微分方程式の自由度を示す．自由度は通常，メッシュ内のノード数に従属変数の数を掛けたものである．自由度は，組み立てられる行列の行と列の数を示す．
定常状態問題については，NDSolve`Iterateを呼び出すことで，LinearSolveを使って方程式系の解を求める．
NDSolve`Iterate[state];解はその後，解データオブジェクト内の状態オブジェクトに保存される．
{forwardSolution} = state["SolutionData"];
NDSolve`SolutionDataComponent[forwardSolution, "DependentVariables"]//Short
u /. NDSolve`ProcessSolutions[state]Plot3D[%[x, y], {x, 0, 1}, {y, 0, 1}]有限要素オプションをNDSolveに渡す
次のセクションでは，NDSolveの内部の有限要素関数について述べる．これらの関数にはオプションを含むものもある．内部の有限要素関数を使用せずに，これらのオプションをNDSolveに渡すことは常に可能である．これらは，NDSolveのMethodオプションで与えられる．
| NDSolve[…,Method{"PDEDiscretization"{"FiniteElement", femopts}}] | 定常状態問題について，オプション femopts とともに有限要素法を使う |
| NDSolve[…,Method{"PDEDiscretization"{"MethodOfLines", "SpatialDiscretization"->{"FiniteElement", femopts}}}] | 非定常問題の空間離散化について，オプション femopts とともに有限要素法を使う |
NDSolve[{Laplacian[u[x, y], {x, y}] == 1, DirichletCondition[u[x, y] == 0, True]}, u, {x, 0, 1}, {y, 0, 1}, Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.1}, "IntegrationOrder" -> 5}}]NDSolve[{D[u[t, x], t] - Laplacian[u[t, x], {x}] == t, DirichletCondition[u[t, x] == 0, True], u[0, x] == 0}, u, {t, 0, 1}, {x, 0, 1}, Method -> {
"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.1}, "IntegrationOrder" -> 5}}}]ToBoundaryMeshとToElementMeshのオプションは，上で示すように，"MeshOptions"を指定することによって，NDSolveに渡すことができる．さらに，InitializePDECoefficientsとPDESolveのオプションおよびそのサブオプションは，それぞれ"InitializePDECoefficientsOptions"と"PDESolveOptions"を通してNDSolveに渡すことができる．また，積分次数と補間次数はオプションを使って指定できる．内部的な有限要素法関数のオプションの用途と目的については，それぞれの参照ページおよびオプションについてのチュートリアル，有限要素のためのNDSolveオプションに記載されている．
ワークフローの概要
次のセクションでは，NDSolveが有限要素モデルを解く方法についてステップごとに示す．NDSolveは基本的にこれらのセクションで見られるように，同じ関数を使う．さらに，NDSolve`FEM`FiniteElementDataオブジェクトのコンテンツが下で紹介する関数を通して作成される．
初期化の段階中に，偏微分方程式と境界条件が分析され，別々の成分に分類され，これらの結果がそれぞれPDECoefficientDataとBoundaryConditionDataのオブジェクトに保存される．有限要素法データが設定され，FEMMethodDataオブジェクトに保存される．
| InitializePDECoefficients | 偏微分方程式の係数を初期化する |
| InitializeBoundaryConditions | 境界条件を初期化する |
| InitializePDEMethodData | 偏微分方程式のメソッドデータを初期化する |
第2段階では，偏微分方程式と境界条件が離散化され，DiscretizedPDEDataとDiscretizedBoundaryConditionDataに保存される．「離散化」は原則的に，偏微分方程式と境界条件がいわゆる系の行列によって表される離散バージョンによって近似されていることを意味する．
| DiscretizePDE | 初期化された偏微分方程式を離散化する |
| DiscretizeBoundaryConditions | 初期化された境界条件を離散化する |
最後の段階では，離散化された偏微分方程式と離散化された境界条件を組み合せて（DeployBoundaryConditions），LinearSolveを使って方程式を解く．
| DeployBoundaryConditions | 離散化された境界条件を離散化された偏微分方程式に配備する |
| LinearSolve | 方程式の線形系を解く |
いったん解が求まると，InterpolatingFunctionを構築することができる．
| ProcessPDESolutions | InterpolatingFunctionオブジェクトと生成する |
2つ目と3つ目の段階の代りに，関数PDESolveを利用することができる．
偏微分方程式問題の設定
偏微分方程式を有限要素法のような数値計算法によって解きやすい状態にするには，3つの成分が必要である．
以下のように，ここからの説明で使うモデル偏微分方程式を設定する．
op = Laplacian[u[x, y], {x, y}] - 1Ω = ImplicitRegion[0 ≤ x ≤ 1 && 0 ≤ y ≤ 1, {x, y}]Subscript[Γ, D] = DirichletCondition[u[x, y] == -1, x == 0]Subscript[Γ, N] = NeumannValue[-1 + u[x, y], x == 1]与えられた境界条件で領域上の偏微分方程式を解くには，NDSolveが直接使える．
NDSolve[{op == Subscript[Γ, N], Subscript[Γ, D]}, u, {x, y}∈Ω]これは，NDSolveのワンステップ機能を使う典型的なワークフローである．
NDSolve`FEM`関数を使ってステップごとに解く場合と比べるために，まずNDSolveで使われるデータを抽出し調べる．
{state} = NDSolve`ProcessEquations[{op == Subscript[Γ, N], Subscript[Γ, D]}, u, {x, y}∈Ω, Method -> {"PDEDiscretization" -> "FiniteElement"}]femdata = state["FiniteElementData"]femdata["Properties"]PDECoefficientDataは，InitializePDECoefficientsへの呼出しによって作成される．BoundaryConditionDataはInitializeBoundaryConditionsへの呼出しによって，FEMMethodDataはInitializePDEMethodDataへの呼出しによって，作成される．これらのデータオブジェクトが続く計算のデータを保つ．これらのデータオブジェクトの作成は，次のセクションで説明される．
定常偏微分方程式
初期化段階
初期化段階では，NDSolve`ProcessEquationsの実行中に作成されるすべてのデータオブジェクトが作成される．このセクションは，そのプロセスについて示す．このセクションを飛ばして，離散化段階に進み，NDSolve`ProcessEquationsが作成する初期化されたデータ構造を利用することも可能である．これで，例えば新しい数値離散化法に，NDSolve`ProcessEquationsを方程式プリプロセッサとして使用することができる．
変数データと解データ
NDSolveが解を求めるには，従属変数(u)と独立変数({x,y})，および領域(Ω)を指定しなければならない．このことは，有限要素関数でも同じである．変数データと解データは，別々の成分に対応する特定のデータ構造を作成するリストに保存される．これらのデータリストは，NDSolve`VariableDataとNDSolve`SolutionDataを使って生成することもできる．データリストの要素は，NDSolve`SetSolutionDataComponentを使って挿入または変更することができる．
vd = NDSolve`VariableData[{"DependentVariables" -> {u}, "Space" -> {x, y}}]変数データに加わるのが解のデータ構造である．解のデータ構造は，基本的に，変数データが数値的に生まれ変わったものである．例えば，解データの"Space"成分の数値バージョンは，ElementMeshも含むNumericalRegionである．「要素メッシュの作成」に領域Ωからのメッシュの作成についての詳細が記載されている．
nr = ToNumericalRegion[Ω]sd = NDSolve`SolutionData[{"Space" -> nr}]InitializePDECoefficients
偏微分方程式は
で定義される．係数 
，
，
，
はスカラーである．
，
，
は，長さ 
のベクトルであり，
は 
×
行列である．InitializePDECoefficientsへの入力係数は，上の偏微分方程式とモデル化される偏微分方程式を対応させることによって求めることができる．
opcoefficients = {"DiffusionCoefficients" -> {{IdentityMatrix[2]}}, "LoadCoefficients" -> {{1}}}明示的に指定されていない係数はすべてゼロであると想定される．
initCoeffs = InitializePDECoefficients[vd, sd, coefficients]初期化された係数は，PDECoefficientDataのデータ構造に保存される．表示形式は，初期化された係数の系の大きさ，つまり従属変数の数(1)，および演算子の大きさ，つまり空間次元(2)を示す．
{initCoeffs["SystemSize"],
initCoeffs["SpatialDimension"]}initCoeffs["All"]初期化注に，係数は定常係数，非定常係数等，さまざまな係数のカテゴリに分類される．このことについては詳しく後で説明する．
initCoeffs["Properties"]InitializeBoundaryConditions
偏微分方程式の初期化と同じように，境界条件も初期化されなければならない．
initBCs = InitializeBoundaryConditions[vd, sd, {{Subscript[Γ, D], Subscript[Γ, N]}}]境界条件の係数は，偏微分方程式の係数と同じようなカテゴリに分類される．
initBCs["Properties"]InitializePDEMethodData
後に続くステップで必要となる，有限要素法データを初期化するには，InitializePDEMethodDataを使う．現在，このフレームワークで使用できる離散化法は，有限要素法だけである．したがって，デフォルトで，InitializePDEMethodDataはFEMMethodDataオブジェクトを生成する．初期化中にメソッドに特定の設定が行われる．例えば，メソッドの補間と積分の次数が他の諸々と一緒に設定され，ElementMeshがNumericalRegionから生成される．
methodData = InitializePDEMethodData[vd, sd]初期化された偏微分方程式のメソッドデータは，FEMMethodDataオブジェクトに保存される．FEMMethodDataの表示形式は，自由度，各従属変数の補間次数，および積分次数を示す．
{methodData["DegreesOfFreedom"],
methodData["InterpolationOrder"],
methodData["IntegrationOrder"]}「要素メッシュの生成」のチュートリアルに説明されているように，メッシュ生成オブションはすべて指定することができる．
methodData = InitializePDEMethodData[vd, sd, Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.1}}]{methodData["DegreesOfFreedom"],
methodData["InterpolationOrder"],
methodData["IntegrationOrder"]}メソッドの初期化中にElementMeshオブジェクトが生成され，NumericalRegionオブジェクトの中に保存される．
mesh = nr["ElementMesh"]方程式系の自由度は，複数の従属変数が与えられている場合には，要素メッシュのノード，従属変数の数，従属変数の補間次数を組み合せたものから求められる．従属変数が1つだけの場合には，自由度は，ElementMeshオブジェクトのノードに対応する．
Length[mesh["Coordinates"]] === methodData["DegreesOfFreedom"]"IntegrationOrder"は，有限要素の演算子を積分するのに使われる確度の次数である．
InitializePDEMethodData[vd, sd, Method -> {"FiniteElement", "IntegrationOrder" -> 5}]["IntegrationOrder"]一次要素メッシュについては，デフォルトの積分次数は，二次の正確な積分点を使用するが，二次要素メッシュについては，デフォルトでは四次の正確な積分点を使用する．サポートされる最大の積分次数は5である．
積分次数を減らすと，解が速く解かれるが，この場合には，あまり正確ではない解になる可能性がある．積分次数を増やすと，有限要素モデルを解くのにより長い時間がかかる．
{s1} = NDSolve`ProcessEquations[{op == Subscript[Γ, N], Subscript[Γ, D]}, u, {x, y}∈Ω, Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.1}}];s1["FiniteElementData"]["FEMMethodData"] === methodData初期化された有限要素データをNDSolve`StateDataから抽出する
有限要素法のデータを手作業で初期化する代りに，偏微分方程式，境界条件，NDSolve`ProcessEquationsをこのために利用することができる．
{state} = NDSolve`ProcessEquations[{op == Subscript[Γ, N], Subscript[Γ, D]}, u, {x, y}∈nr, Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.1}}];
femdata = state["FiniteElementData"]femdata["FEMMethodData"] === methodDatafemdata["PDECoefficientData"] === initCoeffs離散化段階
離散化段階では，連続偏微分方程式と境界条件は，離散変数で近似される．このプロセスは，偏微分方程式を解くプロセスの中心となるものである．
DiscretizePDE
discretePDE = DiscretizePDE[initCoeffs, methodData, sd]DiscretizedPDEDataの表示形式は，系の行列の一次元を示す．
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"]MatrixPlot[stiffness]{discretePDE["LoadVector"], discretePDE["StiffnessMatrix"], discretePDE["DampingMatrix"], discretePDE["MassMatrix"]}行列の組立て過程において，有限要素にアクセスすることは可能である．有限要素の形式は，行列の組立て過程が効率的でなるものであることに注意する．
finiteElements = DiscretizePDE[initCoeffs, methodData, sd, "SaveFiniteElements" -> True];
Short[finiteElements["StiffnessElements"], 3]DiscretizeBoundaryCondition
DiscretizeBoundaryConditionsは，与えられた境界条件の離散化バージョンを計算する．後のステップで，離散化された境界条件が実際の系の行列に配備される．
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd]DirichletConditionとNeumannValueの境界指定は，違った形で表示される．NeumannValueからの寄与は，系の負荷行列と系の剛性行列に対してである．
{discreteBCs["LoadVector"],
discreteBCs["StiffnessMatrix"]}離散化された境界条件の"StiffnessMatrix"成分は，NeumannValueが従属変数に依存する場合には非零の寄与を持ち，それ以外の場合にはゼロである．言い換えれば，ノイマン値が一般化されたノイマン値である場合には，"StiffnessMatrix"成分への寄与がある．
DirichletCondition表示は，位置を含む行列を通して与えられる．これらの位置は，条件にどの自由度が適用されるべきかを指定する．返される2つ目の行列は，適用される値を含む．
{discreteBCs["DirichletMatrix"],
discreteBCs["DirichletValues"]}解の段階
DeployBoundaryConditions
離散化された境界条件が効果を持つためには，系の行列に配備されなければならない．メモリを節約するために，これはその場で行われ，新しい行列は生成されない．
DeployBoundaryConditions[{load, stiffness}, discreteBCs]DeployBoundaryConditionsは属性HoldFirstを持つので，系の行列はコピーをまず作成するということをせずに，その場で変更されることがある．このことはしかし，DeployBoundaryConditionsへの入力が，実際の行列としてではなく値として行列を持つシンボルで与えられなければならないことを意味する．変更されていない系の行列が後で必要となる場合には，DiscretizedPDEDataから再度抽出することができる．言い換えれば，変更されていない系の行列がもはや不要であり，コンピュータのメモリ消費が問題となるような場合には，DiscretizedPDEDataを例えばClearAllを使ってクリアする必要がある．
境界条件の配備中に，NeumannValueからの寄与が剛性行列と負荷行列に加えられる．その後，系の行列は，DirichletConditionが満たされるように変更される．
{load, stiffness}一般に，境界条件の配備の後でのみ，方程式の系を解くことができるということに注意する．
方程式を解く
方程式の系は，LinearSolveを使って解くことが可能である．
Short[solution = LinearSolve[stiffness, load]]後処理
後処理のステップとして，補間関数をLinearSolveから得た解から作成することができる．この解のデータは求めた解でアップデートされ，ProcessPDESolutionsがその後InterpolatingFunctionを構築する．
NDSolve`SetSolutionDataComponent[sd, "DependentVariables", Flatten[solution]];{ifun} = ProcessPDESolutions[methodData, sd]ProcessPDESolutionsは，複数の従属変数についても使えるという点が利点である．この場合には，ProcessPDESolutionsを使う代りに，ElementMeshInterpolationを使って補間関数を直接作成することができる．ここで第1引数はメッシュ，第2引数はメッシュノードにおける値である．
ifun = ElementMeshInterpolation[mesh, solution]Plot3D[ifun[x, y], {x, y}∈mesh]通常，補間関数はその領域外でクエリが行われた場合には，extrapolate外挿する．この動作を変更することは可能である．
ifun[2, 2]
ifun = ElementMeshInterpolation[{mesh}, solution, "ExtrapolationHandler" -> {Function[Indeterminate], "WarningMessage" -> False}]ifun[2, 2]"DefaultExtrapolationHandler"のシステムオプションを設定すると，デフォルトの動作を変更することができる．
SystemOptions["FiniteElementOptions"]補間関数からElementMeshを抽出することも可能である．
ifun["ElementMesh"]非線形の偏微分方程式
以下のセクションでは，非線形の定常偏微分方程式の解について話し合う．これは複数のレベルで行う．まず，トップレベルの例について紹介し，これをどのように使うのか，どのように実装するのかについてさらに深く見ていくことにする．
トップレベルの例
非線形の偏微分方程式を解く過程を説明するために，特定の境界を持つ領域上での最小曲面を求める必要がある．この最小曲面は，次の非線形方程式を解くことによって求めることができる．
係数
は，従属変数 
を含み，特に従属変数 
の勾配を含むので，非線形である．非線形の偏微分方程式を公式化するためには，方程式1に詳しく示されている
の任意の係数は，
と 
および 
のような一次偏微分係数の両方に従属できる．しかしこれらは， 
のような一次より高次の微分係数に従属することはできない．
一般に，解法として有限要素法を宇買う非線形方程式は，NDSolveを普通に使う場合のように設定することができる． NDSolveは，その引数をホールドしないので，入力方程式は，NDSolveに到達する前に評価される．必要以上に早く評価されることを避けるために，偏微分方程式はInactiveの形式で指定することができる．Inactiveでラップされていない非線形方程式は，好ましくない不要に高次の導関数を生み出すことがあり，そのような導関数はNDSolveで解くことはできない．特に，この場合のように，従属変数の導関数が拡散係数 
内に存在する場合には解けない．そのような場合には，偏微分方程式を非アクティブ化しなければならない．偏微分方程式の非アクティブ化は，偏微分方程式を指定する最も一般的な方法であり，どのメソッドを使っていいのか分からない場合には，このメソッドを選ぶべきである．
c = 1 / Sqrt[(1 + Grad[u[x, y], {x, y}].Grad[u[x, y], {x, y}])];
eqn = Inactive[Div][{{c, 0}, {0, c}}. Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0;Ω = Disk[];
bcs = DirichletCondition[u[x, y] == Sin[2π * (x + y)], True];FindRootによく似て，NDSolveは定常非線形偏微分方程式を有限要素法で解く場合には，それぞれの従属変数に対して，初期シード値が必要である．初期シード値は，InitialSeedingで指定することができる．
usol = NDSolveValue[{eqn, bcs}, u, {x, y}∈Ω, InitialSeeding -> {u[x, y] == 0}]初期シード値が指定されていない場合には，ゼロであると見なされる．
Plot3D[usol[x, y], {x, y}∈Ω, PlotRange -> All]ここで方程式が非アクティブな形式ではない場合には，NDSolveはこの偏微分方程式を解くことができない．
NDSolveValue[{eqn//Activate, bcs}, u, {x, y}∈Ω, InitialSeeding -> {u[x, y] == 0}];
非アクティブな偏微分方程式についての詳細は，「形式的な偏微分方程式」についてのセクションを参照されたい．
以下のセクションでは，上のような非線形方程式が与えられたときにNDSolveが内部で行う解の一般的な過程を紹介する．線形方程式の場合と同じ段階が非線形方程式の場合にも使われる．
初期化段階
非線形方程式の初期化段階は，線形の場合とほとんど同じである．非線形の場合に唯一異なるのは，FindRootの場合と同じように，初期シード値が指定されなければならない点である．これは，解データのコンポーネント"DependentVariables"を初期シードに設定することで行う．
nr = ToNumericalRegion[Ω];
vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData[{"DependentVariables", "Space"} -> {{0}, nr}];この場合，初期シード値は，数値でも独立変数である 
と 
に基づく任意式でもよい．かなり非線形な問題の場合には，あまり非線形ではない問題を解いて，結果として得られるInterpolatingFunctionをより高度な非線形問題の初期シード値として使うという方法がある．初期シード値が指定されていない場合には，
が初期シード値として使われる．
methodData = InitializePDEMethodData[vd, sd];
ibcs = InitializeBoundaryConditions[vd, sd, {{bcs}}];
pdeCoeffs = InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{{{c, 0}, {0, c}}}}];離散化と解の段階
非線形方程式を解くことは，インタラクティブな過程である．つまり，最下層レベルでは，インタラクティブな過程にFindRootが使われる．まず最初のステップは，偏微分方程式係数を線形化することである．この過程の詳細は下で説明する．その後初期シード値を使って，線形化された偏微分方程式を評価する．これでLinearSolveで解ける線形方程式系が作成される．この過程の結果は，偏微分方程式に代入し直され，アップデートされた，線形化された偏微分方程式は解の新しい近似値で離散化される．この過程は，理想的には解が収束するまで繰り返される．
関数PDESolveは FindRootの呼出しを調整し，線形化と繰返しの離散化を行い，解を含む解データオブジェクトを返すか，アルゴリズムが収束に失敗した場合にはPDESolveが求めることができた解の最良の近似値を返すかする．
sdNew = PDESolve[pdeCoeffs, ibcs, vd, sd, methodData];後処理
InterpolatingFunctionオブジェクトは，解データから作成できる．
{eif} = ProcessPDESolutions[methodData, sdNew]Plot3D[eif[x, y], {x, y}∈Disk[], PlotRange -> All]線形化過程
FindRootはデフォルトでニュートン・ラフソン(Newton–Raphson)法を使って非線形方程式を解く．このメソッドは， 想定される解，初期シード値で始め，方程式の線形化された形式でこれを改善しようとする．
線形化過程をよりよく理解するためには，スカラー関数のためのニュートン法の微分を詳しく見ると役に立つ．
を 
の非線形関数とし，
を以下が成り立つ未知の解とする．
初期シード値 
を取り，それを未知の解 
から引くと，結果として剰余 
が得られる．
g[x_] := Block[{y1, y2}, {y1, y2} = x;{y1 + 2 * y2 - 2, y1 ^ 2 + 4 * y2 ^ 2 - 4}]
D[g[{y1, y2}], {{y1, y2}}]jacobian[x_] := Block[{y1, y2}, {y1, y2} = x;{{1., 2.}, {2 y1, 8 y2}}]
x0 = {1., 2.};
x1 = x0 - LinearSolve[jacobian[x0], g[x0]]
x2 = x1 - LinearSolve[jacobian[x1], g[x1]]解が求まるまで，同じ過程を続ける．次のステップとして，反復過程をベクトルで与えられるステップを計算する関数で囲む．2つの連続するステップのノルムが閾値より大きい間に，新しい次のステップを計算する．
step[x_] := x - LinearSolve[jacobian[x], g[x]]result = NestWhile[step, x0, Norm[#1 - #2] > 10 ^ -3&, 2, 10]g[result]与えられた例で1点気を付けなければならないことは，
と関数行列式 
の両方が一定の部分 
に従属する部分があることである．有限要素法を使って非線形偏微分方程式を解く場合には，反復過程の前に一定部分を一度計算し，それらを各ステップで計算される非一定部分に加えるとよい．
方程式1で得た係数形式の非線形偏微分方程式の初期化過程も同じステップで行われる．方程式2から始めて，
と 
を設定して時間独立係数形式を得る．
に変える．これで時間独立係数形式は以下のように書き換えられる．
線形化を行うためには，上のスカラーの場合と同じように行う．
および 
は 
の関数であり，以下が得られる未知の解である．
任意のシード値 
を取り，これを未知の解 
から引くと，以下のように剰余 
が得られる．
線形化された偏微分方程式の実際の解は，方程式3で与えられた係数形式を，任意数の従属変数の以下の係数と一緒に利用する．
ここで 
は整数である．
は
から 
（従属変数の数）までである．
の 
個のスカラー方程式に対して 
が得られる．
ニュートンの求根法とそれをどのように偏微分方程式に適用するかについての説明の後は，方程式
をより深いレベルで解くことを続けて行い，「離散化と解の段階」 のセクションで紹介したPDESolveへの呼出しがどのように使えるかを示す．最初のステップでは，線形化された係数が構築される．
depVars = NDSolve`SolutionDataComponent[vd, "DependentVariable"]vars = NDSolve`SolutionDataComponent[vd, "Space"]{s = Length[depVars], n = Length[vars]}dep = (#@@vars)& /@ depVarsddep = Grad[dep, vars]load = pdeCoeffs["LoadCoefficients"][[All, 1]];
dload = pdeCoeffs["LoadDerivativeCoefficients"][[All, 1, All, 1]];
diff = pdeCoeffs["DiffusionCoefficients"];
cconv = pdeCoeffs["ConservativeConvectionCoefficients"][[All, All, All, 1]];
conv = pdeCoeffs["ConvectionCoefficients"][[All, All, 1]];
react = pdeCoeffs["ReactionCoefficients"];
F = load - Plus@@@Table[conv[[r, c]].ddep[[c]], {r, s}, {c, s}] - react.dep
Γ = -Plus@@@Table[diff[[r, c]].ddep[[c]], {r, s}, {c, s}] - Plus@@@Table[cconv[[r, c]] * dep[[c]], {r, s}, {c, s}] + dload[[All, 1]]方程式6の
についての偏微分係数 
は，Dを使って計算することができる．
-Table[D[ Γ[[i, k]], ddep[[j, l]]], {i, s}, {j, s}, {k, n}, {l, n}]//Simplify//MatrixForm方程式7の
についての偏微分係数 
は，Dを使って計算することができる．
-Table[D[ Γ[[i, k]], dep[[j]]], {i, s}, {j, s}, {k, n}]方程式 8 の 
についての偏微分係数 
は，Dを使って計算することができる．
-Table[D[ F[[i]], ddep[[j, l]]], {i, s}, {j, s}, {l, n}]方程式9の 
についての偏微分係数 
は，Dを使って計算することができる．
-Table[D[ F[[i]], dep[[j]]], {i, s}, {j, s}]これらの係数は，ユーティリティ関数LinearizePDECoefficientsを使って，便利なことに有限要素のコンテキストから求めるとこができる．
linearizedPDECoeffs = LinearizePDECoefficients[pdeCoeffs, vd, sd]linearizedPDECoeffs["DiffusionCoefficients"]//Simplify//MatrixFormlinearizedPDECoeffs["LoadDerivativeCoefficients"]//Simplify//MatrixForm線形化された偏微分方程式を解く
一旦偏微分方程式が線形化されると，係数は分割される．これは，求根アルゴリズムの異なる部分で係数が必要となるからである．係数形式
左辺のすべての項は，剛性行列にまとめられる．剛性行列は偏微分方程式の関数行列式である．剛性行列は，上のスカラーの例からの 
の項に類似するものである．
は 
の接線であるので，剛性行列も接線剛性行列と呼ばれる．偏微分方程式の右辺の項は， 
に対応する．PDECoefficientDataの分割は，関数SplitPDECoefficientsで行う．
{linLoadPDEC, linStiffnessPDEC, linDampingPDEC, linMassPDEC} = SplitPDECoefficients[linearizedPDECoeffs, vd, sd]係数 
と 
はこれで変数linLoadPDEC内にあるが，残りはlinStiffnessPDEC内にある．減衰行列と質量行列の係数は，ここで扱っている問題については存在しないので，linDampingPDECおよびlinMassPDECは空である．
sdU = EvaluateInitialSeeding[methodData, vd, sd];seed = NDSolve`SolutionDataComponent[sdU, "DependentVariables"];linear = DiscretizePDE[linearizedPDECoeffs, methodData, sdU];{linearLoad, linearStiffness, linearDamping, linearMass} = linear["SystemMatrices"];linearBCs = DiscretizeBoundaryConditions[ibcs, methodData, sdU]femRHS[u_ ? VectorQ] := Block[{load, nonlinear, nonlinearLoad, nonlinearBCs},
NDSolve`SetSolutionDataComponent[sdU, "DependentVariables", u];
nonlinear = DiscretizePDE[linLoadPDEC, methodData, sdU, "Nonlinear"];
nonlinearLoad = nonlinear["LoadVector"];
nonlinear = Null;
load = linearLoad + nonlinearLoad;
nonlinearLoad = Null;
(* subtract the linear Robin boundary value *)
load -= linearBCs["StiffnessMatrix"].u;
nonlinearBCs = DiscretizeBoundaryConditions[ibcs, methodData, sdU, "Nonlinear"];
DeployPartialBoundaryConditions[{load, Null}, nonlinearBCs];
DeployPartialBoundaryConditions[{load, Null}, linearBCs];
load = -load;
Normal[Flatten[load]]
];ここで与えられた関数femRHSは，ここで取り扱っている例には存在しないが，一般化されたノイマン値も取り扱う．実際，ここで与えられた関数は，実際の実装に大変近いものである．
関数femRHSでは負荷ベクトルのみが利用されることに注意する．これは，SplitPDECoefficientsコマンドで生成されるlinLoadPDECには，負荷係数 
と負荷微分係数 
のみが含まれるからである．
一旦nonlinearおよびnonlinearLoadの変数が必要なくなると，これらはNullに設定され，これらが指すメモリが解放される．
femJacobian[u_ ? VectorQ] := Block[{stiffness, nonlinear, nonlinearStiffness, nonlinearBCs}, NDSolve`SetSolutionDataComponent[sdU, "DependentVariables", u];
nonlinear = DiscretizePDE[linStiffnessPDEC, methodData, sdU, "Nonlinear"];
nonlinearStiffness = nonlinear["StiffnessMatrix"];
nonlinear = Null;
stiffness = linearStiffness + nonlinearStiffness;
nonlinearStiffness = Null;
nonlinearBCs = DiscretizeBoundaryConditions[ibcs, methodData, sdU, "Nonlinear"];
DeployPartialBoundaryConditions[{Null, stiffness}, nonlinearBCs];
DeployPartialBoundaryConditions[{Null, stiffness}, linearBCs];
stiffness
];ここでも，関数femJacobianは実際の実装に大変近いので，もはや必要がない変数はNullに設定される．
関数行列式の構築には，線形化された剛性成分が使われる．これには，拡散係数 
，保守対流係数 
，対流係数 
，反応係数 
が含まれる．これらはすべてSplitPDECoefficientsコマンドで生成されたlinStiffnessPDEC PDECoefficientDataに保存される．
実装では，線形である係数はすべて一度だけ集められることに注目されたい．関数femRHSおよびfemJacobianのどちらでも，非線形係数のみが計算され，先に計算された線形成分に加えられる．これで各反復において変化しない線形係数を無駄に計算せずに済む．最後の非線形と線形の境界条件が配備される．
DeployDirichletConditions[seed, linearBCs]最後に，線形化された偏微分方程式の解を求める．このためには，アフィン共変ニュートン法が使われる．
root = U /. FindRoot[femRHS[U], {U, seed}, Jacobian -> femJacobian[U], Method -> {"AffineCovariantNewton"}];NDSolve`SetSolutionDataComponent[sdU, "DependentVariables", root];InterpolatingFunctionオブジェクトが解データから作成できる．
{eif} = ProcessPDESolutions[methodData, sdU]Plot3D[eif[x, y], {x, y}∈Disk[], PlotRange -> All]前に説明した手順は，PeriodicBoundaryConditionが存在する場合にはうまく作動しないという点において制約されている．
非定常偏微分方程式
非定常偏微分方程式は，時間依存の偏微分方程式である．最も簡単な非定常偏微分方程式は，偏微分方程式の一部として，依存変数の時間導関数を持つ．この場合は，離散化された時間導関数を表す追加の系の行列が作成される．その後，系の行列の時間積分を，例えばNDSolveを使って行うことができる．それ以外の点では，ワークフローは定常偏微分方程式の場合とほぼ同じである．より複雑な例には，依存する係数および／あるいは境界条件が含まれる．後ろのセクションでそれらを時間積分する方法について述べる．
定常係数と定常境界条件を持つ非定常偏微分方程式 — はじめに
最初の例として，1Dの熱方程式を見てみよう．空間領域が0から1で，時間領域が0から1である，次のモデル偏微分方程式を考えてみる．
ufun = NDSolveValue[{D[u[t, x], {t, 1}] - D[u[t, x], {x, 2}] == 1, u[0, x] == 0, DirichletCondition[u[t, x] == 0, True]}, u, {t, 0, 1}, {x, 0, 1}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}}]この場合，メソッドオプションは，空間変数 
が有限要素法で離散化されるべきであり，時間積分が線分法で行われることを指定する．このアプローチは，半離散化と呼ばれる．反対に，完全な離散化は，時間変数 
と空間変数の両方を離散化する．
Plot3D[ufun[t, x], {t, 0, 1}, {x, 0, 1}]偏微分方程式の解を求める過程は，定常の場合に似ている．その大きな理由として，偏微分方程式の係数も境界条件も時間に依存しないことが挙げられる．非定常系の解を求めるために，偏微分方程式を半離散化する方法が取られる．偏微分方程式は空間(
)内で，有限要素法を使って離散化される．そうすると，その有限要素法の離散化の結果は，時間変数(
)における常微分方程式系である．NDSolveを使って，その常微分方程式系が時間積分される．
vd = NDSolve`VariableData[{"DependentVariables" -> {u}, "Space" -> {x}, "Time" -> t}]nr = ToNumericalRegion[FullRegion[1], {{0, 1}}]sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}]initCoeffs = InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-IdentityMatrix[1]}}, "LoadCoefficients" -> {{1}}, "DampingCoefficients" -> {{1}}]initBCs = InitializeBoundaryConditions[vd, sd, {{DirichletCondition[u[t, x] == 0, True]}}]methodData = InitializePDEMethodData[vd, sd]mesh = nr["ElementMesh"]次に離散化を行う．非定常系が考慮されているのだが，偏微分方程式の係数も境界条件も時間に依存しないので，離散化は定常成分のみを持つ．
discretePDE = DiscretizePDE[initCoeffs, methodData, sd]discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd]{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs]次に初期条件が設定され，NDSolveは離散化された方程式系の時間積分を行う．
init = Table[{0.}, {methodData["DegreesOfFreedom"]}];
init[[discreteBCs["DirichletRows"]]] = discreteBCs["DirichletValues"];tufun = NDSolveValue[{damping.u'[ t] + stiffness.u[ t] == load, u[0] == init}, u, {t, 0, 1}, Method -> {"TimeIntegration" -> "IDA"}, AccuracyGoal -> $MachinePrecision / 4, PrecisionGoal -> $MachinePrecision / 4]下の方のセクションで，関数行列式を指定することによって時間積分の過程を速める方法について説明する．
NDSolveが解の過程のすべての面を制御していた上のNDSolveへの呼出しから自動的に選択されたオプションにマッチするように，NDSolveオプションを選ぶ．
ClearAll[fun]
fun[t_ ? NumericQ] := fun[t] = ElementMeshInterpolation[{mesh}, tufun[t]]Plot3D[fun[t][x] - ufun[t, x], {t, 0, 1}, {x, 0, 1}]定常係数と定常境界条件を持つ非定常偏微分方程式
境界条件は，
では 
であり，
では 
である．マーカーについての詳細は，「要素メッシュの生成」のチュートリアルを参照されたい．
偏微分方程式は，一階の時間導関数を持つが，偏微分方程式の係数 
と境界条件は時間依存ではない．偏微分方程式を解く1つの方法は，まず偏微分方程式と境界条件を離散化してから，結果の常微分方程式を時間積分器としてNDSolveで時間積分する方法である．この過程は，NDSolveが内部で非定常偏微分方程式を解く方法であり，これについては次のセクションで見る．
この例では，マーカー付きの境界メッシュが事前定義されている．形状は，Scanning Electrochemical Microscopy (SECM)プローブをモデル化する．このプローブは，7つの電極が円状のグラス構造に埋め込まれているので，七極管である．
bmesh = [image];boundaryMesh = ToBoundaryMesh[bmesh]BoundaryMeshRegionオブジェクトとElementMeshオブジェクトの間の関係についての詳細は， 「要素メッシュの生成」チュートリアルを参照のこと．
事前定義された境界メッシュ領域では，整数マーカーが境界の表面と点に設定されている．これらのマーカーは，境界条件が解のプロセス中に適用される位置として可視化し，使うことができる．
Show[
bmesh,
boundaryMesh["Wireframe"["MeshElement" -> "PointElements", "MeshElementStyle" -> {Directive[White], Directive[PointSize[0.015], Red], Directive[PointSize[0.02], Blue]}]]]この境界条件はディリクレ条件だけを含むので，点の要素とマーカーだけが重要である．ディリクレ条件は常に点要素に作用する．境界メッシュも実装された境界要素についてマーカーを持つ．これらは，ノイマンタイプの境界値を指定するのに必要であり，これらの値は常に境界要素に作用する．
Show[
bmesh,
boundaryMesh["Wireframe"["MeshElement" -> "BoundaryElements", "MeshElementStyle" -> {Directive[EdgeForm[], FaceForm[]], Directive[EdgeForm[Black], FaceForm[Red]],
Directive[EdgeForm[Black], FaceForm[Blue]]}]]]mesh = ToElementMesh[boundaryMesh, "NodeReordering" -> True, "MeshOrder" -> 1]厳密に言うと，ノードの並び替えはオプションである．ノードの並び替えメソッドのオプションがTrueに設定されている場合には，系の行列が持つバンド幅は最小化される．低いバンド幅は，反復ソルバに有益である．ノードの並び替えは1度行うだけでよく，ソルバが各時間ステップに使われるので，非定常問題にノードの並び替えを使うことが役に立つ場合が多い．
vd = NDSolve`VariableData[{"DependentVariables", "Space", "Time"} -> {{u}, {x, y, z}, t}]この場合，ElementMeshオブジェクトがNumericalRegionオブジェクトを作成するために直接使われたかもしれない ImplicitRegionの指定は存在しない．
sd = NDSolve`SolutionData[{"Space", "Time"} -> {ToNumericalRegion[mesh], 0.}]解のデータを追加する．"Time"成分は，初期時間t0であり，積分が始まる．"DependentVariables"と"TemporalDerivatives"の成分を設定すると，初期条件と初期条件の導関数が設定される．
pdeCoefficients = InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-IdentityMatrix[3]}}, "DampingCoefficients" -> {{1}}]bcs = InitializeBoundaryConditions[vd, sd, {{DirichletCondition[u[t, x, y, z] == 1, ElementMarker == 1], DirichletCondition[u[t, x, y, z] == 0, ElementMarker == 2]}}]methodData = InitializePDEMethodData[vd, sd]dpde = DiscretizePDE[pdeCoefficients, methodData, sd]dbcs = DiscretizeBoundaryConditions[bcs, methodData, sd]{l, s, d, m} = dpde["SystemMatrices"]DeployBoundaryConditions[{l, s, d}, dbcs]最初の非定常行列である減衰行列が追加されたことに注意する．これは後で時間積分を行う際に必要となる．
MatrixPlot[s]剛性行列の構造に注意する．"NodeReordering"->Falseをメッシュ生成の過程に使うと，あまり構造化されていない系の行列が返される．
stationarySolution = LinearSolve[s, l];stationarySolution = LinearSolve[s, l, "Method" -> "Pardiso"];ElementMeshSurfacePlot3D[stationarySolution, mesh, Boxed -> False]eufun = ElementMeshInterpolation[{mesh}, stationarySolution, "ExtrapolationHandler" -> {Function[Indeterminate], "WarningMessage" -> False}]Manipulate[ContourPlot[eufun[x, y, z], {y, -12, 12}, {z, -11, 0}, ColorFunction -> "TemperatureMap", Contours -> Range[0, 1, 0.2], PlotRange -> All], {{x, 0}, -12, 12}
, SaveDefinitions -> True]Show[Graphics3D[{Opacity[0.25], EdgeForm[], ElementMeshToGraphicsComplex[mesh, All, "MeshElement" -> "BoundaryElements"]}],
ContourPlot3D[eufun[x, y, z], {x, -12, 12}, {y, -12, 12}, {z, -11, -1}, Mesh -> None, Axes -> None], Boxed -> False]偏微分方程式の離散化が常微分方程式の系で起る．常微分方程式の時間積分は，NDSolveで行われる．これには，境界条件にマッチする初期条件を設定する必要がある．
init = Table[{0.}, {methodData["DegreesOfFreedom"]}];
init[[dbcs["DirichletRows"]]] = dbcs["DirichletValues"];ElementMeshSurfacePlot3D[init, mesh, Boxed -> False]NDSolveの時間積分の速度を上げるために，オプションを使って関数行列式の疎なパターンを指定することが可能である．系が全体として’NDSolveを通して解かれる場合には，これは不要である．ここでは，NDSolveがベクトル変数の構造を読み解こうとはしないため，必要である． 関数行列式の疎なパターンは，減衰行列のパターンを通して使用可能である．
{dampingSparsity = d["PatternArray"],
stiffnessSparsity = s["PatternArray"]}Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
ufun = NDSolveValue[{
d.u'[ t] + s.u[ t] == l
, u[ 0] == init}, u, {t, 0, 1}
, Method -> {"EquationSimplification" -> "Residual"}
, Jacobian -> {"X" -> {Automatic, Sparse -> stiffnessSparsity}, "XP" -> {Automatic, Sparse -> dampingSparsity}}
, EvaluationMonitor :> (currentTime = t;)
]]前のバージョンのNDSolveでは，関数行列式が"Jacobian{Automatic,Sparsesparsity}"と指定されていた．"sparsity" は，減衰行列の疎なパターンであった．バージョン12.2以降では，剛性行列と減衰行列の両方に対して，疎なパターンの指定が優先されるようになった．古い関数行列式の設定も引き続き使用可能であるが，時間積分がゆっくりになるため，剛性行列と減衰行列の両方に対して関数行列式を指定することに切り替えることを強く推奨する．
EvaluationMonitorを使うと，解の過程が少しゆっくりになるが，解の過程をライブで管理する事ができる．
解は，非定常の系がどのように上で計算された定常状態に到達するかを示す．
allGraphics =
ElementMeshSurfacePlot3D[ufun[#], mesh, Boxed -> False, Background -> Gray, ColorFunctionScaling -> False]& /@ {0.01, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1.};GraphicsGrid[Partition[allGraphics, 4], ImageSize -> Large]時間についての単一の補間関数と空間変数を構築することができる．このため，非定常補間が設定された時点のステップが抽出されなければならない．
Short[iftime = ufun["Coordinates"][[1]]]次に，補間の値を設定する．それぞれの時間のステップについて，ノードにおける値のリストが必要である．
data = Transpose[ufun["ValuesOnGrid"], {1, 3, 2}];efun = ElementMeshInterpolation[{iftime, mesh}, data]Plot[efun[t, 0., 0., -5], {t, 0, 1}]定常係数と定常境界条件を持つ非定常偏微分方程式のモデル次数の簡約
系の行列へのアクセスが重要である非定常の場合，モデル次数の簡約を考慮する．モデル次数の簡約が問題とするのは，非定常方程式系において、もとの特性のほとんどを保ちつつ，もとのものより小さい同等の方程式系が存在するかどうかということである．
減衰行列 
および剛性行列 
は 
×
の大きさで，負荷行列は 
×1の大きさである．今度は，
の自由度のずっと小さい2つ目の方程式系を設定する：
2つ目の自由度のため，簡約された系 
は，もとの系 
よりもずっと小さく，時間積分はかなり効率的である．以下のアルゴリズムは，
と 
の間でマップする射影ベクトル 
を設定する：
ここで 
の大きさは 
，
の大きさは 
，そして 
は近似の相違である．
射影ベクトル 
を簡単に説明すると，非零の開始ベクトル 
と行列 
を与えられた場合に，アルゴリズムは以下の式のようになる正規直交 
を求める：
mor[{load_, stiffness_, damping_}, size_Integer] := Module[
{lsf, tmp, vt, v},
lsf = LinearSolve[-stiffness];
vt = Orthogonalize[Transpose[-lsf[damping.load]]];
Do[tmp = lsf[damping.vt[[-1]]];
vt = Orthogonalize[Join[vt, {tmp}]];
, {i, 2, size}];
v = Transpose[vt];
{{vt.load, vt.stiffness.v, vt.damping.v}, v}
]余談だが，プラットフォームによっては，Method"Pardiso"をLinearSolveのオプションとして使うことができるものもある．このソルバは非常に効率的であり，非正則行列に近い行列もうまく処理することができる．
手続きは，簡約された系の次元にもとの系の次元を掛けた射影ベクトル 
を返す．
AbsoluteTiming[{sys, v} = mor[{l, s, d}, 50];]{loadR, stiffnessR, dampingR} = sys;
Dimensions /@ {loadR, stiffnessR, dampingR}initR = Transpose[v].init;
dinitR = LinearSolve[dampingR, -(stiffnessR.initR - loadR)];
Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
{ufunR} = u /. NDSolve[{
dampingR.u'[t] + stiffnessR.u[t] == loadR
, u[0] == initR
, u'[0] == dinitR
}, {u}, {t, 0, 1}
, Method -> {"EquationSimplification" -> "Residual"}
, EvaluationMonitor :> (currentTime = t;)
];
]
nrms = {};
errorGraphics = (
vals = ufun[#];
valsR = v.ufunR[#];
nrms = {nrms, Norm[vals - valsR, Infinity]};
)& /@ Range[0, 1, 0.1];
ListLinePlot[Flatten[nrms], PlotRange -> All]非定常係数を持つ非定常偏微分方程式
今度は，1つあるいは複数の係数が時間依存である非定常偏微分方程式について見てみよう．時間依存負荷 
を持つ次のモデル偏微分方程式を考えてみる．
これは，上の例と同じ定常境界条件を持つ．この例では，以下の負荷関数を使う．
変数と解のデータも上の例と同じである．有限要素法のメソッドデータも再利用できる．
偏微分方程式の初期化は，追加の成分である負荷係数を持つようになった．これは時間変数 
に依存する．
pdeCoefficients = InitializePDECoefficients[vd, sd
, "DiffusionCoefficients" -> {{-IdentityMatrix[3]}}
, "DampingCoefficients" -> {{1}}
, "LoadCoefficients" -> {{Sin[t]}}
]pdeCoefficients["All"]効率的な時間積分を行うために，係数は定常係数と非定常係数に分類される．定常係数は上の例と同じように，1度だけ評価される．時間依存成分を持つ係数は，時間ステップから時間ステップまで再評価され，定常離散化に追加されなければならない．
偏微分方程式の係数のInitializationは，すべての係数の分類を内部で行う．この分類は，チェックして利用することができる．それぞれの分類（例えば，定常係数）は，位置を含む2つのリストで表される．最初のリストにはスカラー係数が含まれ，2つ目のリストにはベクトルと行列の係数が含まれる．これをさらに詳しく見るために，一般的な偏微分方程式の形式を考えてみよう．
係数 
，
，
，
はスカラー，
，
，
はベクトル，
は 
×
行列である．返された2つのリストのうち最初のものは，スカラー係数（
，
，
，
），2つ目はベクトルと行列の係数（
，
，
，
）を表す．
pdeCoefficients["Transient"]Extract[pdeCoefficients["All"], #]& /@ pdeCoefficients["Transient"][[1]]Extract[pdeCoefficients["All"], #]& /@ pdeCoefficients["Transient"][[2]]"Stationary"係数は，指定されず，ゼロであるものも含めて，定常であるすべての係数を表す．
定常と非定常の成分に分類することによって，偏微分方程式の定常と非定常の離散化を別々に効率的に計算することが可能になる．
sdpde = DiscretizePDE[pdeCoefficients, methodData, sd, "Stationary"];
sdpde["SystemMatrices"]上の例では，文字列の引数"Stationary"はデフォルトであるが，必須ではない．離散化された定常偏微分方程式の負荷ベクトルは空である．負荷ベクトルには非定常の成分しか含まれていない．
sd2 = sd;
NDSolve`SetSolutionDataComponent[sd2, "Time", 3.];
DiscretizePDE[pdeCoefficients, methodData, sd2, "Transient"]["SystemMatrices"]離散化された非定常偏微分方程式の係数の負荷ベクトルが追加される．時点 
は，数値3で置き換えられている．その他の系の行列には非定常係数が含まれていないので，加えられていないことに注意する．
偏微分方程式の係数を分類するのは，定常係数はいったん離散化し，保存してから，時間積分のループで再利用することができるからである．つまり，時間変数 
の依存を示した係数だけを各ステップにおいて離散化すればいいということである．これによって，計算にかかる多くの労力を節約することができる．重要な幾何情報が事前に計算され，FEMMethodDataオブジェクトに保存されているので，それぞれの非定常成分の計算時間も最小に抑えることができる．
dbcsrhs[t_ ? NumericQ, u_] := Module[{l, s, d, m, tdpde},
{l, s, d, m} = sdpde["SystemMatrices"];
NDSolve`SetSolutionDataComponent[sd2, "Time", t];
tdpde = DiscretizePDE[pdeCoefficients, methodData, sd2, "Transient"];
{l, s, d, m} += tdpde["SystemMatrices"];
DeployBoundaryConditions[{l, s, d}, dbcs];
lsf[d_] := lsf[d] = LinearSolve[d];
lsf[d][l - s.u]
]右辺の計算には，すでに計算された偏微分方程式の定常係数が利用される．この例では，元帥行列は非定常成分を持たず，因数分解はキャッシュできる．減衰行列が時間依存の成分を持つ場合には，それぞれの時間ステップで，減衰行列がアップデートされ，LinearSolveステップがキャッシュされないので，上の右辺の計算は変更されなければならない．すると，因数分解が各時間ステップにおいて行われなければならない．
非定常境界条件がないので，初期条件は上と同じように設定される．
init = Table[{0.}, {methodData["DegreesOfFreedom"]}];
init[[dbcs["DirichletRows"]]] = dbcs["DirichletValues"];疎なパターンは，また減衰行列から導き出して設定される．この場合は，境界条件は減衰行列に配備されていない．
s = sdpde["SystemMatrices"][[2]];
d = sdpde["SystemMatrices"][[3]];{stiffnessSparsity = s["PatternArray"], dampingSparsity = d["PatternArray"]}Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
ufun = NDSolveValue[{
u'[ t] == rhs[t, u[t]]
, u[ 0] == init}, u, {t, 0, 1}
, Jacobian -> {"X" -> {Automatic, Sparse -> stiffnessSparsity}, "XP" -> {Automatic, Sparse -> dampingSparsity}}
, EvaluationMonitor :> (currentTime = t;)
]]allGraphics =
ElementMeshSurfacePlot3D[ufun[#], mesh, Boxed -> False, Background -> Gray, ColorFunctionScaling -> False]& /@ {0.01, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1.};
GraphicsGrid[Partition[allGraphics, 4], ImageSize -> Large]非線形の非定常係数を持つ非定常偏微分方程式
このセクションでは，1つあるいは複数の係数が非線形である非定常偏微分方程式について見てみる．非線形の減衰行列 
と非線形の拡散係数 
を持つ以下のモデル偏微分方程式を考えてみよう．
非線形の熱方程式が特定の例として使われる．温度 
，密度 
，温度依存の熱容量 
，熱依存の伝導性 
で以下の式を使う．
ここでは，物事を単純に保ち，有限要素法のプログラミング面に出来るだけ焦点を置くために，位相変化は一次元領域でモデル化され，単純性のためにパラメータが選ばれている．
位相変化の物理特性についての詳細は，位相変化がNDSolveレベルでどのようにモデル化できるかを説明する熱伝導の偏微分方程式モデリングのチュートリアルも参照されたい．
上の偏微分方程式を使って，位相変化をモデル化する．1つの位相から次の位相への変化をモデル化するのによく使われるのは，平滑な推移関数を使う方法である．
SmoothedStepFunction[fmin_, fmax_, ts_, m_] :=
Function[t, (fmin * Exp[ts * m] + fmax * Exp[m * t]) / (Exp[ts * m] + Exp[m * t])]
θ = SmoothedStepFunction[1, 0, 0, 50];Show[Plot[θ[T], {T, -1, 1}, ...], Graphics[{...}]]k[T_] := θ[T] * 2 + (1 - θ[T]) * 1 / 2;
Cp[T_] := θ[T] + (1 - θ[T]) * 2;
ρ = 2 * 10 ^ 6;Ω = Line[{{0}, {0.01}}];左の境界には一定流入量のエネルギーが存在し，右の境界には固定の温度が設定されている．時間の経過とともに，ものの初期位相は領域の左から右の端に変化する．
bcLeft = NeumannValue[5000, x == 0];
bcRight = DirichletCondition[T[t, x] == -10, x == 0.01];初期化段階は，上のセクションに示されるのと同じ方法で行われる．
nr = ToNumericalRegion[Ω];
vd = NDSolve`VariableData[{"DependentVariables" -> {T}, "Space" -> {x}, "Time" -> t}];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];initCoeffs = InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-k[T[t, x]] * IdentityMatrix[1]}}, "DampingCoefficients" -> {{ρ * Cp[T[t, x]]}}];
initBCs = InitializeBoundaryConditions[vd, sd, {{bcLeft, bcRight}}];
methodData = InitializePDEMethodData[vd, sd, Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.01 / 50}}];
NDSolveの方程式を設定するために，今度も上のセクションとほぼ同じ方法を取る．ここでは，上のセクションで 
と呼ばれたものを拡張し，同時に非線形係数を計算して非線形の減衰行列を含む．これにより，すべての偏微分方程式係数が1つの関数（ここでは微分代数方程式ソルバの剰余を計算するため，discretizePDEResidualと呼ぶことにする）に含まれるようになる．
sdpde = DiscretizePDE[initCoeffs, methodData, sd];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd];離散化関数は，時間 
，従属変数 
，従属変数 
の導関数を引数として得る．次に，時間従属偏微分方程式の係数と境界条件が計算され，定常成分に加えられる．その後，非線形成分が計算され，これも加えられる．そうすると，関数は離散化されたバージョンを返す．
discretizePDEResidual[t_ ? NumericQ, u_ ? VectorQ, dudt_ ? VectorQ] :=
Module[{l, s, d, tdpde, tbcs, nldpde, nlbcs},
NDSolve`SetSolutionDataComponent[sd, "Time", t];
NDSolve`SetSolutionDataComponent[sd, "DependentVariables", u];
l = sdpde["LoadVector"];
s = sdpde["StiffnessMatrix"];
d = sdpde["DampingMatrix"];
tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, "Transient"];
l += tdpde["LoadVector"];
s += tdpde["StiffnessMatrix"];
d += tdpde["DampingMatrix"];
nldpde = DiscretizePDE[initCoeffs, methodData, sd, "Nonlinear"];
nlbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, "Nonlinear"];
l += nldpde["LoadVector"];
s += nldpde["StiffnessMatrix"];
d += nldpde["DampingMatrix"];
DeployBoundaryConditions[{l, s, d}, nlbcs];
DeployBoundaryConditions[{l, s, d}, tbcs];
DeployBoundaryConditions[{l, s, d}, sbcs];
d.dudt + s.u - l
]
init = EvaluateOnElementMesh[{x}, -10, methodData["ElementMesh"]]["ValuesOnGrid"];疎なパターンを計算するために，ヘルパー関数を利用して，特定の時点における減衰行列と剛性行列を構築する．積分時間全体を通して有効である減衰行列と剛性行列の疎なパターンを作成する必要がある．したがって，任意の時点において構築された行列の疎なパターンを抽出し，その時点において，可能な時間従属係数が確実にアクティブであるようにする．しかし，ここでは，剛性行列と減衰行列のパターンだけが問題であるので，使われた厳密な非零の時間は重要ではない．
sparsityPatterns[t_, u_ ? VectorQ] := Module[
{s, d, tdpde, nldpde},
NDSolve`SetSolutionDataComponent[sd, "Time", t];
NDSolve`SetSolutionDataComponent[sd, "DependentVariables", u];
s = sdpde["StiffnessMatrix"];
d = sdpde["DampingMatrix"];
tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
s += tdpde["StiffnessMatrix"];
d += tdpde["DampingMatrix"];
nldpde = DiscretizePDE[initCoeffs, methodData, sd, "Nonlinear"];
s += nldpde["StiffnessMatrix"];
d += nldpde["DampingMatrix"];
{s["PatternArray"], d["PatternArray"]}
]{stiffnessSparsity, dampingSparsity} = sparsityPatterns[10 ^ -6, init]tEnd = 400;
Monitor[tTfun = NDSolveValue[{discretizePDEResidual[t, T[t], T'[ t]] == 0, T[0] == init}, T, {t, 0, tEnd}
, Method -> {"EquationSimplification" -> "Residual"}
, Jacobian -> {"X" -> {Automatic, Sparse -> stiffnessSparsity}, "XP" -> {Automatic, Sparse -> dampingSparsity}}
, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor]Tfun = ElementMeshInterpolation[{tTfun["Coordinates"][[1]], methodData["ElementMesh"]}, Partition[tTfun["ValuesOnGrid"], 1]]nframes = 30;
frames = Table[Labeled[Plot[Tfun[t, x], {x, 0, 0.01}, PlotRange -> {-20, 50}, PlotTheme -> "Detailed", PlotLabel -> Style["T(t,x)", 14], FrameLabel -> {"x", "Temperature"}, LabelStyle -> Directive[14], PlotLegends -> None, ImageSize -> Medium], Text[Style[Row[{"Time", " = ", N[t, 4]}], 14]]], {t, 0, tEnd, tEnd / nframes}];
ListAnimate[frames, SaveDefinitions -> True]熱方程式と適用可能な境界条件についての詳細は，「熱移動」のチュートリアルおよび「熱移動の偏微分方程式と境界条件」のガイドページに記載されている．
積分係数を持つ非定常偏微分方程式
非線形の非定常係数を持つ非定常偏微分方程式のセクションで提示されたメソッドは拡張して，非定常の積分偏微分方程式（省略形はPIDE）を解くことができる．各時間ステップでElementMeshInterpolationを構築し，NIntegrateを利用して現行の補間関数上で積分するというアイディアである．このアプローチは，積分だけに限らない．実際，数字を与えるどのような演算でも補間関数に使える．
Ω = Line[{{0}, {5}}];bcLeft = DirichletCondition[u[t, x] == 2 / Sqrt[Pi], x == 0];
bcRight = DirichletCondition[u[t, x] == 0, x == 5];初期化段階は，上のセクションで示される方法と同じように行われる．
nr = ToNumericalRegion[Ω];
vd = NDSolve`VariableData[{"DependentVariables" -> {u}, "Space" -> {x}, "Time" -> t}];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];initCoeffs = InitializePDECoefficients[vd, sd,
"DampingCoefficients" -> {{1}}, "DiffusionCoefficients" -> {{-IdentityMatrix[1]}},
"ReactionCoefficients" -> {{-x}}
];
initBCs = InitializeBoundaryConditions[vd, sd, {{bcLeft, bcRight}}];
methodData = InitializePDEMethodData[vd, sd];sdpde = DiscretizePDE[initCoeffs, methodData, sd];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd];discretizePDEResidualと呼ばれる離散化関数で，線形成分をロードし，非定常で非線形の成分を加えてから，現行の補間関数を構築する．次に，構築された関数上で数値積分を行い，積分を含む，見付からない偏微分方程式係数を構築する．見付からない係数を離散化し，その結果を系の行列に加える．そうすると，関数は離散化されたバージョンを返す．
各時間ステップでElementMeshInterpolationを作成するために，メッシュをNumericalRegionから抽出する．
mesh = methodData["ElementMesh"]discretizePDEResidual[t_ ? NumericQ, u_ ? VectorQ, dudt_ ? VectorQ] :=
Module[{l, s, d, tdpde, tbcs, nldpde, nlbcs,
currentU, integral, integralCoeffs, integraldpde},
NDSolve`SetSolutionDataComponent[sd, "Time", t];
NDSolve`SetSolutionDataComponent[sd, "DependentVariables", u];
l = sdpde["LoadVector"];
s = sdpde["StiffnessMatrix"];
d = sdpde["DampingMatrix"];
tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, "Transient"];
l += tdpde["LoadVector"];
s += tdpde["StiffnessMatrix"];
d += tdpde["DampingMatrix"];
nldpde = DiscretizePDE[initCoeffs, methodData, sd, "Nonlinear"];
nlbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, "Nonlinear"];
l += nldpde["LoadVector"];
s += nldpde["StiffnessMatrix"];
d += nldpde["DampingMatrix"];
(* discretize the integral *)
currentU = ElementMeshInterpolation[{mesh}, u];
integral = NIntegrate[x * currentU[x], {x}∈mesh];
integralCoeffs = InitializePDECoefficients[vd, sd, "ReactionCoefficients" -> {{integral}}];
integraldpde = DiscretizePDE[integralCoeffs, methodData, sd];
l += integraldpde["LoadVector"];
s += integraldpde["StiffnessMatrix"];
integraldpde = DiscretizePDE[integralCoeffs, methodData, sd, "Transient"];
l += integraldpde["LoadVector"];
s += integraldpde["StiffnessMatrix"];
integraldpde = DiscretizePDE[integralCoeffs, methodData, sd, "Nonlinear"];
l += integraldpde["LoadVector"];
s += integraldpde["StiffnessMatrix"];
DeployBoundaryConditions[{l, s, d}, nlbcs];
DeployBoundaryConditions[{l, s, d}, tbcs];
DeployBoundaryConditions[{l, s, d}, sbcs];
d.dudt + s.u - l
]厳密には，この場合，積分係数は線形であるので，積分係数の非定常で非線形な成分の離散化は不要である．一般化のためにコードに非定常で非線形な積分係数の計算も含む．
init = EvaluateOnElementMesh[{x}, 2 / (Sqrt[Pi] * Exp[x ^ 2]), mesh]["ValuesOnGrid"];sparsityPatterns[t_, u_ ? VectorQ] := Module[
{s, d, tdpde, nldpde},
NDSolve`SetSolutionDataComponent[sd, "Time", t];
NDSolve`SetSolutionDataComponent[sd, "DependentVariables", u];
s = sdpde["StiffnessMatrix"];
d = sdpde["DampingMatrix"];
tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
s += tdpde["StiffnessMatrix"];
d += tdpde["DampingMatrix"];
nldpde = DiscretizePDE[initCoeffs, methodData, sd, "Nonlinear"];
s += nldpde["StiffnessMatrix"];
d += nldpde["DampingMatrix"];
{s["PatternArray"], d["PatternArray"]}
]{stiffnessSparsity, dampingSparsity} = sparsityPatterns[10 ^ -6, init]tEnd = 5;
Monitor[tufun = NDSolveValue[{discretizePDEResidual[t, u[t], u'[ t]] == 0, u[0] == init}, u, {t, 0, tEnd}
, Method -> {"EquationSimplification" -> "Residual"}
, Jacobian -> {"X" -> {Automatic, Sparse -> stiffnessSparsity}, "XP" -> {Automatic, Sparse -> dampingSparsity}}
, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor]ufun = ElementMeshInterpolation[{tufun["Coordinates"][[1]], mesh}, Partition[tufun["ValuesOnGrid"], 1]]Plot3D[ufun[t, x], {x, 0, 5}, {t, 0, tEnd}]結合偏微分方程式
このセクションでは，結合偏微分方程式がどのように解かれるかについて説明する．基本的には，ワークフローは単一の偏微分方程式と同じである．違いを下に述べる．
負荷による梁の変形
例として，構造力学の問題として，梁の変形について考える．変形は，2つの従属変数における結合偏微分方程式としてモデル化できる．
Clear[u, v, x, y]
op = {Inactive[Div][{{0, -((Y*ν)/(1 - ν^2))}, {-((Y*(1 - ν))/(2*(1 - ν^2))), 0}} .
Inactive[Grad][v[x, y], {x, y}], {x, y}] + Inactive[Div][{{-(Y/(1 - ν^2)), 0}, {0, -((Y*(1 - ν))/(2*(1 - ν^2)))}} .
Inactive[Grad][u[x, y], {x, y}], {x, y}], Inactive[Div][{{0, -((Y*(1 - ν))/(2*(1 - ν^2)))}, {-((Y*ν)/(1 - ν^2)), 0}} .
Inactive[Grad][u[x, y], {x, y}], {x, y}] + Inactive[Div][{{-((Y*(1 - ν))/(2*(1 - ν^2))), 0}, {0, -(Y/(1 - ν^2))}} .
Inactive[Grad][v[x, y], {x, y}], {x, y}]};平面応力の演算子では，YはYoungの係数であり，ν はポアソン(Poisson)の率である．これらは物質の特性である．
あるいは，関数SolidMechanicsPDEComponentを使って，偏微分方程式の成分を生成することもできる．SolidMechanicsPDEComponentの関数ページにもテキスト形式での方程式が含まれている．
Ω = ImplicitRegion[True, {x, y}];
mesh = ToElementMesh[Ω, {{0, 5}, {0, 1}}, "MaxCellMeasure" -> 0.1];
mesh["Wireframe"]境界条件は以下の通りである．左辺では，梁は x と y の両方向に固定されている．右辺の境界では，1単位の下向きの圧力が適用される．NeumannValueは，負の y 方向に梁にかかる圧力を指定することによって境界負荷を強制する．
Subscript[Γ, Du] = DirichletCondition[u[x, y] == 0, x == 0];
Subscript[Γ, Dv] = DirichletCondition[v[x, y] == 0, x == 0];
Subscript[Γ, Nv] = NeumannValue[-1, x == 5];初期化は，上のセクションでの場合と同じ成分を持つ．まず，変数と解のデータが設定される．
vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u, v}, {x, y}}];結合偏微分方程式については，2つ以上の従属変数が与えられなければならない．
sd = NDSolve`SolutionData["Space" -> ToNumericalRegion[mesh]];メソッドデータの初期化は，単一の従属変数の場合と同じようになされる．
methodData = InitializePDEMethodData[vd, sd]FEMMethodDataの表示形式は，自由度，それぞれの従属変数の補間次数，そして積分次数を示す．
方程式系の自由度は，要素メッシュのノード，従属変数の数，および従属変数の補間次数を組み合せた結果である．
Length[mesh["Coordinates"]] * Length[NDSolve`SolutionDataComponent[vd, "DependentVariables"]]
methodData["DegreesOfFreedom"]結合偏微分方程式の場合は，係数がもう少し関係してくる．結合偏微分方程式は，基本的に，単一の偏微分方程式が集まったものである．しかし，それぞれの成分がそれぞれの従属変数に対して存在する．次の2つの従属変数 
と 
の偏微分方程式を考えてみよう．
結合偏微分方程式は
で定義される．係数 
と 
は，それぞれの偏微分方程式でスカラーである．各偏微分方程式において，
， 
，
は長さ 
のベクトルであり，
は 
×
行列である．InitializePDECoefficientsの入力係数は，上の結合偏微分方程式とモデル化される結合偏微分方程式との対応によって求められる．
diffusionCoefficients = "DiffusionCoefficients" -> {{{{-(Y/1 - ν^2), 0}, {0, -(Y (1 - ν)/2 (1 - ν^2))}}, {{0, -(Y ν/1 - ν^2)}, {-(Y (1 - ν)/2 (1 - ν^2)), 0}}}, {{{0, -(Y (1 - ν)/2 (1 - ν^2))}, {-(Y ν/1 - ν^2), 0}}, {{-(Y (1 - ν)/2 (1 - ν^2)), 0}, {0, -(Y/1 - ν^2)}}}} /. {Y -> 10 ^ 3, ν -> 33 / 100};結合偏微分方程式の場合には，各偏微分方程式の各従属変数について拡散係数行列がある．
initCoeffs = InitializePDECoefficients[vd, sd, {diffusionCoefficients}]PDECoefficientDataの表示形式は，初期化された係数について，系の大きさ，つまり従属変数の数（2個）と演算子の大きさ，つまり空間次元（2）を示す．
次に，境界条件を設定する．従属変数のそれぞれには，境界条件のリストが含まれる．これが境界条件を偏微分方程式の系の特定部分に関連付ける．境界条件の最初の集合は，結合系の最初の偏微分方程式に，境界条件の2つ目の集合は結合系の2つ目の偏微分方程式に，という風に関連付けられる．これによって，境界条件を一意的に偏微分方程式の部分に関連付けることが可能である．
initBCs = InitializeBoundaryConditions[vd, sd, {{Subscript[Γ, Du]}, {Subscript[Γ, Dv], Subscript[Γ, Nv]}}]discretePDE = DiscretizePDE[initCoeffs, methodData, sd]DiscretizedPDEDataの表示形式は，一次元の系の行列を示す．
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"]MatrixPlot[stiffness]剛性行列の可視化は，もとになる偏微分方程式の係数が2つの従属変数の系であることを示す．自由度は，2つの従属変数間で分けられる．
split = Span@@@Transpose[{Most[# + 1], Rest[#]}&[methodData["IncidentOffsets"]]]discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd]DeployBoundaryConditions[{load, stiffness}, discreteBCs]方程式系はこれでLinearSolveを使って解くことができるようになる．
Short[solution = LinearSolve[stiffness, load]]それぞれの従属変数の補間関数を作成することができる．解の形式LinearSolveは，
と 
への両方の解を含む1つの解ベクトルであるので，解ベクトルは分けられなければならない．
ufun = ElementMeshInterpolation[{mesh}, solution[[split[[1]]]]]
vfun = ElementMeshInterpolation[{mesh}, solution[[split[[2]]]]]ContourPlot[ufun[x, y], {x, y}∈mesh, ColorFunction -> "Temperature", AspectRatio -> Automatic]
ContourPlot[vfun[x, y], {x, y}∈mesh, ColorFunction -> "Temperature", AspectRatio -> Automatic]dmesh = ElementMeshDeformation[mesh, {ufun, vfun}]Show[{
mesh["Wireframe"],
dmesh["Wireframe"["ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]固体力学と適用可能な境界条件についての詳細は，「固体力学」のモノグラフおよび「固体力学の偏微分方程式と境界条件」のガイドページに記載されている．
揺れる梁—非定常の結合偏微分方程式
非定常の結合偏微分方程式には，これ以上はほとんど何も必要ではない．以下の例では，梁を偏向させ，レイリー(Rayleigh)の減衰で動きを弱める方法を示す．レイリーの減衰では，減衰行列
は，質量行列と剛性行列の線形の組合せである．2つの従属変数を持つ非定常の結合偏微分方程式についての系の行列の構造は，次のように表される．
massCoefficients = "MassCoefficients" -> {{1, 0}, {0, 1}}initCoeffs = InitializePDECoefficients[vd, sd, {diffusionCoefficients, massCoefficients}]discretePDE = DiscretizePDE[initCoeffs, methodData, sd]DiscretizedPDEDataの表示形式は，一次元の系の行列を示す．
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"]減衰行列の疎な配列は空である．これは，質量行列の係数が設定され，ここで使用されるデフォルトの減衰行列の係数がゼロであるからである．
MatrixPlot[mass]質量行列は，ブロック対角が追加され，対角外のブロックがゼロである，偏微分方程式系数の構造を示す．
レイリーの減衰は，質量行列と剛性行列の線形の組合せとして計算される．
rayleighDamping = 0.1 * mass + 0.04 * stiffness;DeployBoundaryConditions[{load, stiffness, rayleighDamping, mass}, discreteBCs]初期条件と初期条件の導関数は，ゼロの設定される．これは，ディリクレ境界条件がすべてゼロであり，梁が静から始まるので妥当である．
dof = methodData["DegreesOfFreedom"];
init = dinit = ConstantArray[0, {dof, 1}];dampingSparsity = ArrayFlatten[{{mass["PatternArray"], mass["PatternArray"]}, {rayleighDamping["PatternArray"], rayleighDamping["PatternArray"]}}]MatrixPlot[dampingSparsity]stiffnessSparsity = ArrayFlatten[{{stiffness["PatternArray"], stiffness["PatternArray"]}, {stiffness["PatternArray"], stiffness["PatternArray"]}}]MatrixPlot[stiffnessSparsity]次のステップは，NDSolveを使って行う実際の時間積分である．シミュレーションのために，単一の記号的変数 
を導入する．記号的変数 
は，
と 
のベクトルを組み合せたものを表し，NDSolveには記号 
が初期条件と系の行列の大きさの両方から，
と 
の長さを組み合せたベクトルを表さなければならないことが分かっている．言い換えれば，それぞれの自由度について変数を導入する必要はない．
Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
tfun = NDSolveValue[{
mass.uv''[ t] + rayleighDamping.uv'[ t] + stiffness.uv[ t] == load
, uv[ 0] == init, uv'[ 0] == dinit}, uv, {t, 0, 45}
, Method -> {"EquationSimplification" -> "Residual"}
, Jacobian -> {"X" -> {Automatic, Sparse -> stiffnessSparsity}, "XP" -> {Automatic, Sparse -> dampingSparsity}}
, EvaluationMonitor :> (currentTime = t;)
]]それぞれのタイムスタンプで同じ量の変形を表示するには，"ScalingFactor"を1に設定する．
graphics = Function[t,
dmesh = ElementMeshDeformation[mesh, Part[tfun[t], #]& /@ split, "ScalingFactor" -> 1];
Show[{
mesh["Wireframe"["MeshElement" -> "BoundaryElements"]],
dmesh["Wireframe"["ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]
}, PlotRange -> {{0, 5.5}, {1.1, -1.}}]] /@ Range[0, 45, 1];ListAnimate[graphics, SaveDefinitions -> True]固体力学と適用可能な境界条件についての詳細は，「固体力学」のモノグラフおよび「固体力学の偏微分方程式と境界条件」のガイドページに記載されている．
揺れる動的に負荷された梁
このセクションでは，時間依存のボディ負荷 
と時間依存の境界条件を加えることによって，上の例をさらに1ステップ発展させる．
そのためには，変数と解のデータは時間変数集合を持つ必要がある．
vd = NDSolve`VariableData[{"Time", "DependentVariables", "Space"} -> {t, {u, v}, {x, y}}];sd = NDSolve`SolutionData[{"Time", "Space"} -> {0., ToNumericalRegion[mesh]}];methodData = InitializePDEMethodData[vd, sd]
loadCoefficients = "LoadCoefficients" -> {{4 Sin[t]}, {0}}次に，NeumannValueを通して圧力を指定することによって，y 方向に適用される時間依存の境界の圧力を加える．
Subscript[Γ, Nv] = NeumannValue[-1 - 1 / 2 * Sin[t], x == 5];initCoeffs = InitializePDECoefficients[vd, sd, {diffusionCoefficients, massCoefficients, loadCoefficients}]initBCs = InitializeBoundaryConditions[vd, sd, {{Subscript[Γ, Du]}, {Subscript[Γ, Dv], Subscript[Γ, Nv]}}]離散化された定常偏微分方程式の係数と境界条件を事前に計算する．
sdpde = DiscretizePDE[initCoeffs, methodData, sd, "Stationary"]sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, "Stationary"]
現行時間 
，
，
，
の値を得るためのヘルパ関数を書く．これらは1つのベクトル 
および 
に結合される．現行時間は解のデータで更新される必要がある．定常の系の行列が離散化された偏微分方程式から抽出される．これらはそれぞれの時間ステップで変化しないので，再度計算されることはないということを覚えておくことが重要である．次に，偏微分方程式の非定常部分と境界条件が離散化される．非定常と定常の系の行列が足し合わされ，レイリー減衰が計算される．定常と非定常の境界条件はどちらも配備される．関数は離散化された 
を返す．
rhs[t_ ? NumericQ, uv_, duv_] := Module[{l, s, d, m, tdpde, tbcs, rayleighDamping},
NDSolve`SetSolutionDataComponent[sd, "Time", t];
{l, s, d, m} = sdpde["SystemMatrices"];
tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, "Transient"];
{l, s, d, m} += tdpde["SystemMatrices"];
rayleighDamping = 0.1 * m + 0.04 * s;
DeployBoundaryConditions[{l, s, rayleighDamping, m}, tbcs];
DeployBoundaryConditions[{l, s, rayleighDamping, m}, sbcs];
l - s.uv - rayleighDamping.duv
]初期条件と初期条件の導関数は零に設定される．ディリクレ条件はすべて零であり張りは残りから始まるため，この設定は妥当である．
dof = methodData["DegreesOfFreedom"];
init = dinit = ConstantArray[0, {dof, 1}];疎なパターンの指定には，偏微分方程式のの定常離散化を使う．定常コンポーネントは負荷ベクトルにのみ存在するので，これは妥当である．言い換えれば，質量行列と剛性行列に時間依存のコンポーネントがない場合にのみ以下は正しい．
mass = sdpde["MassMatrix"];
stiffness = sdpde["StiffnessMatrix"];
rayleighDamping = 0.1 * mass + 0.04 * stiffness;dampingSparsity = ArrayFlatten[{{mass["PatternArray"], mass["PatternArray"]}, {rayleighDamping["PatternArray"], rayleighDamping["PatternArray"]}}]stiffnessSparsity = ArrayFlatten[{{stiffness["PatternArray"], stiffness["PatternArray"]}, {stiffness["PatternArray"], stiffness["PatternArray"]}}]Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
tfun = NDSolveValue[{
mass.uv''[ t] == rhs[t, uv[t], uv'[t]]
, uv[ 0] == init, uv'[ 0] == dinit}, uv, {t, 0, 45}
, Method -> {"EquationSimplification" -> "Residual"}
, Jacobian -> {"X" -> {Automatic, Sparse -> stiffnessSparsity}, "XP" -> {Automatic, Sparse -> dampingSparsity}}
, EvaluationMonitor :> (currentTime = t;)
]]graphics = Function[t,
dmesh = ElementMeshDeformation[mesh, Part[tfun[t], #]& /@ split];
Show[{
mesh["Wireframe"["MeshElement" -> "BoundaryElements"]],
dmesh["Wireframe"["ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]
}, PlotRange -> {{0, 5.5}, {1.5, -1.5}}]] /@ Range[0, 45, 1];ListAnimate[graphics, SaveDefinitions -> True]ボディの負荷による 
方向の動きに，および 
における時間従属の境界負荷により引き起こされる 
方向の動きに収縮と膨張がある．
固体力学と適用可能な境界条件についての詳細は，Solid Mechanics（固体力学）のモノグラフおよびSolid Mechanics PDEModelsのガイドページに記載されている．
大規模な有限要素法分析
NDSolveで実装される有限要素法は，スピードの点で最適化されている．この最適化により失われるものもある．数値計算法では，スピードとメモリの消費が両立することがないことはよく知られている．つまり，実装が速ければ速いほど，手元の問題を解くのにより多くのメモリを必要とする．しかしながら，実装された有限要素法をさまざまな方法を通して調整し，偏微分方程式の解を求めるのに追加の時間を費やすことによって，より少ない量のメモリを使うようにすることは可能である．「有限要素法で偏微分方程式を解く」にはいくつかの方法が提案されている．以下のセクションでは，有限要素解法で必要なメモリ量を減らすためのさまざまなその他のオプションについて述べる．これらの提案には，一般的なものから特定のオプションまでいろいろある．
一般にWolfram言語は，先に計算された結果への参照を保存する．この仕組みにより，将来の計算には不要であるデータが累積されることになってしまう．履歴の長さを制御する仕組みを設定し，前の出力を保存しないようにすることが可能である．
$HistoryLength = 0;- 幾何学データの事前計算(InitializePDEMethodData)．
- 偏微分方程式系を解く(DiscretizePDE)．
- 方程式系を解く(LinearSolve)．
有限要素分析には，メモリの妨げになるものが2つある． 剛性行列の作成(DiscretizePDE)と方程式系の解法(LinearSolve)である．メッシュ生成と幾何学データの事前計算の両方に対するオプションが，離散化と解法の最中にメモリの必要量に影響を与える．さらに，デフォルトで，離散化段階そのものが2つのステップで行われる．まず，特定の系の行列についてすべての要素が計算され，その後，それらの要素から行列が組み立てられる．メモリの消費を少なくするために取られたアクションが，計算のために余分な時間を必要とする，あるいはその特定のアクションが解の品質にも影響を与えるという場合には，追加の区別が必要である．
オプション（NDSolveに直接渡すこともできる）は，例を使って説明するのが一番である．大量のメモリを必要とする典型的な例に，3Dの適用例および／もしくは数多くの結合偏微分方程式を必要とする例がある．
op = SolidMechanicsPDEComponent[{{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}}, <|"PoissonRatio" -> 1 / 3, "YoungModulus" -> 5|>];Subscript[Γ, DX] = DirichletCondition[u[x, y, z] == 0, z == 0];
Subscript[Γ, DY] = DirichletCondition[v[x, y, z] == 0, z == 0];
Subscript[Γ, DZ] = DirichletCondition[w[x, y, z] == 0, z == 0];
Subscript[Γ, NZ] = NeumannValue[1, z == 1];重要な最初のステップは，要素メッシュの生成である．さまざまな点が，メッシュ状で偏微分方程式を解くのに必要なメモリに影響する．"MeshOrder"は，どのくらいのメモリが必要かということに対して，大きく影響する．例えば，一次四面体は4つのノードを持つのに対し，二次四面体は，10個のノードを持つ．"MeshOrder"を減らすと，「有限要素法で偏微分方程式を解く」でも述べたように，期待される解の品質も低下する．
vd = NDSolve`VariableData[{"DependentVariables" -> {u, v, w}, "Space" -> {x, y, z}}]bmesh = ToBoundaryMesh[FullRegion[3], {{0, 1}, {0, 1}, {0, 1}}, "MaxBoundaryCellMeasure" -> Infinity]次の離散化段階における行列の組立ての際に，メッシュからのすべての要素がまず計算されてから，次に組み立てられる．これはつまり，ある時点で，要素と系の行列がメモリにあるということである．これは，オプション"MeshElementBlocks"を正の整数に設定することによって，避けることができる．その後，行列の要素はたくさんのブロックに分割される．要素メッシュのブロックが計算されると，これが系の行列に組み立てられる．その後ブロックは破棄され，次のブロックが系の行列に組み立てられ，すべてのブロックが計算，組み立てられるまで続けられる．このスキームでは，最大のメモリ消費は，完全に組み立てられた系の行列であり，最大の要素メッシュブロックの大きさ大きさとその計算された要素である．"MeshElementBlocks"を大きく設定しすぎると，組立て過程に時間がかかる．
AbsoluteTiming[
MaxMemoryUsed[mesh = ToElementMesh[bmesh, "MaxCellMeasure" -> 0.00005, "MeshElementBlocks" -> 10]] / 1024. ^ 2]meshByteCount[mesh] / 1024. ^ 2sd = NDSolve`SolutionData[{"Space" -> ToNumericalRegion[mesh]}]方程式の処理中に，InitializePDEMethodDataへの呼出しがなされる．この呼出しの間，有限要素法に必要な幾何学データが事前計算される．NDSolve`ProcessEquationsが方程式の処理のみに使われる場合には，幾何学的な事前処理はオフにすることができる．
{nds} = NDSolve`ProcessEquations[{op == {0, 0, Subscript[Γ, NZ]}, Subscript[Γ, DX], Subscript[Γ, DY], Subscript[Γ, DZ]}, {u, v, w}, {x, y, z}∈mesh, Method -> {"FiniteElement", "PrecomputeGeometryData" -> False}];dfc = -nds["FiniteElementData"]["PDECoefficientData"]["All"][[2, 1]];
ClearAll[nds];initCoeffs = InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> dfc]initBCs = InitializeBoundaryConditions[vd, sd, {{Subscript[Γ, DX]}, {Subscript[Γ, DY]}, {Subscript[Γ, DZ], Subscript[Γ, NZ]}}]離散化段階において，事前計算された幾何学データが，系の行列への各要素の貢献を計算するのにひっようである．非定常系においては，幾何学的な事前計算をオフにすると，このデータが各時間ステップごとに再計算されなければならない．これによって，系を解くのにより長い時間がかかるが，非常に大規模な系の場合には，解を得るためにこの方法を取らなければならないかもしれない．
AbsoluteTiming[methodData = InitializePDEMethodData[vd, sd]]MemoryInUse[] / 1024. ^ 2AbsoluteTiming[MaxMemoryUsed[discretePDE = DiscretizePDE[initCoeffs, methodData, sd]] / 1024. ^ 2]MemoryInUse[] / 1024. ^ 2load = discretePDE["LoadVector"];stiffness = discretePDE["StiffnessMatrix"]ByteCount[stiffness] / 1024. ^ 2discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd]次に，非常に細かい点について考慮する．系の行列はDiscretizedPDEDataのデータオブジェクトに保存される．これらは上で抽出され，load と stiffness の記号に保存される．境界条件の配備は，これらの系の行列を変更する．しかし，これらはDiscretizedPDEData（もう一度未変更の系の行列を抽出する可能性がある）に保存されるので，複製バージョンが作成される．これを避け，メモリ消費量を節約するために，系の行列について1度だけの参照が存在するように，DiscretizedPDEDataオブジェクトを削除することができる．
ClearAll[discretePDE];
ClearAll[methodData];MemoryInUse[] / 1024. ^ 2DeployBoundaryConditions[{load, stiffness}, discreteBCs]MemoryInUse[] / 1024. ^ 2AbsoluteTiming[MaxMemoryUsed[solution = LinearSolve[stiffness, load, Method -> "Pardiso"];] / 1024. ^ 2]