微分代数方程式の数値解法

はじめに	微分代数方程式の指数減少
微分代数方程式の指数	微分代数方程式の矛盾のない初期化
微分方程式の取扱い	微分代数方程式の例題
微分代数方程式の解法

はじめに

一般に，常微分方程式系は以下の標準形

で表すことができる．従属変数の導関数は，過渡的な独立変数と従属変数について明示的に表される．関数に十分な連続性がある限り，従属変数の値が独立変数の特定の値において与えられるという初期値問題に対する一意解は常に見付けることができる．

微分代数方程式では，一般に導関数は明示的に表されない．実際，従属変数の中にはその導関数が通常方程式に現れないものもある．例えば，以下の方程式

にはの導関数が明示的には含まれていない．このような変数は代数的変数と呼ばれることがよくある．

微分代数方程式系の一般形式は以下の通りである．

微分代数方程式系を解くには，多数のステップが必要となることが多い．以下のフローチャートはNDSolveで微分代数方程式を解くことに関する一般的なプロセスを示すものである．

NDSolveで微分代数方程式系を解くのに関わるプロセスを示すフローチャート

通常，微分代数方程式系は独立変数について微分することで常微分方程式系に変換することができる．微分代数方程式の指数とは，常微分方程式系を得るために微分代数方程式系を微分する回数である．

NDSolveに組み込まれている微分代数方程式ソルバは，指数1の系に使うものなので，高指数の系の場合は解を得るために指数減少法が必要になる場合がある．NDSolveは指数減少法を実行することができる．系に1より大きい指数が含まれていると，NDSolveは微分代数方程式を解くために必要なメッセージとステップ数を出力する．

高指数の微分代数方程式での第1ステップは，系の指数を減少させることである．指数を減少させるために微分するプロセスは，指数減少法と呼ばれ，NDSolveの記号的プロセスで実行できる．

指数減少法のプロセスによって，新しい同等の系になる．1つのアプローチでは，この新しい系は微分された変数の代りに新しい（ダミー）変数を代入することにより構築される．これだと，一意に解くことのできる拡張された系になる．別のアプローチでは，すべての変数に対する微分が構築されるまで，もとの系を回微分する．先行する回微分された系は，不変量として扱われる．

この新しい指数減少系を解くためには，矛盾のない初期条件集合を見付けなければならない．標準形の常微分方程式系は常に，開始時にの値を与えることにより初期化することができる．しかし．微分代数方程式の場合は残差方程式を満足する初期条件を見付けることは難しいことがある．これは結局，変数のいくつかだけが独立して指定されている非線形の代数系を解くことになるのである．このことに関しては後のセクションで詳しく説明する．さらに，初期化は一貫していなければならない．つまり，残差方程式の微分も満足されなければならないのである．通常，高指数の系になるほど初期化が難しい．この場合，一般にNDSolveではの要素がどのように相互作用を及ぼし合うか見ることができないので，指数減少が自動でできない．しかし，NDSolveには，オプションを介して利用できる，指数減少を実行し，微分代数方程式の矛盾のない初期化を見付けるためのさまざまなメソッドがある．

このドキュメントの残りの部分は次のようになっている．まず微分代数方程式の指数に関連した用語を紹介し，次に指数1の微分代数方程式の解き方を説明する．その次に高指数の微分代数方程式の指数減少法を説明する．微分代数方程式の矛盾のない初期化も扱う．ここまでの例題に加え，最後のセクションではさらに例題を提示する．

この例で使用するパッケージをロードする：

Wolfram Language code:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

微分代数方程式の指数

一般的な微分代数方程式系では，系の指数とはおよびについて一意にを解くために必要な，系に含まれる方程式の一部あるいはすべてを微分する最小回数である．文献では指数の考え方について多少異なるものもあるが，このドキュメントでは，上記の考えを使う．

次の微分代数方程式を考える．

の微分方程式を得るためには，最初の方程式を1回微分しなければならない．これで，新しい系ができる．

微分された方程式系には，説明の必要なが含まれている．の方程式を得るために，2つ目の方程式は3回微分しなければならない．これにより以下が得られる．

系の一部はの導関数方程式を得るために3回微分しなければならないので，この微分代数方程式の指数は3である．微分された最終的な系は，指数0の系と言われる．微分によって得られた常微分方程式系を見るためには，適切な代入を行わなければならない．この場合は最初の方程式から2つ目の方程式を引くことである．結果の方程式は次のように書くことができる．

同様に，指数1の系は2つ目の方程式を2回微分することにより見付けることができ，以下の系となる．

系の指数を1に減少させるというのは一般的なことである．内在する積分ルーチンがそれを効率的に扱えるからである．

微分代数方程式の指数は明らかでない場合がある．これは初期条件が何であるかにより，指数1あるいは指数2の動作を見せる可能性のある系を考える．

次は3つの方程式を含む微分代数方程式系であるが，微分項は1つしかない．

Wolfram Language code:

DAE = (|                                     |
| ----------------------------------- |
| Derivative[1][x1][t] == x3[t]        |
| x2[t] (1 - x2[t]) == 0               |
| x1[t] x2[t] + x3[t] (1 - x2[t]) == t |);

この微分代数方程式では，2つ目の方程式でが0か1であることが要求されているため初期条件は無ではない．

ならば，残りの2つの方程式は指数1の系になる．3つ目の方程式の微分が2つの常微分方程式系を与えるからである．

初期条件である一方，の初期条件は変更することができる．

以下は

に指定された初期条件を持ち，

で始めた微分代数方程式系を解く：

Wolfram Language code:

sol = NDSolve[{DAE, Subscript[x, 1][0] == -1, Subscript[x, 2][0] == 0}, {Subscript[x, 1], Subscript[x, 2], Subscript[x, 3]}, {t, 0, 1}]

Wolfram Language code: Plot[Evaluate[{Subscript[x, 1][t], Subscript[x, 2][t], Subscript[x, 3][t]} /. First[sol]], {t, 0, 1}]

一方，ならば，残りの2つの方程式は指数2の系になる．3つ目の方程式の微分がを与え，それを最初の方程式に代入するとになる．これを微分して常微分方程式を得る．

およびであるため，初期条件には追加の自由度はない．

以下は

から始めて微分代数方程式系を解く：

Wolfram Language code: sol = NDSolve[{DAE, Subscript[x, 2][0] == 1}, {Subscript[x, 1], Subscript[x, 2], Subscript[x, 3]}, {t, 0, 1}]

Wolfram Language code: Plot[Evaluate[{Subscript[x, 1][t], Subscript[x, 2][t], Subscript[x, 3][t]} /. First[sol]], {t, 0, 1}]

この例は，必要な初期条件の数と微分代数方程式の指数の両方が，実際に見付けられる解に影響を及ぼすことを示している．

微分方程式の取扱い

微分方程式を解くために，NDSolveはユーザが指定した系を3つの形式のいずれかに変換する．系がどのように構築されるかにより．異なる積分法が選ばれるので，このステップは非常に重要である．

オプションMethod->{"EquationSimplification"->simplification}を使うと，どの形式で方程式を表すかを指定することができる．次の簡約化メソッドは，"EquationSimplification"オプションに対して指定することができる．

Automatic	系を形式で記号的に解こうとする．解が見付からない，あるいは解くのに時間がかかりすぎる場合は，"MassMatrix"あるいは"Residual"メソッドを使って簡約化を試みる．これがデフォルトである
"Solve"	可能であれば，系を形式で記号的に解こうとする
"MassMatrix"	可能であれば，系を形式に簡約する
"Residual"	方程式の左辺から右辺を引いて，残差関数を作る

方程式の簡約化のオプション

次の系を考える

デフォルトでは独立変数と従属変数に関して明示的に導関数について解くことが試みられる．効率を考えて，NDSolveはまずLinearSolveを使って明示的な系を見つけようとする．しかし，それに失敗すると，Solveを使って記号的に系が解かれる．

NDSolveは自動的に系を明示的な形式にしようとする：

Wolfram Language code:

state = First@NDSolve`ProcessEquations[{{x'[t] + y'[t] == Sin[t] - z'[t], y'[t] + z'[t] == Sin[t] + x[t], x'[t] + z'[t] == y[t] + Cos[t]}, {x[0] == 1, y[0] == 1, z[0] == 1}}, {x, y, z}, {t}];
state["NumericalFunction"]["FunctionExpression"]

"Solve"メソッドを使うことは記号解が見付けるときのデフォルトと同じである．

オプション"MassMatrix"を使って，形式

の系を生成する：

Wolfram Language code:

state = First@NDSolve`ProcessEquations[{{x'[t] + y'[t] == Sin[t] - z'[t], y'[t] + z'[t] == Sin[t] + x[t], x'[t] + z'[t] == y[t] + Cos[t]}, {x[0] == 1, y[0] == 1, z[0] == 1}}, {x, y, z}, {t}, Method -> {"EquationSimplification" -> "MassMatrix"}];
state["NumericalFunction"]["FunctionExpression"]

Wolfram Language code: state["MassMatrix"] // MatrixForm

オプション"Residual"を使って，形式

の系を生成する：

Wolfram Language code:

state = First@NDSolve`ProcessEquations[{{x'[t] + y'[t] == Sin[t] - z'[t], y'[t] + z'[t] == Sin[t] + x[t], x'[t] + z'[t] == y[t] + Cos[t]}, {x[0] == 1, y[0] == 1, z[0] == 1}}, {x, y, z}, {t}, Method -> {"EquationSimplification" -> "Residual"}];
state["NumericalFunction"]["FunctionExpression"]

系が残差形式になると，導関数は一意のシンボルとなるよう生成された新しい変数集合によって表される．NDSolveStateDataの"Variables"プロパティと"WorkingVariables"プロパティを使うと，これらの作業変数と指定された変数の間の対応が分かる．

作業変数と指定された変数との間の対応を表示する：

Wolfram Language code: Flatten[MapThread[Thread[Rule[##]]& , {state["Variables"], state["WorkingVariables"]}]]

明示的な常微分方程式系を生成する過程は，システムの記号的取扱いにより高価になる場合がある．このため，方程式を得るのに1秒という時間制限がある．この時間を超えると，NDSolveはメッセージを生成し，啓を微分代数方程式として解こうとする．

非線形の系の変数と方程式を定義する：

Wolfram Language code:

vars = {x, y, z};
eqns = { x'[t]  + y'[t] ^ 2 + z'[t] ^ 2 == Sin[t], x'[t] y'[t] + z'[t] == Sin[t] + x[t], x'[t] + z'[t] == y[t] + Cos[t]};
ics = {x[0] == 1, y[0] == -1, z[0] == 1};

デフォルトの方程式簡約化ストラテジーを使って解く：

Wolfram Language code: NDSolve[{eqns, ics}, vars, {t, 0, 1}]

次のメソッドオプションは，簡約化メソッドに指定して"EquationSimplification"の動作をよりよく制御することができる．

オプション名	デフォルト値
"TimeConstraint"	1	Solveが導関数を明示的に解くのに許された最大時間（秒）
"SimplifySystem"	False	できるだけ多くの従属変数について解析解を得ることで簡約化するかどうか

メッセージに示されているように，Solveが1秒というデフォルトの時間制限を超えていたのでNDSolveはこの導関数について明示的に解かなかった．NDSolveが明示的な式を得るのにかける時間は，オプション"TimeConstraint"で制御することができる．

オプション"TimeConstraint"を使うと，NDSolveが明示的な式を得るのにかける時間を秒単位で制御することができる：

Wolfram Language code: NDSolve[{eqns, ics}, vars, {t, 0, 1}, "Method" -> {"EquationSimplification" -> {Automatic, "TimeConstraint" -> 10}}]

上の例では，最初の方程式に二次項があるため，可能性のある4つの解が得られる．微分代数方程式として解く場合，どのように初期化されるかによって，数値解はこれらの解のいずれか一つに従う．

メソッドオプション"SimplifySystem"の使い方を示すために，次の例題を考える．

上記の系の解析解は，の解を求めるためにを最初の方程式に代入することで得られる．次にこれをの解を得るために第2の方程式に代入する．最終的な系の解は次のようになる．

NDSolveは，まず系の指数減少を実行して等価となる新しい常微分方程式系を形成してからでないと，この系を解くことはできない．これは高価になることを意味し得る．サブオプションの"SimplifySystem"を有効にすると，NDSolveは解析式または解析解が見付かる変数をすべて検出し，もとの系に繰り返し代入し続ける．この結果，もとの系は簡約化されるか，場合によっては（上記のように）解析解に戻る．

上の微分代数方程式を"SimplifySystem"->Trueで解く：

Wolfram Language code:

NDSolve[{x'[t] + y[t] == Sin[t], y'[t] + z[t] == Sin[t], x[t] == Cos[t]}, {x, y, z}, {t, 0, 1}, Method -> {"EquationSimplification" -> {Automatic, "SimplifySystem" -> True}}]

ある変数に対して解析解が存在しない場合は，NDSolveは解として補間関数を返す．

解析解が従属変数のうちのいくつかについてのみ見付けられる系を解く：

Wolfram Language code:

NDSolve[{x'[t] + y[t] == Sin[t], y'[t] + z[t] == Sin[t], y[t] == Cos[t], x[0] == 1}, {x, y, z}, {t, 0, 1}, Method -> {"EquationSimplification" -> {Automatic, "SimplifySystem" -> True}}]

サブオプションを有効にすると，定数パラメータは系に直接代入することができる．

系の一部として定数パラメータを含む常微分方程式を解く：

Wolfram Language code:

NDSolve[{y''[x] + Sin[a * y[x]]y[x] == b, y[0] == b, 
	y'[0] == a, a + b == 1, 2a - b == 3}, {y, a, b}, {x, 0, 30}, Method -> {"EquationSimplification" -> {Automatic, 
	"SimplifySystem" -> True}}]

微分代数方程式の解法

NDSolveには，微分代数方程式を解くためのさまざまな解法が組み込まれている．

指数1の微分代数方程式の一般的な残差形式には2つの方法が使える．

IDA—（SUNDIALSパッケージからの暗示的な微分代数ソルバ）後退微分公式(BDF)に基づく．

StateSpace—導関数について暗示的にを解き，内在する常微分方程式ソルバを使う．

系が指数2以上である場合，これらのソルバはどちらも失敗する．そのような系を正確に解くためには，指数減少法が必要である．指数1の系は常微分方程式に簡約できるが，上のソルバのどちらかを使った方が効率がよい．

NDSolveは，特殊形式に簡約することのできる微分代数方程式に使えるソルバも備えている．

MassMatrix—形式の微分代数方程式用

Projection—不変量の常微分方程式用

次のサブセクションでは，これらのメソッドを微分代数方程式に使うことについて説明する．

不変量のある常微分方程式

以下のような形式の微分代数方程式に出会うことはよくある．

（ここでは微分方程式に矛盾しない不変量）

微分代数方程式系が指数減少法によって常微分方程式に変換されると，常微分方程式を得るために微分された方程式は，その常微分方程式に矛盾しない．解は積分されるので，これらの方程式が満足されることを確実にすることは通常大切なことである．

このような系は未知数よりも多い方程式を持ち，制約条件が常微分方程式と矛盾していれば過剰決定系である．NDSolveは同数の方程式と従属変数を持つ系を想定しているため，そのような微分代数方程式は扱わない．そのような系を解くために，NDSolveに組み込まれている"Projection"メソッドが各ステップごとに計算された解をに投影することによって不変量を扱う．これにより，解が進化するにつれて確実に代数方程式は満足する．

このような制約条件付きの系の例には，振り子の動作をモデル化する非線形振動子がある．

振り子の動作のシミュレーションのために，方程式，不変制約条件，初期条件を定義する：

Wolfram Language code:

equation = x''[t] + Sin[x[t]] == 0;
invariant = x'[t]^2 - 2 Cos[x[t]];
start = {x[0] == 0, x'[0] == 1};

微分方程式は不変量の導関数であるため，この方程式を解く一つの方法は，不変量を使うことである．

メッセージが示すとおり，NDSolveはこの系を直接解かない：

Wolfram Language code: NDSolveValue[{equation, invariant == -1, start}, x, {t, 0, 1000}]

時間積分メソッド"Projection"を使うことで，時間刻みで不変量を満足する解を得る：

Wolfram Language code:

psol = NDSolveValue[{equation, start}, x, {t, 0, 1000}, Method -> {"TimeIntegration" -> {"Projection", "Invariants" -> invariant}}]

時間刻みにおける不変量をプロットする．これはほぼ機械精度まで満足する：

Wolfram Language code:

tsteps = First[InterpolatingFunctionCoordinates[psol]];ListPlot[Transpose[{t, invariant + 1 /. x -> psol} /. t -> tsteps]]

質量行列を持つ系を解く

導関数が線形の場合，微分系は次のように書くことができる．

ここでは質量行列と言われることの多い行列である．は減衰行列と呼ばれることもある．行列が非特異ならば，転置して，系を常微分方程式として解くことができる．しかし，特異行列がある場合は，系は微分代数方程式として解くことができる．どちらの場合も特殊形式を使うと利点が多い．

円柱の分子の動きを記述する以下の系

は，次のように質量行列で表すことができる．

系の変数，方程式，初期条件を定義する：

Wolfram Language code:

vars = {x[t], y[t], z[t]};
eqns = {x'[t] + y'[t] == y[t] - x[t], x'[t] - y'[t] == x[t] + y[t], z'[t] == 1 / (1 + t)};
init = {x[0] == 1, y[0] == 0, z[0] == 0};

オプション"MassMatrix"を使って，系を形式

で解く：

Wolfram Language code: sol = NDSolveValue[{eqns, init}, vars, {t, 0, 50}, Method -> {"EquationSimplification" -> "MassMatrix"}];

Wolfram Language code: ParametricPlot3D[Evaluate[sol], {t, 0, 50}, BoxRatios -> 1, PlotStyle -> Thick]

微分代数方程式の指数減少

NDSolveの微分代数方程式用の組込みソルバは現在完全に自動的に指数1の微分代数方程式を扱う．高指数の系では，解を得るために指数減少が必要である．この指数減少はNDSolveに適切なオプションを与えることで実行することができる．

NDSolveは記号的手法を使って指数減少を実行する．つまり，微分代数方程式系がという形式で表されるならば，ここでNDSolveはの要素がどのように相互作用を及ぼし合い，記号的導関数を計算するのかを見ることができず，指数減少は行えない．このためNDSolveはデフォルトでは指数減少を実行しないのである．

例として，振り子の動きを記述する古典的な微分代数方程式を考えてみる．

ここで紐の長さによって制約される点における質量がある．は実質的に紐の張力であるラグランジュ乗数である．指数減少の説明を簡単にするために，としよう．以下の図は振り子の系の略図である．

振り子に働く力の図

系の指数が方程式から明確でない場合，できることはNDSolveで試してみることである．ソルバで解けないときは，多くの場合指定された初期条件に対してどのような指数が現れるかを示すメッセージが生成される．

制約条件付きの振り子の系に対する方程式と初期値の完全な集合を定義する：

Wolfram Language code:

cpdae = {x''[t] + λ[t] x[t] == 0, y''[t] + λ[t] y[t] == -1, x[t] ^ 2 + y[t] ^ 2 == 1};
cpinit = {x[0] == 1, x'[0] == 0, x''[0] == 0, y[0] == 0, y'[0] == 0, y''[0] == -1, λ[0] == 0};

制約条件付きの振り子の系を解いてみる：

Wolfram Language code: NDSolve[{cpdae, cpinit}, {x, y, λ}, {t, 0, 100}]

このメッセージは指数が3であるということを示し，系を解くためには指数減少が必要であるといっている．ここでは，初期化で起こり得る問題を避けるために，完全に矛盾のない初期化が使われる．解法のこの部分は，「微分代数方程式の矛盾のない初期化」で扱う．

NDSolveが，一般的な指数減少オプション iropts で指数減少法 irmeth を使うように指定するためには，オプション設定Method->{"IndexReduction"->{irmeth,iropts}}を使う．

None	指数減少なし（デフォルト）
Automatic	自動的にメソッドを選択する
"Pantelides"	グラフベースのPantelidesアルゴリズム
"StructuralMatrix"	構造行列ベースのアルゴリズム

指数減少法

指数減少法すべてに次のオプションを与えることができる．

"ConstraintMethod"	Automatic	指定された系の制約条件をどのように扱うか
"IndexGoal"	Automatic	減少させる指数（1か0）

指数減少法のオプション

指数減少に自動的に選ばれたメソッドを使って解く：

Wolfram Language code:

sol = First[NDSolve[{cpdae, cpinit}, {x, y, λ}, {t, 0, 100}, Method -> {"IndexReduction" -> Automatic}, MaxSteps -> ∞]]

Wolfram Language code: Plot[Evaluate[{x[t], y[t]} /. %], {t, 0, 100}]

指数減少は，微分代数方程式系の方程式を微分することによって実行される．deqn=eqn が eqn の代りに系に含まれるように，指数減少法の過程で，方程式 eqn が微分されるとする．微分が終了すると，微分された方程式は完全決定系を構成し，自分で解くことができる．しかし，系の中に制約条件としてもとの方程式を含ませるのは通常大切なことである．これは指数減少オプションの"ConstraintMethod"で制御される．

None	もとの方程式を制約条件として維持せず，差分方程式だけで解く
"DummyDerivatives"	MattsonとSoderlind [MS93]の方法を使い，いくつかの導関数を代数変数（ダミー導関数）で置き換える
"Projection"	指数を0に減少させて常微分方程式を得，不変量として微分されたもとの方程式で"Projection"時間積分メソッドを使う

制約条件メソッド

指数減少法を使い，微分された方程式で解く：

Wolfram Language code:

ucsol = First[NDSolve[{cpdae, cpinit}, {x, y, λ}, {t, 0, 100}, Method -> {"IndexReduction" -> {Automatic, "ConstraintMethod" -> None}}, MaxSteps -> ∞]];

制約条件

が2つの解をどれだけ満足しているかを比較する：

Wolfram Language code: LogPlot[Evaluate[Abs[x[t] ^ 2 + y[t] ^ 2 - 1] /. {sol, ucsol}], {t, 0, 100}]

微分された方程式だけの数値解では，振り子の長さが1から離れていく．もとの方程式を制約条件として解に含ませることは，指数減少の大切な一面である．

次の2つのセクションでは指数減少アルゴリズムと，NDSolveで使うことのできる制約メソッドについて説明する．

指数減少法

NDSolveには指数減少を行うアルゴリズムとしてPantelides法と構造行列法の2つがある．これらの方法は両方とも方程式の記号形式を使ってどの方程式を微分するかを決定し，記号的微分を使って微分された系を得る．

どちらの方法も方程式の構造的な接続の概念を使っている．特に，変数が方程式に明示的に現れたら，はのtermed incidentである．構造的な接続行列とは，がの接続ならば，要素を持ち，それ以外の場合はを持つ行列である．対角要素に1が存在するように構造的接続行列を並べ替えることができるならば，微分代数方程式系は構造的指数1を持つ．このような並替えができるということは，各変数に対して1つの方程式があるように，方程式と変数の間の対応が存在するということを意味している．

Pantelides法は，非常に大きな系でも効率的であり得る接続に基づいたグラフに使える．

構造的な指数は，その系の実際の指数よりは小さい場合があることに注意する．実際の系では，構造的に存在する項がいくつかの解とともに消滅するか，ヤコビ行列が特異であるかする可能性がある．構造行列法は，ヤコビアンに特異性がある可能性を考慮し，系によってはよりよい指数減少を実行するが，計算費用は大きくなる．

次のセクションでは，Pantelides法の説明に便利な構造的接続行列を得るコマンドを紹介し，その次のセクションでは構造行列法について説明する．

構造的接続行列

NDSolve`StructuralIncidenceArray[eqns,vars]

方程式 eqns の変数 vars に対する構造的接続行列のパターンを持つSparseArrayを返す

構造的接続行列を得る

変数

，

に対する系

，

の構造的接続行列を得る：

Wolfram Language code: S = NDSolve`StructuralIncidenceArray[{x'[t] == y[t], y[t] == 2 x[t]}, {x'[t], y[t]}]

Wolfram Language code: MatrixForm[S]

SparseArrayはメモリ使用量が少なくてすむ非零のパターンだけを含んでいる．これは接続のパターンを表すことを意味する．ArrayRulesを使うと，1のSparseArrayに変換することができる．

パターン配列を1と0に変換するコマンドを定義する：

Wolfram Language code: convertPatternArray[S_SparseArray] := SparseArray[Drop[ArrayRules[S], -1] /. Verbatim[_] -> 1, Dimensions[S]]

Wolfram Language code: ArrayRules[S]

Wolfram Language code: S1 = convertPatternArray[S]

Wolfram Language code: MatrixForm[S1]

ここでは対角要素の零を避けるのに並替えが必要ないので，この系は構造的指数1を持つ．

導関数変数についての構造的並替えで，対角要素に零を避けるよう並べ替えられるならば，この系は指数0である．

変数

，

に対する系

，

の構造的接続行列を得る：

Wolfram Language code: MatrixForm[NDSolve`StructuralIncidenceArray[{x'[t] == y[t], y[t] == 2 x[t]}, {x'[t], y'[t]}]]

対角要素で零を避ける並替えはない．しかし，2つ目の方程式が微分されるとこれは可能になる．

2つ目の方程式を変数

，

について微分して，系

，

の構造的接続行列を得る：

Wolfram Language code: MatrixForm[NDSolve`StructuralIncidenceArray[{x'[t] == y[t], D[y[t] == 2 x[t], t]}, {x'[t], y'[t]}]]

微分された系が構造的指数0を持つので，もとの系は指数1であることが証明される．

もう一つの例として制約条件付きの振り子(2)を考える．

現れる最高階数の導関数について，制約条件付きの振り子の系の構造的接続行列を得る：

Wolfram Language code: MatrixForm[NDSolve`StructuralIncidenceArray[cpdae, {x''[t], y''[t], λ[t]}]]

この行列は対角要素の零要素を避けるために並べ替えることができないので，指数は1より大きい．

関数NDSolve`StructuralIncidenceArrayは接続をチェックするために，式のFullFormで文字のマッチングを行うだけである．この関数は，方程式や変数が適切な形であるかどうかのチェックはしない．

次は，それぞれの方程式の中のシンボル x および y の文字の存在をテストする：

Wolfram Language code: MatrixForm[NDSolve`StructuralIncidenceArray[{x'[t] == y[t], y[t] == 2 x[t]}, {x, y}]]

x'[t]のFullFormが x を含むDerivative[1][x][t]なので，これは完全行列である．

Pantelides法

Pantelidesによって提唱された方法[P88]は，もとは微分代数方程式の矛盾のない初期化を求めるために提唱されたグラフ理論的なメソッドである．これは従属変数と方程式の二部グラフに使うことができ，アルゴリズムがグラフを走査すると，系は最大で指数1に簡約される．ここでの走査とは，指数行列の対角線上にゼロがないような変数と方程式の順序を指す．

二部グラフに完全走査がない場合，アルゴリズムは方程式を微分することで経路を拡大し，その過程で新しい変数（前の変数の導関数）を導入する．もとの系が構造的に特異でなければ，アルゴリズムは走査で終了する．通常アルゴリズムは走査に必要な微分だけを行うが，これは貪欲法なので，常に最小数とは限らない．

このアルゴリズムの詳細は，このドキュメントでは割愛する．Wolfram言語に組み込まれている実装は，[P88]で説明されているアルゴリズムに密接に従っている．これはGraphを使ってグラフ計算を実装し，微分が求められたときはDを使って方程式を記号的に微分する．

系の走査があるときは，指数行列がブロック下三角形式になるような変数と方程式の順序を見付けることが可能である．ブロック下三角形式は制約条件を維持するためのダミー導関数メソッドの設定において，また"StateSpace"時間積分メソッドにおいても使われる．ブロック下三角行列の順序に関する詳細は，「微分代数方程式の状態空間法」をご覧いただきたい．

振り子系の指数減少法(2)を始めるために，変数を考える．通常，方程式に現れる最高階の導関数から開始する．振り子方程式は以下の通りである．

また，最高階数の導関数についてのグラフ接続行列は以下である．

接続行列は方程式の変数の存在をチェックすることによって構築される．方程式に変数が存在する場合，重み1が与えられ，存在しない場合は0が指定される．したがって，最初の方程式の変数のチェックでが得られる（接続行列の最初の行に見られるとおりである）．

最初の接続行列から，走査がないことが分かる．最後の方程式は2回微分される．

これで接続行列は次のようになる．

これは，結果として接続行列になる最初と最後の方程式を入れ替えることによって対角要素に0がないように並べ替えることができる．

結果の系の指数は1（がまだ代数的変数として現れるため）であり，微分された系に導関数を含まない．もう一度系を微分することによっての微分系が与えられる．ゆえに，もとの系の指数は3となる．

最初のアプローチとして2回微分された方程式を持つ系は直接解くことができる．

2度微分された方程式を含む系を解く：

Wolfram Language code:

d2sol = NDSolve[{x''[t] + λ[t] x[t] == 0, y''[t] + λ[t] y[t] == -1, D[x[t] ^ 2 + y[t] ^ 2 == 1, {t, 2}], cpinit}, {x, y, λ}, {t, 0, 100}, MaxSteps -> ∞]

最初に指定された方程式に基づいた不変量である数量をプロットする：

Wolfram Language code: Row[Table[Plot[Evaluate[D[x[t] ^ 2 + y[t] ^ 2 - 1, {t, i}] /. d2sol], {t, 0, 100}, ImageSize -> Small], {i, 0, 1}]]

もとの振り子の長さの制約条件はあまりよく満足されておらず，実際には伸びている．新しい系がもとの系を方程式を完全に等しく表していないので，制約条件は十分満足しないのである．これは新しい系にもとの制約条件ではなく，微分された方程式が含まれているという事実から分かる．

新しい系を持つもとの系の力学をモデル化するためには，追加の方程式が満足されなければならない．しかし，追加の方程式を含めることは未知数よりも方程式が多くなる．この問題を解決するために，指数減少された系が平衡となるよう，ダミー導関数の方法が使われる．

Pantelidesアルゴリズムは，非常に大きい問題に対してもよく制御された複雑性を備えたグラフアルゴリズムを使うことができるので，効率的である．しかし，これは事象に基づいているため，特異のヤコビアンになる系では問題があるかもしれない．構造行列法ならばこのような場合でも使えるかもしれないが，大きい系では計算複雑度が大きくなる可能性がある．

構造行列法

構造行列アルゴリズムはPantelidesアルゴリズムの代りとなる．構造行列法はUnger [UKM95]とChowdhary et al. [CKL04]の研究に従っている．Pantelidesアルゴリズムのようなグラフベースのアルゴリズムは指数減少法を実行するのに走査に依存するが，特定の系では変数の打消しがある場合があるということを考慮しない．そのため，アルゴリズムは系の指数を少なく見積ることになる．微分代数方程式が指数0あるいは1の系に正しく簡約されないならば，数値積分は失敗するか不正確な結果を出す可能性がある．

Pantelides法で考慮されない変数打消しの例として，次の微分代数方程式を考える．

指数減少法を実行する最初のステップとして，2つ目の方程式を1度微分すると，以下の新しい系が得られる．

Pantelides法は最初の微分の後に変数の導関数と走査を見付けるので，そこで中断する．最初の微分の後すべての導関数が見付かるので，Pantelides法は系の指数を0と推定する．指数0の系は常微分方程式系と等価である．微分された系の導関数のヤコビアンは以下の

であり，これは明らかに特異行列である．代入すると，微分された系は次と等価であることが分かる．

これは，常微分方程式系を得るためには，2つ目の方程式に2回目の微分を行う必要があることを意味している．これで最終的に以下の系が得られる．

この系は導関数に非特異ヤコビアンを持つ．これは正確に指数0の系になっており，これから常微分方程式系を構成することができる．2回の微分が必要だったので，この系の実際の指数は1ではなく2である．

Pantelides法と異なり，構造行列法は線形項に関連する打消しすべてを考慮する．この方法のステップは，一階の系として記述された振り子の例(2)を使って示すことができる．

まず2つの行列およびが構築される．これらはそれぞれ導関数と変数に関連付けられた接続行列である．行列は以下で与えられる．

上の行列の項目としてを含むものがある．は方程式の変数（行）が非線形の方程式で表示されていることを示すプレースホルダとして使われている．2つ目の方程式では，項とは積として現れている．それゆえ，2行目の最初と最後の列にの印がついているのである．

行列が構築されると，その目的は微分代数方程式系を常微分方程式系に変換することになる．そのためには，行列は完全階数でなければならない．これを実現するために，まず微分する必要のある方程式が見付けられる．この場合は最後の方程式である．微分の後，構造行列は次のようになる．

行列の最後の行は微分したために変化した．これをもっとはっきりと見るために，(3)の最後の方程式を考える．方程式を微分するとが得られる．これで項が方程式に存在するようになり，非線形で表されている．したがって，導関数に関連付けられた接続行列の行はであり，変数に関連付けられた接続行列の行はとなる．

次のステップでは，行列に対して，ガウス消去法の形式での行列因子分解が関わる．の最終行の要素は，行1と3を使うことで消去できる．この消去は行列の行にも影響を及ぼす．の最初の行にを掛け，それに最後の行を代入すると次が得られる．

3列目に*を掛け，それを最後の列から引くと次のようになる．

行列はまだ階数落ちであり，上記のステップを再び繰り返さなければならない．2回目の繰返しで以下の系が得られる．

3回目の反復では次が得られる．

行列はこれで完全階数の行列になったので，ここで反復を終了する．の完全階数行列を得るために3回の反復が必要であったので，この系の指数は3である．Pantelides法とは異なり，構造行列法は接続行列に明示的に作用する．

次の例ではPantelides法よりも構造行列法が優れている点，および両者の差を示す．次の線形系を考える．

方程式と初期条件を設定する：

Wolfram Language code:

eqns = {x'[t] + y'[t]== y[t] -Sin[t] + 2 Cos[t], x[t]+y[t]==-Sin[t]};
ic = {y[0]==-3};

Pantelides法を使って系を解いてみる：

Wolfram Language code:

NDSolve[{eqns, ic}, {x, y}, {t, 0, 1}, Method -> {"IndexReduction" -> "Pantelides", "EquationSimplification" -> "Residual"}]

これが失敗するのは，特異のヤコビ行列があるためである．微分代数方程式の指数が1より大きい場合，メッセージNDSolve::nderrあるいはNDSolve::ndcfが出力されるのは珍しくない．

指数減少に構造行列法を使ってこの系を解く：

Wolfram Language code: sol = NDSolve[{eqns, ic}, {x, y}, {t, 0, 2 Pi}, Method -> {"IndexReduction" -> "StructuralMatrix"}]

厳密解はとである．厳密解と数値解を比較すると，計算された解が非常に正確だということが分かる．

Wolfram Language code:

GraphicsRow[Plot[Evaluate[# /. sol], {t, 0, 2 Pi}]& /@ {Norm[x[t] - (3 Cos[t] - 2 Sin[t])], Norm[y[t] - (Sin[t] - 3 Cos[t])]}, ImageSize -> 500]

系の指数を過小／過大評価した際の影響を考えるため，以下の線形系を考える．

方程式，初期条件，変数を設定する：

Wolfram Language code:

Clear[x1,x2,x3,x4,x5,x6,ics,vars];

eqns = (⁠|                                                        |
| ------------------------------------------------------ |
| x1'[t] == Cos[t] - x2[t] + x3[t]                        |
| x2'[t] == Cos[t] - x3[t] + x4[t]                        |
| x3'[t] == Cos[t] - x4[t] + x5[t]                        |
| x4'[t] == Cos[t] - x5[t] + x6[t]                        |
| x5'[t] == Cos[t] + Sin[t] - x6[t]                       |
| 0 == x1[t] + x2[t] + x3[t] + x4[t] + x5[t] - 5 * Sin[t] |⁠);

ics = {x1[0]==x2[0]==x3[0]==x4[0]==x5[0]==x6[0]==0};
vars = {x1,x2,x3,x4,x5,x6};

すべての変数の厳密解はである．これは指数4の系である．Pantelides法は指数を2だと推定している．

指数減少にPantelides法を使って系を解いてみる：

Wolfram Language code: NDSolve[{eqns, ics}, vars, {t, 0, 4 Pi}, Method -> {"IndexReduction" -> "Pantelides"}] ;

ヤコビアンが特異であるため，ソルバは収束の失敗に遭遇する．

指数減少に構造行列法を使って系を解いてみる：

Wolfram Language code: solS1 = NDSolve[{eqns, ics}, vars, {t, 0, 4 Pi}, Method -> {"IndexReduction" -> "StructuralMatrix"}]

どちらの指数減少法を使ったらよいかを決めるためにNDSolveが実行するテストがある．この場合，オプションAutomaticを使うと，NDSolveは構造行列法を選ぶ．

"IndexReduction"->Automaticに設定して系を解く：

Wolfram Language code: solS2 = NDSolve[{eqns, ics}, vars, {t, 0, 4 Pi}, Method -> {"IndexReduction" -> Automatic}]

もし，この系が指数0として扱われるよう強制されるなら，積分を実行することができる．しかし指数の推定が不正確ではない場合，解は全く異なる多様体について計算されるので，結果は大きくずれる可能性がある．

Pantelidesアルゴリズムを使って，この系を指数0の系として解く：

Wolfram Language code:

solP0 = NDSolve[{eqns, ics}, vars, {t, 0, 4 Pi}, Method -> {"IndexReduction" -> {"Pantelides", "IndexGoal" -> 0}, "EquationSimplification" -> "Residual"}]

この問題の厳密解は既知なので，数値解と厳密解の違いを見るのが望ましい．

Wolfram Language code:

Row[Plot[Evaluate[(#[t] - Sin[t]) /. solP0], {t, 0, 4 Pi}, PlotLabel -> #, PlotRange -> All, ImageSize -> Small]& /@ vars]

系の指数を低く見積ったため，必要な微分が実行されず，解は誤った多様体上で近似される．一方，構造行列法は系の指数を正しく見付けることができる．

"StructuralMatrix"を使って得られた系の結果をプロットする：

Wolfram Language code:

Row[Plot[Evaluate[(#[t] - Sin[t]) /. solS1], {t, 0, 4 Pi}, PlotLabel -> #, PlotRange -> All, ImageSize -> Small]& /@ vars]

構造行列法を使ってこの系を指数0の系として解く：

Wolfram Language code:

solSM0 = NDSolve[{eqns, ics}, vars, {t, 0, 4 Pi}, Method -> {"IndexReduction" -> {"StructuralMatrix", "IndexGoal" -> 0}}];

Wolfram Language code:

Row[Plot[Evaluate[(#[t] - Sin[t]) /. solSM0], {t, 0, 4 Pi}, PlotLabel -> #, PlotRange -> All, ImageSize -> Small]& /@ vars]

構造行列法は，このような系をずっとうまく扱う．しかし，打消しに対して支払う代償もあることに注意することが大切である．構造行列法は，Pantelidesアルゴリズムと比べると格段の計算オーバーヘッドを有する行列の生成，維持，操作に依存する．したがってこのアルゴリズムは小規模から中規模の問題にのみ適している．

制約法

ダミー導関数

ダミー導関数の目的は，もとの方程式と指数減少法で微分された方程式を組み合せた系が過剰決定系にならないように代数的変数を導入することである．

振り子の系(2)を再び見てみる．

次はNDSolveに組み込まれた，ダミー導関数を含む指数減少法を使って系を解く：

Wolfram Language code:

ddsol = First[NDSolve[{cpdae, cpinit}, {x, y, λ}, {t, 0, 100}, MaxSteps -> ∞, Method -> {"IndexReduction" -> {Automatic, "ConstraintMethod" -> "DummyDerivatives"}}]]

この例では，"ConstraintMethod"->"DummyDerivatives"が強調のために含まれているが，ダミー導関数はデフォルトなので実際には必要ない．

もともと指定された方程式に基づいた不変量である数量をプロットする：

Wolfram Language code:

GraphicsRow[Table[Plot[Evaluate[D[x[t] ^ 2 + y[t] ^ 2 - 1, {t, i}] /. ddsol], {t, 0, 100}, PlotRange -> All], {i, 0, 1}], ImageSize -> 500]

これで制約条件は十分満足された．

MattsonとSöderlindのダミー導関数法の考え方[MS93]は，導関数を表す新しい変数を導入することで，過剰決定系を得ることなく制約方程式を再び導入することができるようにすることである．

系(2)に対して指数減少された方程式を思い出してみる．

ここでは，3つ目の方程式からとその最初の導関数からの2つの制約条件がある．

とをそれぞれ，に置き換わる代数的変数とする．制約条件が含まれると，これらの変数の系全体は以下で与えられる．

また，変数の接続行列は次のようになる．

これは非零対角成分を持つよう並べ替えることができる．

接続行列を並べ替える：

Wolfram Language code:

incidenceMatrix = {{1, 0, 1, 0, 1}, {0, 0, 0, 1, 1}, {1, 1, 1, 1, 0}, {1, 1, 0, 0, 0}, {1, 0, 0, 0, 0}};
{row, col} = Flatten /@ SparseArray`BlockTriangularOrdering[incidenceMatrix];
incidenceMatrix[[row]][[col]] // MatrixForm

一般にアルゴリズムはどの指数減少された微分代数方程式にも適用することができるが，系を積分する間に，ソルバが特異性に遭遇する可能性がある．

と

を，手動で導入したダミー導関数で置き換え，系を積分してみる．ダミー導関数を使うことはデフォルトなので，オプション設定は必要ない：

Wolfram Language code:

ddxsol = First[NDSolve[{xd2[t] + x[t] λ[t] == 0, y''[t] + y[t] λ[t] == -1, 2 xd1[t]^2 + 2 xd2[t] x[t] + 2 Derivative[1][y][t]^2 + 2 y[t] Derivative[2][y][t] == 0, 2 xd1[t] x[t] + 2 y[t] Derivative[1][y][t] == 0, x[t]^2 + y[t]^2 == 1, cpinit /. {x' -> xd1, x'' -> xd2}} , {x, xd1, xd2, y, λ}, {t, 0, 100}, MaxSteps -> Infinity]]

上の系はのときに特異になるので，積分を終了することができない．

以下は計算された分だけ解をプロットする：

Wolfram Language code: Plot[Evaluate[{x[t], y[t]} /. ddxsol], {t, 0, 1.85}]

このプロットから，が0に近付くとNDSolveが動かなくなることが分かる．が0のとき，いくつかの項がなくなり，接続行列は特異になり，非零の対角要素は見付からなくなる．そのため，ダミー導関数の系は実質的にを通過する解に沿っての指数を持つ．

ダミー導関数の別の置換，つまりとをとで置き換えることを考えてみる．この置換を持つ系はで特異になる．

この場合，が0から離れているときはとをとで置き換えるような動的状態選択を使い，が小さすぎるようになったらとの代りにとを使って動的にこの系に切り換えるようにするとよい．NDSolveは必要に応じて，ダミー導関数を使った動的状態選択を扱う状態置換規則で自動的にWhenEventを使う．動的状態選択の詳細はこのチュートリアルでは触れない．詳細はMattsonとSoderlindの論文[MS93]を参照のこと．

投影法

ダミー導関数を使う以外の方法として指数を0に減少させて常微分方程式系を得て，もとの制約条件が確実に満足するよう，投影法を使うというものがある．

振り子の系(2)では，系の指数は3である．つまりすべての変数に対する常微分方程式を得るために，最初と2番目の方程式が1回ずつ微分され，3つ目の方程式は3回微分されなければならないということである．この結果，以下が得られる．

これは導関数，，について解くことができる（NDSolveは系を自動的に一階まで簡約する）．

系の力学を正確に表現するためには，最初の回微分された方程式も考慮に入れる必要がある．これらの方程式は下の通りである．

これらの方程式は各時間刻みで満足しなければならない不変方程式として扱われる．常微分方程式はこれらの方程式を微分することにより導出されたので，制約条件は常微分方程式と矛盾のないことが保障されており，制約条件を不定量として"Projection"時間積分メソッドを使ってこの系を解くことができる．

"Projection"メソッドを使って常微分方程式を解く：

Wolfram Language code:

NDSolve[{λ[t] Derivative[1][x][t] + x[t] Derivative[1][λ][t] + Derivative[3][x][t] == 0, λ[t] Derivative[1][y][t] + y[t] Derivative[1][λ][t] + Derivative[3][y][t] == 0, 6 Derivative[1][x][t] Derivative[2][x][t] + 6 Derivative[1][y][t] Derivative[2][y][t] + 2 x[t] Derivative[3][x][t] + 2 y[t] Derivative[3][y][t] == 0, cpinit}, {x, y, λ}, {t, 0, 100}, Method -> {"TimeIntegration" -> {"Projection", "Invariants" -> Apply[Subtract, Flatten[{x''[t] + λ[t] x[t] == 0, y''[t] + λ[t] y[t] == -1, Table[D[x[t] ^ 2 + y[t] ^ 2 == 1, {t, i}], {i, 0, 2}]}], {1}]}}]

これはNDSolveに組み込まれている指数減少法を直接使って自動的に行うことができる．

これはNDSolveに組み込まれている投影法による指数減少法を使って系を解く：

Wolfram Language code:

psol = First[NDSolve[{cpdae, cpinit}, {x, y, λ}, {t, 0, 100}, MaxSteps -> ∞, Method -> {"IndexReduction" -> {True, "ConstraintMethod" -> "Projection"}}]];

このメソッドでは，もとの制約条件はオプションPrecisionGoalとAccuracyGoalによって指定されたNDSolveの局所的許容誤差まで満足する．

この系には満足しなければならない不変量がまだ他にある．そのうちの一つはで表すことができるエネルギー保存である．投影法は時間刻みの最後だけで投影するので，比較をするのに最も適している．

次は時間刻みにおける解に対するエネルギーを得て，ダミー導関数および投影解のエネルギーをプロットする：

Wolfram Language code:

ddsol = First[NDSolve[{cpdae, cpinit}, {x, y, λ}, {t, 0, 100}, MaxSteps -> ∞, Method -> {"IndexReduction" -> {True, "ConstraintMethod" -> "DummyDerivatives"}}]];

Wolfram Language code:

energy[sol_] := Module[{times = Part[x /. sol, 3, 1], es}, 
	es[t_] = 1 + y[t] + (1/2) (-y[t] Derivative[1][x][t] + x[t] Derivative[1][y][t])^2 /. sol;Transpose[{times,   es[times] - es[0]}]];

Wolfram Language code: ListPlot[{energy[ddsol], energy[psol]}]

偏微分代数方程式の指数減少

以下の偏微分方程式系を考える．

この系では2番目の偏微分方程式には，どの変数に対しても時間導関数要素が含まれていない．したがって，結果となる質量行列は特異行列となり，系は指数1を持つため，NDSolveはこの系を線の方法で解くことはできない．この系の指数を減少させるためには，上の系は次のような行列ベクトル形式で表される．

項は空間導関数を離散化することで得られた行列を，項は変数を離散空間区間で離散化することで得られたベクトルを表す．上記の指数減少法を使って新しい系に指数減少を施す．指数減少された結果の系は以下の通りである．

空間導関数は再び系に導入されて以下のようになる．

結果となる系はNDSolveで利用できる従来の方法を使って解くことができる．

系を定義する：

Wolfram Language code:

eqns = {D[u[x, t], t] - D[v[x, t], t] == D[u[x, t], x, x] + u[x, t] D[v[x, t], x] - v[x, t] D[u[x, t], x] + (2 Pi t) ^ 2 v[x, t], D[v[x, t], x, x] + (2 Pi t) ^ 2 u[x, t] == 0};
bcs = {u[0, t] == v[0, t] == 0, u[1, t] == Sin[2 Pi t], Derivative[1, 0][v][1, t] == 2 Pi t Cos[2 Pi t]};
ics = {u[x, 0] == v[x, 0] == 0};

偏微分代数方程式に指数減少を実行することにより系を解く：

Wolfram Language code:

sol  = NDSolveValue[{eqns, bcs, ics}, {u, v}, {x, 0, 1}, {t, 0, 3}, Method -> {"IndexReduction" -> Automatic, "EquationSimplification" -> "Residual", "PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 200, "MaxPoints" -> 200}}}]

定数初期条件に基づいた空間グリッドのデフォルトスペースでは，解が時間に伴って生成する変形を扱うには不十分であるため，空間離散には200個の点が使われたことに注目のこと．

解をプロットする：

Wolfram Language code: Plot3D[Evaluate[#[x, t]], {x, 0, 1}, {t, 0, 3}]& /@ sol

指数減少法は時間導関数で実行されるので新しい境界条件は必要なく，系に加える必要もないという点に注意する．しかし，偏微分代数方程式が高階ならば，追加の初期条件が必要になる場合がある．

微分代数方程式の矛盾のない初期化

微分代数方程式系は最も一般的な形式

で表すことができる．これには微分方程式と代数的制約条件を含むことがある．の解を得るためには，積分を開始するためのおよびの矛盾のない初期条件集合が必要である．一貫性に対して必要な条件は初期条件がを満足するということである．しかし，この条件はだけでは十分ではない．もとの方程式を微分すると，また初期条件により満足される必要のある新しい方程式ができるだけだからである．したがって，必要なことは，必要となる矛盾のない条件すべてを満足する初期条件を見付けることである．矛盾のない初期条件を見付けるという問題は，微分代数方程式系を解く上で最も難しい部分の一つであり得る．

次の線形微分代数方程式問題を考える．

微分代数方程式系はについての導関数を持つ．これが意味するのは，この問題を一意に解くためには，に対する初期条件を指定しなければならないということである．しかし，矛盾のない一意の初期条件集合があるため，変数の初期条件を任意に指定することはできない．代数方程式を1回微分するとが得られる．これを2つ目の方程式に代入するとの解が得られる．を微分して最初の方程式に代入するとの解が得られる．これは，この微分代数方程式系の解が決まっているため初期条件も決まっており，に対する矛盾のない初期条件の唯一の集合はであることを意味する．初期条件の制約条件は一般に明らかではないので，矛盾のない初期条件を見付けるという問題が非常に困難になるのである．

場合によっては，微分代数方程式問題の矛盾のない初期化だけを求めて，解は後で計算することもあろう．これをNDSolveで行う場合は，時間積分の数点を初期条件が指定されている時間と同じにすることで実行できる．このようなステップは，完全な積分を実行する必要がなく，初期条件だけを調べたいときに便利である．

次は微分代数方程式系の初期条件を計算する：

Wolfram Language code:

init = NDSolve[{Derivative[1][Subscript[x, 1]][t] == Subscript[x, 3][t], Subscript[x, 2][t] (1 - Subscript[x, 2][t]) == 0, Subscript[x, 1][t] Subscript[x, 2][t] + Subscript[x, 3][t] (1 - Subscript[x, 2][t]) == t, Subscript[x, 2][2] == 1}, {Subscript[x, 1], Subscript[x, 2], Subscript[x, 3]}, {t, 2}]

Wolfram Language code: First[{Subscript[x, 1][t], Subscript[x, 2][t], Subscript[x, 3][t]} /. init /. t -> 2]

NDSolveは微分代数方程式の初期化にいくつかの方法を提供する．メソッド m はMethod->{"DAEInitialization"->m}を使って指定することができる．

Automatic	使用するメソッドを自動的に決定する
"Collocation"	選点法を使う．このアルゴリズムは指数の高い系の初期化に使うことができる
"QR"	指数1の系にQR分解ベースのアルゴリズムを使う
"BLT"	指数1の系にブロック下三角アプローチを使う

微分代数方程式の初期化メソッド

選点法は残差のブラックボックスとして微分代数方程式系を扱うよう設計されており，初期化の点付近の小さい区間上で残差を満足するという目的と達成しようと試みる．この機能により，高指数の系の初期化を扱うためにこのメソッドを使うことができる．しかし，アルゴリズムは微分代数方程式系を分析しないので，系の構造は計算の効率性のためにはつかない．このためこのアルゴリズムは他の方法より遅いのである．

QR法とBLT法は，特に指数1と指数0の系に使うために設計されている．QRメソッドは変数とその導関数の初期条件が反復的に計算されるよう，2つの分離した部分系を生成するために導関数のヤコビアンを分解することに依存する．これによりこのメソッドは非常に効率よくロバストになっている．

BLT法は微分代数方程式系の構造を検査することに依存しており，もとの系を多数の小さい部分系に分割してその各々が効率よく解けるようにする．このメソッドは全体の系と比較して，ずっと小さいヤコビ行列を扱うことになり，そのため大規模の系には計算的に効率がよい．

以下のセクションでは，NDSolveが異なるアルゴリズムを使ってどのように代数微分方程式の矛盾のない初期条件を得るかを説明する．

選点法

基本的な選点法は基底関数に関する従属変数を

として展開することを利用するものである．とを求めるために，が個の選ばれた点において満足するという条件を実行しようとする．これを行うために，陰的な微分代数方程式は，

として線形化される．ここでとはそれぞれとに関するのヤコビアンを表す．

時間内にこの線形化を個の選点に適用すると，係数に対する線形方程式系が得られる．これはニュートン法を使って反復的に解かれる．

ここでは直線探索法を使って得られる．個の係数と個の選点がある場合，(5)の × の線形方程式系を扱うことになる．

指数減少を行っていない高指数の系の場合，NDSolveは自動的に選点アルゴリズムを使って系を初期化しようとする．振り子の系(2)を考えてみる．

系の積分を開始するためには，の初期条件が分かっていなければならない．このような系の場合，はある物理的特性を表し，合理的な値を割り当てることができる．しかしの値は見付けやすくはないかもしれない．ここで適切な値を決めるのに，NDSolveを使うことができる．

振り子の微分代数方程式に対する矛盾のない初期条件集合を得る：

Wolfram Language code:

sol = First@NDSolve[{x''[t] + λ[t] x[t] == 0, y''[t] + λ[t] y[t] == -1, x[t] ^ 2 + y[t] ^ 2 == 1, x[0] == 1, y'[0] == 1}, {x, y, λ}, {t, 0, 0}]

Wolfram Language code: {x[0], y[0], λ[0]} /. sol

Wolfram Language code: {x'[0], y'[0], λ'[0]} /. sol

部分的な初期条件だけが与えられると，アルゴリズムが必要な初期条件すべてを正しく計算しているのに注目のこと．

選点法のデフォルト設定は，多様な系の初期化を扱うよう選ばれているが，場合によっては矛盾のない初期化を見付けるために設定を変更する必要がある．アルゴリズムを調整するのには以下のオプションが使える．

選点法のオプション名	デフォルト値
"CollocationDirection"	Automatic	選点が置かれる方向を選ぶ．可能な値は"Forward"，"Backward"，"Centered"，"Automatic"である
"CollocationPoints"	Automatic	使用される選点の数．サブオプションの"ExtraCollocationPoints"を使うと，近似の順序を固定したまま格子点が加わる
"CollocationRange"	Automatic	選点が分布する範囲
"DefaultStartingValue"	Automatic	非線形の反復における指定されていない要素に使用する初期値．可能な値は1つの要素の値 ig/{ig}，あるいは2つの数のリスト{ig1,ig2}である．ここで ig1,ig2はそれぞれ変数の初期値とその導関数である
"MaxIterations"	100	実行される繰返しの最大回数

"Collocation"メソッドオプション

次のセクションでは，これらのオプションを詳しく説明し，これらがどのような状況で矛盾のない初期条件を見付けるのに役立つかを示す．

選点の方向

初期化を実行するためには，選点の集合が選ばれる．点はオプション"Centered"，"Forward"，"Backward"を使うことで，関心点の両側，前，後ろのどこに置くこともできる．次の表はこれらのオプションを選んだときの効果を示すものである．

"Forward"	範囲上で実行される選点法
"Backward"	範囲上で実行される選点法
"Centered"	範囲上で実行される選点法

"CollocationDirection"サブオプション

上記の項は初期条件が計算される時間を指している．項はオフセット値である．デフォルトでは，NDSolveは自動的に選点法を実行する方向を選択する．

選点の方向は，不連続性を扱う際に重要になる．

Wolfram Language code:

DAE2 = (|                       |
| --------------------- |
| x'[t]   == Abs[u[t]]   |
| u'[t] == y[t] + Sin[t] |
| u[t] == t              |);

この場合，時間においては不連続になる．

NDSolveは"CollocationDirection"のオプションを自動的に設定する．

Wolfram Language code:

sol2 = First@NDSolve[{DAE2, x[0] == 1}, {x, y, u}, {t, 0, 10}, Method -> {"DAEInitialization" -> {"Collocation", "CollocationDirection" -> Automatic}}]

Wolfram Language code: GraphicsRow[Plot[Evaluate[#[t] /. sol2], {t, 0, 10}, PlotLabel -> #]& /@ {x, y, u}, ImageSize -> 500]

"Centered"の"CollocationDirection"を適用すると，解は発散する．

Wolfram Language code:

NDSolve[{DAE2, x[0] == 1}, {x, y, u}, {t, 0, 1}, Method -> {"DAEInitialization" -> {"Collocation", "CollocationDirection" -> "Centered"}}]

しかし，"Forward"か"Backward"の方向を使うと不連続の導関数を避けることで，正しい収束が得られる．

Wolfram Language code:

sol2 = First@NDSolve[{DAE2, x[0] == 1}, {x, y, u}, {t, 0, 10}, Method -> {"DAEInitialization" -> {"Collocation", "CollocationDirection" -> "Forward"}}]

ほとんどの場合，NDSolveは方向がどうあるべきかを自動的に決定する．しかし，特定の方向に不連続性がある可能性があることが分かっているならば，この設定を使ってそれを避けることができる．これがNDSolveの計算時間を削減する助けとなる．

選点

オプション"CollocationPoints"は，(5)で示されるように，級数近似(4)で使用される選点の数を指定する．選点によって，(4)の級数の次数が決まる．級数の次数はである．選点は選点範囲上に分布したチェビシェフ点である．

"CollocationPoints"の働きを見るために，次の指数3の例題を考える．

Wolfram Language code:

Clear[x1, x2, x3, vars, DAE3];
DAE3 = (⁠|                         |
| ----------------------- |
| x2'[t] + x1[t] == Sin[t] |
| x3'[t] + x2[t] == Sin[t] |
| x3[t] == Cos[t]          |⁠);
vars = {x1, x2, x3};

前の微分代数方程式の解は，解析的に求められる．矛盾のない厳密な初期条件は以下の通りである．

Wolfram Language code:

TableForm[{{-2, 1}, {0, 2}, {1, 0}}, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

メソッド"CollocationPoints"のデフォルト設定で正しい初期化が見付かる：

Wolfram Language code: sol3 = First@NDSolve[DAE3, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> "Collocation" }];

Wolfram Language code:

TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol3, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

選点は手動で設定することができる．"CollocationPoints"->2という設定は，(4)の設定と等価である．つまり，解は定数と線形の基底関数を使って近似されるということである．

定数と線形の基底関数を使って初期化を実行する：

Wolfram Language code:

sol3 = First@NDSolve[DAE3, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation" , "CollocationPoints" -> 2}}];

Wolfram Language code:

TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol3, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

たとえ反復が収束しても，結果は矛盾のない正しい初期条件ではない．この結果は指定された方程式を十分満足しているので，メッセージは生成されない．

指定された方程式の初期条件をテストする：

Wolfram Language code: Apply[Subtract, Flatten[DAE3], {1}] /. sol3 /. t -> 0

高指数の方程式では矛盾のない初期条件は微分された系にも使えることを思い出してほしい．

選点の数を増やし高階の基底関数を許可することで，アルゴリズムは想定された結果に収束することができる．

選点の範囲のある初期化問題を解くよう試みる：

Wolfram Language code:

initialConditions = Table[
	sol3 = First@NDSolve[DAE3, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation" , "CollocationPoints" -> m}}];
	{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0} /. sol3, {m, 2, 10}
	];

結果を選点の関数としてプロットすると，解がどのように収束するかが分かる．Red，Green，Blueはそれぞれ変数 x1，x2，x3に対応する：

Wolfram Language code:

valuePlot = ListLinePlot[Transpose[initialConditions[[All, 1]]], PlotMarkers -> Automatic, PlotStyle -> {StandardRed, StandardGreen, StandardBlue}, Frame -> True, FrameLabel -> {Style["Collocation Points (m)", 18], Style["Values", 18]}, DataRange -> {2, 10}]; 
	
derivativePlot = ListLinePlot[Transpose[initialConditions[[All, 2]]], PlotMarkers -> Automatic, PlotStyle -> {StandardRed, StandardGreen, StandardBlue}, Frame -> True, FrameLabel -> {Style["Collocation Points (m)", 18], Style["Derivatives", 18]}, DataRange -> {2, 10}]; 
	
GraphicsRow[{valuePlot, derivativePlot}, ImageSize -> 500]

近似の次数は選点により決まるので，(4)で得られた結果の線形化された系は平方行列になり，線形方程式の結果の系はLinearSolveを使って解くことができる．しかし，結果の線形系が悪条件であったり特異であったりする可能性もある．そのような系は高感度といわれ，負担邸名反復となり得る．この問題を解決するためには，多くの場合LeastSquaresを使って方程式の過剰決定系を解くことが望ましい．

サブオプションの"ExtraCollocationPoints"を使うと，級数近似の次数は(4)のと同じままだが，(5)におけるの値は（は追加の選点）に修正される．それゆえ，これは(5)の方程式の過剰決定系になる．これは最小二乗法を使って解くことができる．

過剰決定系で初期化する：

Wolfram Language code:

First@NDSolve[DAE3, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation", "CollocationPoints" -> {Automatic, "ExtraCollocationPoints" -> 3}}}];
TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol3, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

選点の範囲

指数減少を行っていない非常に高指数の系の場合，選点法を実行する範囲の指定の他，ずっと高い選点の次数を使うことが必要かもしれない．デフォルトでは，選点の範囲は作業精度と，展開で使用される基底関数に基づいて計算される．

例として，厳密解のあるい指数6の系を使う．

代数微分方程式の指数6の線形系に対する変数と方程式を定義する：

Wolfram Language code:

Clear[x1, x2, x3, x4, x5, x6];
DAE4 = (⁠|                         |
| ----------------------- |
| x6'[t] + x5[t] == Sin[t] |
| x5'[t] + x4[t] == Sin[t] |
| x4'[t] + x3[t] == Sin[t] |
| x3'[t] + x2[t] == Sin[t] |
| x2'[t] + x1[t] == Sin[t] |
| x6[t] == Cos[t]          |⁠);
vars = {x1, x2, x3, x4, x5, x6};

この場合の厳密解はDSolveを使って求めることができる：

Wolfram Language code: exactSol = First@DSolve[DAE4, Through[vars[t]], t]

初期条件の表を作る：

Wolfram Language code: TableForm[Transpose[{exactSol, D[exactSol, t]}] /. t -> 0]

初期化にデフォルトのオプションを使うことで，項のいくつかに誤った結果を得る：

Wolfram Language code:

sol3 = First@NDSolve[DAE4, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation"}}];
	
TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol3, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

通常，系の指数より高次の近似を使う必要がある．これは指数6の系なので，(4)の高次近似に依存しなければならない．前のサブセクションで言及した通り，これはオプション"CollocationPoints"を使って実行することができる．

次数15の近似に対して十分な点を使って初期化する：

Wolfram Language code:

sol3 = First@NDSolve[DAE4, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation", "CollocationPoints" -> 15}}];
	
TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol3, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

解は見付かるが，矛盾のない厳密な初期条件からの重大なエラーがある．この大きなエラーの原因は，この例題ではデフォルトオプションが近すぎるために目立った数値丸め誤差があるのである．

数値丸め誤差を避けるために，長さ1の拡張された選点範囲上に広がる15個の選点を使って初期化する：

Wolfram Language code:

sol3 = First@NDSolve[DAE4, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation", "CollocationPoints" -> 15, "CollocationRange" -> 1}}];
	
TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol3, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

"CollocationPoints"と"CollocationRange"の間には密接な関係がある．選点範囲は，過剰な丸め誤差を避け，高精度の計算を可能にするために選点の数に基づいて計算される．明示的に指定されていない限り， "CollocationRange"はとして計算される．ここではWorkingPrecisionオプションの設定であり，は"CollocationPoints"メソッドオプションの設定である．

初期条件の指定

常微分方程式と微分代数方程式の主な違いの一つに，常微分方程式系を積分するためには，すべての変数に対する初期条件を与えなければならないというものがある．その一方，微分代数方程式では，これらの初期条件のいくつかは固定されており，系により直接計算されなければならない．詳しく見るために，以下の常微分方程式を考える．

上記の系に対して，初期条件は任意に指定することができる．これで同じ系の修正を以下のように考えることができる．

この系は微分代数方程式であるが，力学は常微分方程式系と同じである．大きな違いは，初期条件とは固定されていないが代数方程式によって決定されるということである．

上の例はどの変数を固定するかについての決定ができるほど簡単で，初期条件は手動で計算される．しかし，慶賀もっと複雑で大きくなると，別のツールに依存する必要がある．

ある場合，初期条件は方程式系によって完全に決定される．この場合を例示するために，次の線形代数微分方程式系を考える．

初期条件の指定されていない系を初期化する：

Wolfram Language code:

eqns = {Subscript[x, 2]'[t] + Subscript[x, 1][t] == Sin[t], Subscript[x, 3]'[t] + Subscript[x, 2][t] == Sin[t], Subscript[x, 3][t] == Cos[t]} ;
vars = {Subscript[x, 1], Subscript[x, 2], Subscript[x, 3]};

Wolfram Language code: sol = First@NDSolve[eqns, vars, {t, 0}]

矛盾のない初期条件を得る：

Wolfram Language code:

vt = Through[vars[t]];
xic = Thread[vt -> (vt /. sol) ] /. t -> 0

Wolfram Language code: Thread[D[vt, t] -> ((D[vt, t] /. t -> 0) /. sol) ] /. t -> 0

通常，どの条件が固定されておりどの条件が自由であるかを知ることは不可能であるということに注意しなければならない．正しい初期条件および初期条件の桁数を察する方法の一つとして，初期条件なしで微分代数方程式を直接NDSolveに入れ，メソッド"DAEInitialization"のオプションを操作するというものがある．

一階の系として表される振り子の系(2)を考える．系の状態を一意に定義するためには，2つの初期条件だけが必要である．

Wolfram Language code:

Clear[x, xp, y, yp, λ, vars, DAEPendulum];
DAEPendulum = (⁠|                        |
| ---------------------- |
| x'[t] == xp[t]          |
| xp'[t] == λ[t] x[t]     |
| y'[t] == yp[t]          |
| yp'[t] == λ[t] y[t] - 1 |
| x[t]^2 + y[t]^2 == 1    |⁠) ; 
vars = {x, xp, y, yp, λ};

まず，この微分代数方程式系に対して与えられる初期条件がない場合を考える．デフォルトでは，最初の推定は指定されていないすべての変数に対して1が取られる．指数減少が実行されていないので，NDSolveは自動的に選点法を選んだ．

Wolfram Language code: sol1 = First@NDSolve[DAEPendulum, vars, {t, 0, 0}]

Wolfram Language code:

TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol1, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

アルゴリズムが無限数の可能な集合から可能な初期条件集合を一つ見付けることができたことが分かる．この時点で，初期の推定に非常に高感度である．また，反復の収束は保障されないことも言っておかなければならない．

最初の推定は選点方のオプション"DefaultStartingValue"を使って変更することができる．初期の推定-1が使われると，有効な解の異なる集合が見付かる．

Wolfram Language code:

sol1 = First@NDSolve[DAEPendulum, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation", "DefaultStartingValue" -> -1}}]; 
	
TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol1, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

先に述べたように，2つの初期条件は系の状態を固定するのに十分である．振り子の場合，値とで系が固定されている．

2つの状態を固定して，振り子の系の初期条件を計算する：

Wolfram Language code:

sol1 = First@NDSolve[{DAEPendulum, x[0] == 1, yp[0] == 1}, vars, {t, 0, 0}]; 
	
TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol1, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

状態が本当に固定されていると，理論的には系は最初の推定に対して感度が低くなる．

異なる初期値を持つ初期条件を計算する：

Wolfram Language code:

Table[
	sol1 = First@NDSolve[{DAEPendulum, x[0] == 1, yp[0] == 1}, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation", "DefaultStartingValue" -> dsv}}]; 
	
	TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol1, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]
	, {dsv, {1, -1}}]

ある系では，初期条件の計算は初期推定に敏感である．それを示すために以下の系を考える．

Wolfram Language code:

DAE = (|                                     |
| ----------------------------------- |
| Derivative[1][x1][t] == x3[t]        |
| x2[t] (1 - x2[t]) == 0               |
| x1[t] x2[t] + x3[t] (1 - x2[t]) == t |);

x₂[t]0のこの問題では，一般解は{x₁[t]x₁[0]+t²/2,x₂[t]0,x₃[t]t}であり，解を決定するために x₁[0]を与えることが必要である．初期推定値としてデフォルト値の1を使うことで，アルゴリズムは指定された初期条件を使うよう強制される．

Wolfram Language code:

vars = {Subscript[x, 1], Subscript[x, 2], Subscript[x, 3]};
sol3 = NDSolve[{DAE, Subscript[x, 1][0] == 1}, vars, {t, 0, 0}, Method -> {"DAEInitialization" -> "Collocation"}]

このコードは異なる解の分岐に切り替わる．

Wolfram Language code: Through[vars[0]] /. sol3

これは，状態に対するデフォルトの初期推定が1であるために起こる．ニュートン反復の間に，2つ目の方程式は完全に満足され，その結果反復は発散するので，系は特異になる．完全に失敗するのを避けるため，このアルゴリズムはの指定された条件を修正するよう強制される．

この動作を避けるためによりよい初期推定値を与えることができる．初期推定値に-1を使うと，反復は想定された矛盾のない初期状態に収束する．

Wolfram Language code:

sol3 = NDSolve[{DAE, Subscript[x, 1][0] == 1}, vars, {t, 0, 1}, Method -> {"DAEInitialization" -> {"Collocation", "DefaultStartingValue" -> -1}}] ; 
Through[vars[0]] /. sol3

Wolfram Language code:

Plot[Evaluate[{Subscript[x, 1][t], Subscript[x, 3][t]} /. First[sol3]], {t, 0, 1}, 
	PlotStyle -> {StandardRed, StandardBlue}]

QR法

QR分解はShampineの研究[S02]に基づいている．このアルゴリズムは指数1および指数0の系を意図したものである．このアルゴリズムは効率的なQR分解法を利用して系の変数とその導関数を分離する．そのようにすることにおいて，各反復ステップで1つの系の解が2つの小さな方程式系に分解される．個の変数を持つ微分代数方程式系では，選点法が使われたら大きさの行列であるのに対し，QR法は最大で大きさの行列を扱う必要がある．ほとんどの場合，このアルゴリズムはよりずっと小さい行列を扱う必要しかない．このことで，このメソッドは非常に大きい系を扱うのに非常に効率的なものとなっている．このセクションでは，もとになっているアルゴリズムを簡単に説明する．

形式の指数が減少した系がある場合，まずそれを初期推定値およびについて

として線形化する．QR分解はとなるようにヤコビ行列に施される．ここでは直交行列，は上三角行列，は置換行列を表す．微分代数方程式の場合，行列は完全行階数を持たない．行列の階数により，線形化された系は以下のように書くことができる．

変換された系は，2段階で解かれる．

一旦要素が見付かると，変換された変数は置換行列を使ってもとの形式に戻される．このプロセスは収束まで繰り返される．

QR法は，メソッド"DAEInitialization"のMethod->{"DAEInitialization"->{"QR",qropts}}（qropts はQR法のオプション）を使って呼び出すことができる．

"DefaultStartingValue"	Automatic	非線形反復で指定されていない成分の初期値．可能な値は1要素の値 ig/{ig}，2数のリスト{ig1,ig2}である．ここで ig1,ig2はそれぞれ変数とその導関数の初期値である
"MaxIterations"	100	最大反復回数

QR初期化のオプション

振り子の例を考える．

Wolfram Language code:

Clear[x, y, λ, t, pendulumDAE];
pendulumDAE = (⁠|                        |
| ---------------------- |
| x''[t] == x[t] λ[t]     |
| y''[t] == y[t] λ[t] - 1 |
| x[t]^2 + y[t]^2 == 1    |⁠) ;

これが高指数の微分代数方程式だとすると，系に対して指数減少法を施さなければならない．系の指数が1か0に減少されたら，変数とその導関数すべてが初期化されなければならない．

初期化に自動の指数減少法とQRアルゴリズムを使って，振り子の微分代数方程式を解く：

Wolfram Language code:

sol = NDSolve[{pendulumDAE, x[0] == 1, y'[0] == 1}, {x, y, λ}, {t, 0, 10}, Method -> {"IndexReduction" -> Automatic, "DAEInitialization" -> "QR"}]

Wolfram Language code: GraphicsRow[Plot[Evaluate[#[t] /. sol], {t, 0, 10}, PlotLabel -> #]& /@ {x, y, λ}]

上記の例題では，矛盾のない初期値の集合となるように初期条件が与えられている．初期条件が矛盾する場合，アルゴリズムは矛盾のない初期値を見付けようと試みる．

矛盾のある初期条件を含むQRアルゴリズムを使って，振り子の微分代数方程式を解く：

Wolfram Language code:

solQR1 = NDSolve[{pendulumDAE, x[0] == 2, y'[0] == 1}, {x, y, λ}, {t, 0, 10}, Method -> {"IndexReduction" -> Automatic, "DAEInitialization" -> "QR"}]

Wolfram Language code: GraphicsRow[Plot[Evaluate[#[t] /. solQR1], {t, 0, 10}, PlotLabel -> #]& /@ {x, y, λ}, ImageSize -> 500]

修正された初期条件を調べる：

Wolfram Language code: {x[0], y[0], λ[0]} /. solQR1

Wolfram Language code: {x'[0], y'[0], λ'[0]} /. solQR1

オプション"DefaultStartingValue"を指定することで，反復に使われる初期推定値を変更することができる．矛盾のある初期条件の場合，これにより矛盾のない異なる集合が導かれる．

QRとオプション"DefaultStartingValue"を使って，振り子の微分代数方程式を解く：

Wolfram Language code:

solQR2 = NDSolve[{pendulumDAE, x[0] == 2, y'[0] == 1}, {x, y, λ}, {t, 0, 1}, Method -> {"IndexReduction" -> Automatic, "DAEInitialization" -> {"QR", "DefaultStartingValue" -> -1}}];

Wolfram Language code: GraphicsRow[Plot[Evaluate[#[t] /. solQR1], {t, 0, 10}, PlotLabel -> #]& /@ {x, y, λ}, ImageSize -> 500]

異なる初期推定を使って得られた修正初期条件を調べる：

Wolfram Language code: {x[0], y[0], λ[0]} /. solQR2

Wolfram Language code: {x'[0], y'[0], λ'[0]} /. solQR2

BLT法

BLT法は，指数減少した系をブロック下三角形式に変換して，その系の部分集合を反復的に解いていくことを基本としている．部分集合自体はガウス・ニュートン法で解かれる．指数減少した系を生み出す指数減少アルゴリズムは，すでにブロック下三角形式を持つ傾向にあり，より小さい部分系をそれに対応する小さい行列計算で操作することにより計算の複雑性を減少させるため，この方法は便利なことがある．また，このアルゴリズムは初期条件の矛盾を検証して修正を試みる．

BLT法はメソッド"DAEInitialization"でMethod->{"DAEInitialization"->{"BLT",bltopts}}を使うことで呼び出すことができる．ここで bltopts はBLT法のオプションである．

"DefaultStartingValue"	Automatic	非線形の反復において指定されていない成分の初期値．可能な値は1要素の値 ig/{ig}か2数のリスト{ig1,ig2}である．ここで ig1,ig2 はそれぞれ変数とその導関数の初期値である．
"MaxIterations"	100	最大反復回数

BLT法のオプション

BLT法を使う利点を見るためには，指数減少アルゴリズムが，指数1か0である拡張された方程式系に至るということを認識することが大切である．発熱反応を記述する以下の化学反応系を考えてみる．

Wolfram Language code:

eqns = (⁠|                                              |
| -------------------------------------------- |
| c'[t] + c[t] + r[t] == 4 + t + t ^ 2          |
| T'[t] + 2 T[t] + r[t] + Tc[t]  == 1 + Exp[-t] |
| 1 + T[t]Log[r[t] / c[t]] == 0                 |
| c[t] == Cosh[t - 1]                           |⁠) ; 
vars = {c, r, T, Tc};

Wolfram Language code:

state = First@NDSolve`ProcessEquations[eqns, vars, t, Method -> {"IndexReduction" -> True}];
residual = state["NumericalFunction"]["FunctionExpression"];

下の微分代数方程式は4つの方程式から成っているが，指数減少された系は8つの方程式からできている．指数減少された方程式は以下の通りである．

Wolfram Language code: MatrixForm[residual[[2]]]

拡張された系には，ダミーの導関数が含まれたことに注目する．この系をさらに分析すると，系の構造を見ることができる．

指数減少された系に対する接続行列を生成する：

Wolfram Language code:

neq = residual[[2]];
nvars = residual[[1, 2 ;; 9]];
incidenceMatrix = NDSolve`StructuralIncidenceArray[neq, nvars];
MatrixForm[incidenceMatrix]

その系のブロック三角順序を見付ける：

Wolfram Language code:

{bltrow, bltcol} = SparseArray`BlockTriangularOrdering[incidenceMatrix];
TableForm[Table[{nvars[[bltcol[[indx]]]], neq[[bltrow[[indx]]]]}, {indx, 1, Length[bltrow]}], TableHeadings -> {None, {"solve for", "using Equation"}}, TableAlignments -> "Center"]

上記は，拡張された系が実際には逐次的に解くことができるということを示している．それゆえ，9つの方程式から成る1つの系を持つのではなく，方程式を1つずつ逐次解いていくことができるのである．初期化の見地から，t->0と設定することで変数の初期化が可能になる．

Wolfram Language code:

{bltrow, bltcol} = SparseArray`BlockTriangularOrdering[incidenceMatrix];
TableForm[Table[{nvars[[bltcol[[indx]]]], neq[[bltrow[[indx]]]] /. t -> 0.}, {indx, 1, Length[bltrow]}], TableHeadings -> {None, {"solve for", "using Equation"}}, TableAlignments -> "Center"]

これがBLT法を使う初期化に潜む主要原理である．BLT法は，下三角形式を持つ大きい系を扱うときは特に便利で計算的に効率がよいとされている．

初期化にBLTを使って化学反応系を解く：

Wolfram Language code: sol = NDSolve[eqns, vars, {t, 0, 3.5}, Method -> {"IndexReduction" -> True, "DAEInitialization" -> "BLT"}] ;

Wolfram Language code: Row[Plot[Evaluate[#[t] /. sol], {t, 0, 3.5}, PlotLabel -> #, ImageSize -> Small]& /@ vars]

初期条件が矛盾している場合，BLT法は初期条件を適切に修正することができる．しかしながら，BLT法で得られる結果とQR法で得られる結果は異なる可能性があることに注意する必要がある．初期条件が矛盾した振り子の例を考えてみる．

Wolfram Language code:

Clear[x, y, λ, t];
PendulumDAE = (⁠|                        |
| ---------------------- |
| x''[t] == x[t] λ[t]     |
| y''[t] == y[t] λ[t] - 1 |
| x[t]^2 + y[t]^2 == 1    |⁠) ;

QR法と矛盾した初期条件を使って，振り子の微分代数方程式を解く：

Wolfram Language code:

solQR = NDSolveValue[{PendulumDAE, x[0] == 2, y'[0] == 1}, {x, y, λ}, {t, 0, 1}, Method -> {"IndexReduction" -> True, "DAEInitialization" -> "QR"}];

Wolfram Language code: {Through[solQR[0]], D[Through[solQR[t]], t] /. t -> 0}

BLTアルゴリズムと矛盾した初期条件を使って，振り子の代数微分方程式を解く：

Wolfram Language code:

solBLT = NDSolveValue[{PendulumDAE, x[0] == 2, y'[0] == 1}, {x, y, λ}, {t, 0, 1}, Method -> {"IndexReduction" -> True, "DAEInitialization" -> "BLT"}];

Wolfram Language code: {Through[solBLT[0]], D[Through[solBLT[t]], t] /. t -> 0}

矛盾した初期条件

微分代数方程式系の定式化とその解において最も難しいことに一つに，最初の初期条件が矛盾のないものでなければならないということがある．これは単純に見えるが，高指数の微分代数方程式系では，初期条件も満足しなければならない隠れた条件／方程式に簡単に遭遇する．NDSolveが矛盾した初期条件に遭遇すると，それを修正しようと試みる．

隠れた制約条件と矛盾のある初期条件についての要点を説明するために，振り子の系(2)を考えよう．

最初の場合として，初期条件が系の代数的制約条件に反するときを考える．これはの初期条件が制約条件に常に反するように意図的に設定することで行える．

矛盾した初期条件を定義する：

Wolfram Language code: ics = {x[0] == 2, yp[0] == 1};

系がNDSolveに渡されると，NDSolve::ivresというメッセージが生成され初期条件が矛盾していることを表す．その後NDSolveは矛盾した値をその反復の初期推定として使って，矛盾のない初期条件を見付けようとする．

矛盾した値を初期推定として使って，矛盾のない初期条件の集合を見付ける：

Wolfram Language code:

sol3 = First@NDSolve[{x'[t] == xp[t], xp'[t] == λ[t] x[t], y'[t] == yp[t], yp'[t] == λ[t] y[t] - 1, x[t] ^ 2 + y[t] ^ 2 == 1, ics}, {x, xp, y, yp, λ}, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation"}}]; 
	
vars = {x, xp, y, yp, λ};
TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol3, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

指定された条件が矛盾しているかどうかがはっきりしない場合がある．初期条件は制約条件を満足してるように見えても，隠れた制約条件は満足していないかもしれない．

振り子の系の初期条件を指定する：

Wolfram Language code: ics = {x[0] == 1, y[0] == 0, yp[0] == 1, xp[0] == 0.1};

これらの初期条件は方程式を満足するが，暗示的な制約条件を満足することはできない．この制約条件は隠れており指数減少によってのみ見付かるので，明白ではない．選点法は指数減少を実行しないで，もとの方程式に作用する．指数減少した系がなくても，選点法は初期条件の矛盾を検出することができ，それを修正しようと試みる．

矛盾のない初期条件集合を得る：

Wolfram Language code:

sol3 = First@NDSolve[{x'[t] == xp[t], xp'[t] == λ[t] x[t], y'[t] == yp[t], yp'[t] == λ[t] y[t] - 1, x[t] ^ 2 + y[t] ^ 2 == 1, ics}, {x, xp, y, yp, λ}, {t, 0, 0}, Method -> {"DAEInitialization" -> {"Collocation"}}]; 
	
vars = {x, xp, y, yp, λ};
TableForm[Transpose[{Through[vars[0]], D[Through[vars[t]], t] /. t -> 0}] /. sol3, TableHeadings -> {Through[vars[0]], {"Value", "Derivative"}}, TableAlignments -> "Center"]

イベントを使った微分代数方程式の再初期化

状態の変化がで起こる，離散変数を持った2つの方程式から成る系を考える．

系の変数と方程式を定義する：

Wolfram Language code:

Clear[x, y, s, t];
vars = {x[t], y[t], s[t]};
eqns = {x'[t] == y[t], y[t] == s[t] x[t], y[0] == 1, x[0] == 1};

状態変化を行うために，イベントのある系を解く：

Wolfram Language code:

sol = NDSolveValue[{eqns, WhenEvent[t == 1, s[t] -> 2]}, vars, {t, 0, 2}, DiscreteVariables -> s]  ;
Plot[sol, {t, 0, 2}, PlotStyle -> {StandardRed, StandardBlue, StandardGreen}]

この例ではとに対する初期条件が与えられている．初期化コードが矛盾のない初期条件を自動的に計算する．時間で，離散変数は1から2に変更しなければならないというイベントが発生する．において変数の再初期化が必要になる．変数は修正された変数と呼ばれる．再初期化中に，以下のようなストラテジーが適用される．

1. すべての状態変数，離散変数，修正された変数は固定しておく．状態変数の導関数が固定されたら（あるいは変更するよう指定されたら），その変数は変更することができる．

2. ストラテジー1が再初期化に失敗したら，状態変数と修正された変数を固定して，他の変数とその導関数は変更できるようにする．

3. ストラテジー2が失敗したら，修正された変数と，代数微分方程式の導関数を持つ状態変数を固定して，その他すべての変数とその導関数は変更できるようにする．

4. ストラテジー3が失敗したら，修正された変数だけを固定して，状態変数とその導関数はすべて変更できるようにする．

再初期化の直前，，，の値は以下の通りである．

におけるイベントの直前の値を得る：

Wolfram Language code: sol /. t -> 0.999999999

参照されたストラテジーをもとに，時間でストラテジー1と2を使うと，，，は変更してはならない．だとすると，これは明らかに代数方程式に反する．したがって，ストラテジー1と2は失敗に終わる．

ストラテジー3では，導関数を持つ変数は変更してはならない．よって，はとともに固定したままであるが，は変更してもよい．これによりの値はプロットで示されている（青い線）ように修正される．

におけるイベントの直後の再初期化された値を得る：

Wolfram Language code: sol /. t -> 1.000000001

WhenEventは，状態変化を含むイベント動作を同時評価することも逐次評価することもできる．

各状態変化の後，同じイベントでも，微分代数方程式は上のストラテジーを使って再初期化される．そのため，可能ならば，必要な状態変化に対して1つの規則だけを与えた方が効率できである．同じイベント時において実行される複数の状態変化があるときの最後の状態は，異なる可能性がある．なぜなら，各後続の評価と設定で新しい値が使われ，微分代数方程式ではこれらの値は再初期化に依存するかもしれないからである．追加のイベントを含む前の問題を考えてみる．

順に評価される状態変化に対して与えられる3つの規則を持つ問題を解く：

Wolfram Language code:

sol = NDSolveValue[{eqns, WhenEvent[t == 1, {s[t] -> 2, x[t] -> y[t], y[t] -> x[t]}]}, vars, {t, 0, 2}, DiscreteVariables -> s]  ; 
Plot[sol, {t, 0, 2}, PlotStyle -> {StandardRed, StandardBlue, StandardGreen}]

なぜこの解が上の解とは異なっているかを見るためには，一連の設定を追ってみる必要がある．まずが2に設定され，微分代数方程式は他の設定を評価する前に再初期化される．再初期化によっての値が2に変わる．そのようにして，上記のストラテジーが残差と矛盾のない解を見付けるからである．2つ目の動作の実行でが2（の現行値）に設定されるが，これは残差と矛盾する．したがってが代数的変数なので，その値は4()に設定される．最後にをに設定するとが2にリセットされ，新規に設定された代数的変数を修正することを避けるため，の値は1にリセットされる．

変化の結果を見る一つの方法は，設定と設定の間にPrintコマンドを入れることである．

設定と設定の間の変数値を出力する：

Wolfram Language code:

sol = NDSolveValue[{eqns, WhenEvent[t == 1, {Print["Before ", vars], s[t] -> 2, Print["After s[t]->2 ", vars], x[t] -> y[t], Print["After x[t]-> y[t] ", vars], y[t] -> x[t], Print["After y[t]-> x[t] ", vars]}]}, vars, {t, 0, 2}, DiscreteVariables -> s]  ; 
Plot[sol, {t, 0, 2}, PlotStyle -> {StandardRed, StandardBlue, StandardGreen}]

最後の2つの設定を入れ替えて解く：

Wolfram Language code:

sol = NDSolveValue[{eqns, WhenEvent[t == 1, {Print["Before ", vars], s[t] -> 2, Print["After s[t]->2 ", vars], y[t] -> x[t], Print["After y[t]-> x[t] ", vars], x[t] -> y[t], Print["After x[t]-> y[t] ", vars]}]}, vars, {t, 0, 2}, DiscreteVariables -> s]  ; 
Plot[sol, {t, 0, 2}, PlotStyle -> {StandardRed, StandardBlue, StandardGreen}]

設定の順序を替えると，異なった解になる．

規則が1つの場合初期化は1つだけなので，すべての変数が設定されると，与えられる値は矛盾がないものでなければならない．さもないと，再初期化ストラテジーは失敗するか，想定された値を与えないようになる．

単独の規則を使って，すべての変数を微分代数方程式に矛盾しない値に設定しようとする：

Wolfram Language code:

sol = NDSolveValue[{x'[t] == y[t], y[t] == s[t] x[t], y[0] == 1, x[0] == 1, WhenEvent[t == 1, {s[t], x[t], y[t]} -> {2, y[t], x[t]}]}, {x[t], y[t], s[t]}, {t, 0, 2}, DiscreteVariables -> s]

ここでは，，が固定されている．これは非常に厳しい制約条件であり，すべてのストラテジーが失敗する．

有効な単独規則内で変数を現行の値に設定すると，同じままであるよう強制することになる．

を2に変更し，

が

で同じままであるよう強制する：

Wolfram Language code:

sol = NDSolveValue[{x'[t] == y[t], y[t] == s[t] x[t], y[0] == 1, x[0] == 1, WhenEvent[t == 1, {s[t], y[t]} -> {2, y[t]}]}, {x[t], y[t], s[t]}, {t, 0, 2}, DiscreteVariables -> s]  ;
Plot[sol, {t, 0, 2}, PlotStyle -> {StandardRed, StandardBlue, StandardGreen}]

が同じままであるため，は変化するよう強制される．

微分代数方程式の例題

トーラスの問題

これはCampbell & Moore [CM95]からのテスト問題である．微分代数方程式は粒子の動きを3Dで表す．これは3つの位置変数，3つの速度変数，1つの係数を持つ指数3の微分代数方程式である．解多様体は四次元であり，厳密解は以下で与えられる．

解として使うパラメータと関数を定義する：

Wolfram Language code:

Clear[x, y, z, r, ρ, c1, x1, x2, x3, u1, u2, u3, λ, eqns, vars, exactTorusSol, torusSol, ics];
r = 10;
ρ = 5;
c1 = (1 - (r/Sqrt[x1[t]^2 + x2[t]^2]));

系の方程式と変数を定義する：

Wolfram Language code:

eqns = (⁠|                                                                           |
| ------------------------------------------------------------------------- |
| x1'[t] == u1[t]                                                            |
| x2'[t] == u2[t]                                                            |
| x3'[t] == u3[t]                                                            |
| u1'[t] == u3[t] Cos[t] - x3[t] Sin[t] - u2[t] + 2 c1 x1[t]  λ[t]           |
| u2'[t] == u3[t] Sin[t] + x3[t] Cos[t] + u1[t] + 2 c1 x2[t] λ[t]            |
| u3'[t] == -x3[t] + 2 x3[t] λ[t]                                            |
| x1[t]^2 + x2[t]^2 + x3[t]^2 - 2 r Sqrt[x1[t]^2 + x2[t]^2] + r^2 - ρ^2 == 0 |⁠);
vars = {x1, x2, x3, u1, u2, u3, λ};

厳密解は既知である：

Wolfram Language code:

exactTorusSol  =   {(ρ * Cos[2 Pi  - t] + r) * Cos[t], (ρ Cos[2 Pi  - t] + r) Sin[t], (ρ Sin[2 Pi - t])};
	
ics = Thread[Through[vars[[1 ;; -2]][0]] == Flatten[{exactTorusSol, D[exactTorusSol, t]}] /. t -> 0]

系を解き，可視化する：

Wolfram Language code:

torusSol = NDSolve[{eqns, ics}, vars, {t, 0, 10 Pi}, Method -> {"IndexReduction" -> {Automatic, "IndexGoal" -> 0}, "EquationSimplification" -> "Residual"}]

Wolfram Language code: ParametricPlot3D[Evaluate[{x1[t], x2[t], x3[t]} /. torusSol], {t, 0, 10 Pi}]

解析解と計算された解を比較する：

Wolfram Language code: errors = Transpose[{({x1[t], x2[t], x3[t]} /. First@torusSol), exactTorusSol}];

Wolfram Language code: Plot[Evaluate[#], {t, 0, 10 Pi}]& /@ errors

タンクからの圧力の解放

これは石油化学工学の動的プロセスを簡約したモデルである．このモデルの詳細は，Hairer et. al. [HLR89]をご参照いただきたい．これは指数2の系である．

系のパラメータを定義する：

Wolfram Language code:

Clear[m, M, μ, P, p, c, k, eqns, vars, sol];
α = 3.35; β = 0.075; γ = 0.0001; a = 49.58;b = 1.2;

方程式と変数を定義する：

Wolfram Language code:

eqns = (⁠|                                                                          |
| ------------------------------------------------------------------------ |
| (1/20) (-c[t] + k[t]) + Derivative[1][k][t] == 0                          |
| (1/75) (-99 + p[t]) + Derivative[1][c][t] + (Derivative[1][p][t]/15) == 0 |
| m[t] - μ[t] + Derivative[1][M][t] == 0                                    |
| -α - β m[t] - γ m[t]^2 + (p[t]/P[t]) == 0                                 |
| -a^2 + P[t]^2 - (μ[t]^2/b^2 k[t]^2) == 0                                  |
| M[t] - 20 P[t] == 0                                                       |
| -15 + 5 Tanh[10 - t] + μ[t] == 0                                          |⁠);
vars = {μ, P, M, m, c, k, p};

解を求めて可視化する：

Wolfram Language code:

sol = First@Quiet@NDSolve[{eqns, k[0] == 0.25, p[0] == 99.1}, vars, {t, 0, 40}, Method -> {"IndexReduction" -> {Automatic, "IndexGoal" -> 0}}]

Wolfram Language code: Row[Plot[Evaluate[#[t] /. sol], {t, 0, 40}, PlotRange -> All, PlotLabel -> #[t], ImageSize -> 125]& /@ vars]

非伸縮性スライダークランク

次の例題は制約条件付きの簡単な多体システムである．このシステムは固定面に固定された質量m1のクランク，質量m3のスライダー，質量m2の連結棒で構成されている．このスライダーは水平面でのみ移動が可能である．各部品の質量のため，このシステムの構築中に慣性モーメントを考慮する必要があった．問題の概略は以下に示す通りである．このシステムは2つの角度を使って完全に定義することができる．

Wolfram Language code: Clear[J1, J2, m1, m2, m3, L1, L2, LM, g, Fs, s1, LHS, RHS, eqns, contraints, vars, SliderCrankSol];

スライダクランクメカニズムの状態は2つの角度と原点からのスライダの距離で完全に表すことができる．

変数を定義する：

Wolfram Language code: vars = {α1, α2, z, λ1, λ2};

スライダー上で働く力を定義する：

Wolfram Language code: Fs[t] := -1;

運動の方程式は，力を分散し，ニュートンの法則と .を適用することによって導かれる．

微分系を定義する：

Wolfram Language code:

s1 = Sin[α1[t] - α2[t]];
c1 = Cos[α1[t] - α2[t]];
LM = (L1 L2 m2) / 2;
eqns = {(J1 + m2 L1^2) α1''[t] + LM c1 α2''[t] == -LM s1 α2'[t]^2
	 - (L1 Cos[α1[t]] λ1[t] + L1 Sin[α1[t]] λ2[t]), 
	LM c1 α1''[t] + J2 α2''[t] == LM s1 α1'[t]^2
	 - (L2 Cos[α2[t]] λ1[t] + L2 Sin[α2[t]] λ2[t]), 
	m3 z''[t] == -Fs[t] - λ2[t]};

この代数方程式は，系の幾何学を定義する：

Wolfram Language code:

constraints = (⁠|                                          |
| ---------------------------------------- |
| L1 Sin[α1[t]] + L2 Sin[α2[t]] == 0        |
| z[t] - L1 Cos[α1[t]] - L2 Cos[α2[t]] == 0 |⁠);

系の物理的パラメータを定義する．ここで

と

は慣性モーメントである：

Wolfram Language code: paramsFSC = {J1 -> 4.5, J2 -> 5.5, m1 -> 1, m2 -> 1, m3 -> 1, L1 -> 1, L2 -> 2, g -> 9.81};

系を解き，可視化する：

Wolfram Language code:

{sliderCrank} = NDSolve[{eqns, constraints, α1[0] == Pi / 4, α1'[0] == -1} /. paramsFSC, vars, {t, 0, 20}, Method -> {"IndexReduction" -> Automatic}];

Wolfram Language code: Row[Plot[Evaluate[#[t] /. sliderCrank], {t, 0, 20}, PlotRange -> All, PlotLabel -> #, ImageSize -> Small]& /@ vars]

スライダクランクをアニメーション化するための関数を定義する：

Wolfram Language code:

bar[loc_ : {0, 0}, rotate_ : 0, scale_ : 1, length_ : 1, width_ : .1, color_ : Hue[.6, .2, .9]] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	color, EdgeForm[Directive[LightDarkSwitched[Black, White], AbsoluteThickness[.5]]], 
	(* polygon instead of rounded rectangle in case texture is needed *)
	Polygon[ Flatten[{Table[loc + scale{width Sin[ a Pi], width Cos[a Pi]}, {a, 1, 2, 1 / 16}], 
	Table[{length, 0} + loc + scale{width Sin[ a Pi], width Cos[a Pi]}, {a, 0, 1, 1 / 16}]}, 1]], 
	AbsolutePointSize[5], LightDarkSwitched[Black, White], Point[loc], Point[loc + scale{length, 0}]
	}, rotate Degree, loc]
	}]
	]
	
barFlange[loc_ : {0, 0}, rotate_ : 0, scale_ : 1, length_ : 1, width_ : .1, size1_ : 1, size2_ : 1, color_ : Hue[.6, .2, .9]] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	color, EdgeForm[Directive[LightDarkSwitched[Black, White], AbsoluteThickness[.5]]], 
	Polygon[ {loc + scale{0, -width size1 }, loc + scale{0, width size1 }, loc + scale{length, width size2 }, loc + scale{length, -width size2 }}], 
	Lighter[color, .5], Disk[loc, scale width size1], Disk[loc + scale{length, 0}, scale width size2], 
	(* pins *)
	StandardGray, Disk[loc, .4scale width size1], Disk[loc + scale{length, 0}, .4scale width size2]
	}, rotate Degree, loc]
	}]
	]
	
dimensionLine[loc_ : {0, 0}, rotate_ : 0, scale_ : 1, text_ : "1", length_ : 1, offset_ : {0, 0}] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Line[{offset + loc + scale{0, 0}, offset + loc + scale{0, .1}}], 
	Arrowheads[{-Small, Small}], Arrow[{offset + loc + scale{0, .05}, offset + loc + scale{length, .05}}], 
	Line[{offset + loc + scale{length, 0}, offset + loc + scale{length, .1}}], 
	Text[text, offset + loc + scale{length / 2, .05}, Background -> LightDarkSwitched[White, Black]]
	}, rotate Degree, loc]
	}]
	]
	
fixedHinge[loc_ : {0, 0}, rotate_ : 0, scale_ : 1] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	GrayLevel[.85], EdgeForm[Directive[LightDarkSwitched[Black, White], AbsoluteThickness[.5]]], Polygon[ Flatten[{{loc + scale{-.25, -.5}}, Table[loc + scale{.15Sin[ a Pi], .15Cos[a Pi]}, {a, 1.55, 2.45, .9 / 16}], {loc + scale{.25, -.5}}}, 1]], 
	GrayLevel[.7], EdgeForm[], Rectangle[loc + scale{-.5, -.75}, loc + scale{.5, -.5}], 
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Line[{loc + scale{-.5, -.5}, loc + scale{.5, -.5}}], 
	AbsolutePointSize[5], LightDarkSwitched[Black, White], Point[{0, 0}]
	}, rotate Degree, loc]
	}]
	]
	
slider[loc_ : {0, 0}, rotate_ : 0, scale_ : 1, length_ : 1, width_ : .1, fixedLength_ : 2, pos_ : {0, 0}, color_ : Hue[0, .2, .9]] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	color, EdgeForm[{LightDarkSwitched[Black, White], AbsoluteThickness[.5]}], 
	Rectangle[loc + pos + scale{-length, -width}, loc + pos + scale{length, width}], 
	GrayLevel[.7], EdgeForm[], Rectangle[loc + scale{-fixedLength length, -width - .25}, loc + scale{fixedLength length, -width}], 
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Line[{loc + scale{-fixedLength length, -width}, loc + scale{fixedLength length, -width}}], 
	GrayLevel[.7], EdgeForm[], Rectangle[loc + scale{-fixedLength length, width + .25}, loc + scale{fixedLength length, width}], 
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Line[{loc + scale{-fixedLength length, width}, loc + scale{fixedLength length, width}}]
	}, rotate Degree, loc]
	}]
	]

システムの動きをアニメーション化する：

Wolfram Language code:

Animate[
	rotate = (α1[t] /. sliderCrank) * (180 / Pi);
	
	f1 = Graphics[{Red, Line[Map[Function[Evaluate[{#, λ1[#]}  /. sliderCrank]], Range[0, t, 0.025]]]}];
	f2 = Graphics[{Blue, Line[Map[Function[Evaluate[{#, λ2[#]}  /. sliderCrank]], Range[0, t, 0.025]]]}];
	
	{bLength, bfLength} = ({L1, L2} /. paramsFSC) * 0.3;
	Column[{
	Graphics[{
	(* parts *)
	bfAngle = ArcSin[bLength Cos[(90 - rotate) °] / bfLength] * 180 / Pi;
	sliderLoc = {bLength Cos[rotate Pi / 180] + bfLength Cos[bfAngle °], 0};
	slider[{.7, 0}, 0, .2, .7, .3, 3.5, sliderLoc - {.7, 0}, Hue[.05, .3, .9]], 
	bar[{0, 0}, rotate, 1, bLength, .04, Hue[.6, .2, .7]], 
	fixedHinge[{0, 0}, 0, .2], 
	bfLoc = {bLength Sin[(90 - rotate) °], bLength Cos[(90 - rotate) °]};
	barFlange[bfLoc, -bfAngle, 1, bfLength, .04, 2, 1.2, Hue[.6, .2, .9]], 
	(* axis *)
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Arrowheads[Medium], 
	Arrow[{{-.25, 0}, {.5, 0}}], (* x *)
	Arrow[{{0, -.05}, {0, .5}}](* y *)
	(* dimensions *)
	, 
	(* force arrow *)
	 Opacity[1], Hue[0, 1, .8], AbsoluteThickness[2], Arrowheads[Medium], Arrow[{sliderLoc + {.35, 0}, sliderLoc + {.15, 0}}], 
	LightDarkSwitched[Black, White], Text[Style["F(t)", 12, Italic], sliderLoc + {.37, 0}, {-1, 0}]
	}, PlotRange -> {{-.5, 1.5}, {-.5, .5}}, ImageSize -> 400]
	, 
	Show[f1, f2, Frame -> True, PlotRange -> {{0, 10}, {-2, 4}}, PlotLabel -> Row[{Style["Tension in slider (λ1)", 14, StandardRed], "        ", Style["Tension in rod (λ2)", 14, Blue]}], ImageSize -> 300]}, Alignment -> Center]
	, {t, 0, 10}, SaveDefinitions -> True, AnimationRunning -> False]

接続軸にバネ制振装置を装備したスライダクランク

これは簡単なスライダクランク機構の拡張である．クランクと接続軸の間にバネと制振装置を装備する．このばねはクランクの動きの範囲を制約し，制振装置はゆっくりとシステムを安定した状態にする．系の概略は以下の通りである．

Wolfram Language code:

Clear[c, L0, Lp, a, J1, J2, m1, m2, k, L, Lp, f, F, LHS, RHS1, RHS2, y1, y2, y3, λ1, λ2, eqns, constraints, vars, SliderCrankSol2, sol];

スライダクランクメカニズムの状態は2つの角度と原点からのスライダの距離で完全に記述することができる．

変数を定義する：

Wolfram Language code: vars = {y1, y2, y3, λ1, λ2};

スライダ上で働く力を定義する：

Wolfram Language code: F[t_] := 2;

運動の方程式は，力を分解し，ニュートンの法則とを適用することによって導かれる．

微分系を定義する：

Wolfram Language code:

L = Sqrt[y2[t]^2 - y2[t](Cos[y3[t]] + a Cos[y1[t]]) + ((1 + a^2)/4) + a (Cos[y3[t] - y1[t]]/2)] ;
Lp  = (1/2 L)(2 y2[t] y2'[t] - y2'[t] (Cos[y3[t]] + a Cos[y1[t]]) + y2[t] (Sin[y3[t]] y3'[t] + a Sin[y1[t]] y1'[t]) - a Sin[y3[t] - y1[t]] ((y3'[t] - y1'[t])/2));
f = k * (L - L0) + c Lp;
	
eqns = {J1 y1''[t] == -(a f (Sin[(1/2) (-y1[t] + y3[t])] + Sin[y1[t]] y2[t])/2 L) - (1/2) Sin[y1[t]] λ1[t] - (1/2) Cos[y1[t]] λ2[t], m2 y2''[t] == F[t] + (f ((1/2) a Cos[y1[t]] + (1/2) Cos[y3[t]] - y2[t])/L) - λ1[t], J2 y3''[t] == -F[t] Sin[y3[t]] - (f ((1/2) a Sin[y1[t] - y3[t]] + Sin[y3[t]] y2[t])/2 L) - (1/2) Sin[y3[t]] λ1[t] - Cos[y3[t]] λ2[t]};

代数的制約条件を定義する：

Wolfram Language code:

constraints = (⁠|                                          |
| ---------------------------------------- |
|  y2[t] - a Cos[y1[t]] - a Cos[y3[t]] == 0 |
| a Sin[y1[t]]  +   Sin[y3[t]] == 0         |⁠);

系のパラメータを定義する：

Wolfram Language code: paramsFlexSC = {a -> 1 / 2, J1 -> 1, J2 -> 2, m1 -> 1, m2 -> 1, L0 -> 1, k -> 1, c -> 4};

系を解き，可視化する：

Wolfram Language code:

SliderCrankSol2 = First@NDSolve[{eqns, constraints, y1[0] == Pi / 4, y1'[0] == -1} /. paramsFlexSC, vars, {t, 0, 40}, Method -> "IndexReduction" -> Automatic];
Row[Plot[Evaluate[#[t] /. SliderCrankSol2], {t, 0, 40}, PlotRange -> All, PlotLabel -> #, ImageSize -> Small]& /@ vars]

アニメーションを生成するための関数を定義する：

Wolfram Language code:

bar[loc_ : {0, 0}, rotate_ : 0, scale_ : 1, length_ : 1, width_ : .1, color_ : Hue[.6, .2, .9]] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	color, EdgeForm[Directive[LightDarkSwitched[Black, White], AbsoluteThickness[.5]]], 
	(* polygon instead of rounded rectangle in case texture is needed *)
	Polygon[ Flatten[{Table[loc + scale{width Sin[ a Pi], width Cos[a Pi]}, {a, 1, 2, 1 / 16}], 
	Table[{length, 0} + loc + scale{width Sin[ a Pi], width Cos[a Pi]}, {a, 0, 1, 1 / 16}]}, 1]], 
	AbsolutePointSize[5], LightDarkSwitched[Black, White], Point[loc], Point[loc + scale{length, 0}]
	}, rotate Degree, loc]
	}]
	]
	
barFlange[loc_ : {0, 0}, rotate_ : 0, scale_ : 1, length_ : 1, width_ : .1, size1_ : 1, size2_ : 1, color_ : Hue[.6, .2, .9]] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	color, EdgeForm[Directive[LightDarkSwitched[Black, White], AbsoluteThickness[.5]]], 
	Polygon[ {loc + scale{0, -width size1 }, loc + scale{0, width size1 }, loc + scale{length, width size2 }, loc + scale{length, -width size2 }}], 
	Lighter[color, .5], Disk[loc, scale width size1], Disk[loc + scale{length, 0}, scale width size2], 
	(* pins *)
	Gray, Disk[loc, .4scale width size1], Disk[loc + scale{length, 0}, .4scale width size2]
	}, rotate Degree, loc]
	}]
	]
	
barMidFlange[loc_ : {0, 0}, rotate_ : 0, scale_ : 1, length_ : 2, width_ : .1, size1_ : 1, size2_ : 1, mid_ : 1, color_ : Hue[.6, .2, .9]] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	color, EdgeForm[Directive[LightDarkSwitched[Black, White], AbsoluteThickness[.5]]], 
	Polygon[ Flatten[{{loc + scale{0, -width size1 }, loc + scale{0, width size1 }}, Table[{length, 0} + loc + scale{width Sin[ a Pi], width Cos[a Pi]}, {a, 0, 1, 1 / 16}]}, 1]], 
	Lighter[color, .5], Disk[loc, scale width size1], Disk[loc + scale{mid length, 0}, scale width size2], 
	(* pins *)
	Gray, Disk[loc, .4scale width size1], Disk[loc + scale{mid length, 0}, .4scale width size2]
	}, rotate Degree, loc]
	}]
	]
	
coilSpring[loc_ : {0, 0}, force_ : 0, segments_ : 4, height_ : .1, endLength_ : .1, thickness_ : .015, ghost_ : None] := 
	Module[{f, s, e}, 
	{
	f = force + 1;
	s = (f - 2endLength) / (segments + .5);
	e = endLength + s / 4;
	GraphicsGroup[{
	LightDarkSwitched[Black, White], Flatten@{
	Thickness[thickness], JoinForm["Round"], CapForm["Butt"], 
	If[ghost === "Ghost", GrayLevel[0, .2], LightDarkSwitched[Black, White]], Line[Flatten[{{loc + {0, 0}, loc + { endLength, 0}}, 
	Flatten[Table[{loc + { e + a s, -height}, loc + {e + a s + .5s, height}}, {a, 0, segments - 1}], 1], 
	{loc + {f - e, -height}, loc + {f - endLength, 0}, loc + {f - endLength, 0}, loc + {f, 0}}}, 1]], 
	If[ghost =!= "Ghost", {White, CapForm["Round"], Thickness[.5thickness], Line[{loc + {0, 0}, loc + { endLength, 0}}], 
	Table[Line[{loc + { e + a s, -height}, loc + {e + a s + .5s, height}}], {a, 0, segments - 1}], 
	Line[{loc + {f - e, -height}, loc + {f - endLength, 0}}], Line[{loc + {f - endLength, 0}, loc + {f, 0}}]}]}
	}]
	}]
	
dimensionLine[loc_ : {0, 0}, rotate_ : 0, scale_ : 1, text_ : "1", length_ : 1, offset_ : {0, 0}] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Line[{offset + loc + scale{0, 0}, offset + loc + scale{0, .1}}], 
	Arrowheads[{-Small, Small}], Arrow[{offset + loc + scale{0, .05}, offset + loc + scale{length, .05}}], 
	Line[{offset + loc + scale{length, 0}, offset + loc + scale{length, .1}}], 
	Text[text, offset + loc + scale{length / 2, .05}, Background -> White]
	}, rotate Degree, loc]
	}]
	]
	
fixedHinge[loc_ : {0, 0}, rotate_ : 0, scale_ : 1] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	GrayLevel[.85], EdgeForm[Directive[LightDarkSwitched[Black, White], AbsoluteThickness[.5]]], Polygon[ Flatten[{{loc + scale{-.25, -.5}}, Table[loc + scale{.15Sin[ a Pi], .15Cos[a Pi]}, {a, 1.55, 2.45, .9 / 16}], {loc + scale{.25, -.5}}}, 1]], 
	GrayLevel[.7], EdgeForm[], Rectangle[loc + scale{-.5, -.75}, loc + scale{.5, -.5}], 
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Line[{loc + scale{-.5, -.5}, loc + scale{.5, -.5}}], 
	AbsolutePointSize[5], LightDarkSwitched[Black, White], Point[{0, 0}]
	}, rotate Degree, loc]
	}]
	]
	
slider[loc_ : {0, 0}, rotate_ : 0, scale_ : 1, length_ : 1, width_ : .1, fixedLength_ : 2, pos_ : {0, 0}, color_ : Hue[0, .2, .9]] := 
	Module[{}, 
	GraphicsGroup[{
	Rotate[{
	color, EdgeForm[{LightDarkSwitched[Black, White], AbsoluteThickness[.5]}], 
	Rectangle[loc + pos + scale{-length, -width}, loc + pos + scale{length, width}], 
	GrayLevel[.7], EdgeForm[], Rectangle[loc + scale{-fixedLength length, -width - .25}, loc + scale{fixedLength length, -width}], 
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Line[{loc + scale{-fixedLength length, -width}, loc + scale{fixedLength length, -width}}], 
	GrayLevel[.7], EdgeForm[], Rectangle[loc + scale{-fixedLength length, width + .25}, loc + scale{fixedLength length, width}], 
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Line[{loc + scale{-fixedLength length, width}, loc + scale{fixedLength length, width}}]
	}, rotate Degree, loc]
	}]
	]
	
springDamper[loc2_ : {0, 0}, scale_ : 1, force2_ : 0, height2_ : .1, endLength2_ : 1, thickness_ : 2] := 
	Module[{f, s, e, loc = {endLength2, .5}, force = 1, segments = 4, height = .4, endLength = .5}, 
	{
	loc = scale loc + loc2;
	f = Dynamic@force2 + force + 1;
	s = (f - 2endLength) / (segments + .5);
	e = endLength + s / 4;
	GraphicsGroup[{
	LightDarkSwitched[Black, White], Flatten@{
	(* spring *)
	AbsoluteThickness[thickness], JoinForm["Round"], CapForm["Butt"], 
	LightDarkSwitched[Black, White], Line[Flatten[{{loc, loc + scale{ endLength, 0}}, 
	Flatten[Table[{loc + scale{ e + a s, -height}, loc + scale{e + a s + .5s, height}}, {a, 0, segments - 1}], 1], 
	{loc + scale{f - e, -height}, loc + scale{f - endLength, 0}, loc + scale{f - endLength, 0}, loc + scale{f, 0}}}, 1]], 
	White, CapForm["Round"], AbsoluteThickness[.5thickness], Line[{loc, loc + scale{ endLength, 0}}], 
	Table[Line[{loc + scale{ e + a s, -height}, loc + scale{e + a s + .5s, height}}], {a, 0, segments - 1}], 
	Line[{loc + scale{f - e, -height}, loc + scale{f - endLength, 0}}], Line[{loc + scale{f - endLength, 0}, loc + scale{f, 0}}], 
	(* cylinder *)
	Gray, Rectangle[loc + scale{.5f + .5, -1.3}, loc + scale{f - .25, -.7}], 
	Table[{If[t == thickness, LightDarkSwitched[Black, White], White], AbsoluteThickness[t], 
	Line[{loc + scale{f - 1.75, -.7}, loc + scale{f - .25, -.7}, loc + scale{f - .25, -1.3}, loc + scale{f - 1.75, -1.3}}], Line[{loc + scale{.5f + .5, -.7}, loc + scale{.5f + .5, -1.3}}], 
	Line[{loc + scale{0, -1}, loc + scale{.5f + .5, -1}}], 
	Line[{loc + scale{f, -1}, loc + scale{f - .25, -1}}]}, 
	{t, {thickness, .5thickness}}], 
	(* frame *)
	LightDarkSwitched[Black, White], AbsoluteThickness[thickness], Line[{loc, loc + scale{0, -1}}], 
	Line[{loc + scale{f, 0}, loc + scale{f, -1}}], 
	Line[{loc + scale{-endLength2, -.5}, loc + scale{0, -.5}}], Line[{loc + scale{f, -.5}, loc + scale{f + endLength2, -.5}}], 
	AbsolutePointSize[5], Point[loc + scale{-endLength2, -.5}], Point[loc + scale{f + endLength2, -.5}]}
	}]
	}]

バネとダンパのあるスライダクランクシステムの動きをアニメーション化する．ダンパとバネがあるため，動きは最終的に水平平面上で止まる：

Wolfram Language code:

Animate[
	rotate = (y1[t] /. SliderCrankSol2) * (180 / Pi);
	scale = 0.7;
	Module[{bLength = scale * 0.5, bfMidLength = scale * 0.5, bfLength = scale, bfAngle, bfLoc, sliderLoc}, 
	Graphics[{
	(* parts *)
	bfAngle = ArcSin[bLength Cos[(90 - rotate) °] / bfMidLength] * 180 / Pi;
	sliderLoc = {bLength Cos[rotate °] + bfMidLength Cos[bfAngle °], 0};
	slider[{.6, 0}, 0, .15, .4, .2, 3, sliderLoc - {.6, 0}, Hue[.05, .3, .9]], 
	bar[{0, 0}, rotate, 1, bLength, .03, Hue[.6, .2, .7]], 
	fixedHinge[{0, 0}, 0, .2], 
	bfLoc = {bLength Sin[(90 - rotate) °], bLength Cos[(90 - rotate) °]};
	barMidFlange[bfLoc, -bfAngle, 1, bfLength, .03, 1, 1, bfMidLength, Hue[.6, .2, .9]], 
	(* damper *)
	springDamper[bfLoc / 4, .08, 9sliderLoc[[1]] - 4.5, .1, 1.35, 2.5], 
	(* axis *)
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Arrowheads[Medium], 
	Arrow[{{-.25, 0}, {1.1, 0}}], (* x *)
	Arrow[{{0, -.05}, {0, .5}}](* y *)
	}, PlotRange -> {{-.15, 1.2}, {-.3, .5}}, ImageSize -> Medium]
	], 
	{t, 0, 20}, SaveDefinitions -> True, AnimationRunning -> False]

2リンク平面ロボットアーム

2リンク平面ロボットアームは，無摩擦の制約された機械システムである．このシステムは，2つ目のリンクの一端の軌道が指定されている一方で，最初のリンクの一端が固定されたままとなるようにジョイントで互いに接続された2つのリンクからできている．以下の図はこのシステムの概略である．

Wolfram Language code: Clear[m1, m2, g, constraints, vars, θ1, θ2, λ, vt];

このシステムはリンク間の角度を使って完全に定義することができる．このシステムの詳細はAscher et al. [ACPR94]でご覧いただきたい．

変数を定義する：

Wolfram Language code: vars = {θ1, θ2, λ};

ロボットアームの終点が従わなければ成らない軌道を定義する．これは

の位置を定義することに相当する．この動きは以下によって制約される：

Wolfram Language code: constraints = {x2 + y2 ^ 2 == 1.5};

運動の方程式は，力を分解しニュートンの法則を適用することで導かれる．微分方程式は角度に対して導かれる．したがって，デカルト座標およびは，角度で表される．

微分方程式を定義する：

Wolfram Language code:

x1 = L1 Cos[θ1[t]]; y1 = L1 Sin[θ1[t]];
x2 = x1 + L2 Cos[θ1[t] + θ2[t]]; y2 = y1 + L2 Sin[θ1[t] + θ2[t]];
eqns = {
	((L1^2 m1/3) + m2 (L1^2 + (L2^2/3) + L1 L2 Cos[θ2[t]])) θ1''[t] + m2 ((L2^2/3) + (1/2) L1 L2 Cos[θ2[t]]) θ2''[t] == -(1/2) g L1 m1 Cos[θ1[t]] - g m2 (L1 Cos[θ1[t]] + (1/2) L2 Cos[θ1[t] + θ2[t]]) - (-L1 Sin[θ1[t]] - L2 Sin[θ1[t] + θ2[t]] + 2 (L1 Cos[θ1[t]] + L2 Cos[θ1[t] + θ2[t]]) (L1 Sin[θ1[t]] + L2 Sin[θ1[t] + θ2[t]])) λ[t] + (1/2) L1 L2 m2 Sin[θ1[t] + θ2[t]] (2 θ1'[t] θ2'[t] + θ2'[t]^2), 
	m2 ((L2^2/3) + (1/2) L1 L2 Cos[θ2[t]]) θ1''[t] + (1/3) L2^2 m2 θ2''[t] == -(1/2) g L2 m2 Cos[θ1[t] + θ2[t]] - (-L2 Sin[θ1[t] + θ2[t]] + 2 L2 Cos[θ1[t] + θ2[t]] (L1 Sin[θ1[t]] + L2 Sin[θ1[t] + θ2[t]])) λ[t] - (1/2) L1 L2 m2 Sin[θ1[t] + θ2[t]] θ1'[t]^2
	};

系に関連付けられたパラメータを定義する：

Wolfram Language code: robotParams = {m1 -> 36, m2 -> 36, L1 -> 1, L2 -> 1, g -> 9.81};

系を解き，可視化する：

Wolfram Language code:

ics = {θ1[0] == 70. Degree, θ1'[0] == 0, θ2'[0] == 0};
roboticArmSol = First@NDSolve[{eqns, constraints, ics} /. robotParams, vars, {t, 0, 10}, Method -> {"IndexReduction" -> {Automatic}}, MaxSteps -> 100000]

Wolfram Language code: Row[Plot[Evaluate[#[t] /. roboticArmSol], {t, 0, 10}, PlotLabel -> #[t], PlotRange -> All, ImageSize -> Small]& /@ vars]

指定された曲線上を動くよう制約されたロボットアームの動きをアニメーション化する：

Wolfram Language code:

paramplot = ParametricPlot[Evaluate[{x2, y2} /. robotParams /. roboticArmSol], {t, 0, 10}];
Animate[
	{Subscript[θ, 1], Subscript[θ, 2]} = {θ1[t], θ2[t]} * (180 / Pi) /. roboticArmSol;
	Module[{bLength = 1, b2Length = 1, b2Angle, b2Loc, endLoc}, 
	Show[Graphics[{
	(* parts *)
	b2Angle = ArcSin[bLength Cos[(90 - Subscript[θ, 1]) °] / b2Length] * 180 / Pi;
	b2Loc = {bLength Sin[(90 - Subscript[θ, 1]) °], bLength Cos[(90 - Subscript[θ, 1]) °]};
	endLoc = b2Loc + {b2Length Sin[(Subscript[θ, 2] - (90 - Subscript[θ, 1]) + 180)°], -b2Length Cos[(Subscript[θ, 2] - (90 - Subscript[θ, 1]) + 180)°]};
	bar[{0, 0}, Subscript[θ, 1], 1, bLength, .05, Hue[.6, .2, .9]], 
	fixedHinge[{0, 0}, 0, .2], 
	bar[b2Loc, Subscript[θ, 1] + Subscript[θ, 2], 1, b2Length, .05, Hue[.6, .2, .9]], 
	(* axis *)
	LightDarkSwitched[Black, White], AbsoluteThickness[.5], Arrowheads[Medium], 
	Arrow[{{0, 0}, {.5, 0}}], (* x *)
	LightDarkSwitched[Black, White], Text[Style["x", 12, Italic], {0, .54}], 
	Arrow[{{0, 0}, {0, .5}}], (* y *)
	LightDarkSwitched[Black, White], Text[Style["y", 12, Italic], {.54, 0}]
	(* dimensions *)
	}, ImageSize -> Medium, PlotRange -> {{-2, 2}, {-2, 2}}], paramplot]
	], 
	{t, 0, 10}, AnimationRunning -> False, SaveDefinitions -> True]

n 振り子問題

振り子系を生成する関数を構築する：

Wolfram Language code:

Clear[x, y, yp, c, eqns, vars, ics, i];
NPendulumDAESystem[np_, length_, gravity_, t_] := Module[{eqns, vars, ics, i}, 
	x[0][t] = y[0][t] = 0; c[np + 1][t] = 0;
	eqns = Flatten[Table[(⁠|                                                                                               |
| --------------------------------------------------------------------------------------------- |
| x[i]''[t] == c[i][t] (x[i][t] - x[i - 1][t]) - c[i + 1][t] (x[i + 1][t] - x[i][t])            |
| y[i]''[t] == c[i][t] (y[i][t] - y[i - 1][t]) -  	c[i + 1][t] (y[i + 1][t] - y[i][t]) - gravity |
| (x[i][t] - x[i - 1][t])^2 + (y[i][t] - y[i - 1][t])^2 == length^2                              |⁠), 
	{i, np}]];
	vars = Flatten[Table[{x[i], y[i], c[i]}, {i, np}]];
	ics = Flatten[Table[{x[i][0] == i, y[i]'[0] == 1.1}, {i, np}]];
	{eqns, vars, ics}]

2つの振り子系を生成する：

Wolfram Language code:

Clear[x, y, c];
{eqns, vars, ics} = NPendulumDAESystem[2, 1, 1, t];

振り子系を解き，可視化する：

Wolfram Language code:

sol = Quiet@NDSolve[{eqns, ics}, vars, {t, 0, 30}, Method -> {"IndexReduction" -> {Automatic, "IndexGoal" -> 0}, "EquationSimplification" -> "Residual"}, MaxSteps -> 100000];
Row[Plot[Evaluate[#[t] /. sol], {t, 0, 30}, PlotLabel -> #, ImageSize -> Small]& /@ vars]

2つの振り子の棒のカオス軌道を可視化する：

Wolfram Language code:

Show[ParametricPlot[Evaluate[{x[1][t], y[1][t]} /. sol], {t, 0, 30}, PlotStyle -> StandardBlue], 
	ParametricPlot[Evaluate[{x[2][t], y[2][t]} /. sol], {t, 0, 30}, PlotStyle -> StandardRed], 
	Graphics[{{LightDarkSwitched[Black, White], CapForm["Round"], Thickness[0.03], Line[{{0, 0}, {x[1][10], y[1][10]}}]}, {Orange, Thickness[0.03], CapForm["Round"], Line[{{x[1][10], y[1][10]}, {x[2][10], y[2][10]}}]}}] /. sol, 
	PlotRange -> All, ImageSize -> 300]

トランジスタ増幅回路

トランジスタ増幅回路のパラメータを定義する：

Wolfram Language code:

Clear[ve, params, vb, r0, r1, r2, r3, r4, r5, α, β, c1, c2, c3, v23, node1, node2, node3, node4, node5, ics, sol, v1, v2, v3, v4, v5];
ve[t_] := (4 / 10) Sin[200 Pi t];
params = {vb -> 6, r0 -> 1000, r1 -> 9000, r2 -> 9000, r3 -> 9000, r4 -> 9000, r5 -> 9000, α -> 99 / 100, β -> 10 ^ -6, c1 -> 10 ^ -6, c2 -> 2 10 ^ -6, c3 -> 3 10 ^ -6};
v23[t_] = β (Exp[(v2[t] - v3[t]) / .026] - 1);
node1 = c1 (v2'[t] - v1'[t]) == v1[t] / r0 - ve[t] / r0;
node2 = c1(v1'[t] - v2'[t]) == (1 - α)v23[t] + v2 [t] / r1 + v2[t] / r2 - vb / r2;
node3 = c2 v3'[t] == v23[t] - v3[t] / r3;
node4 = c3(v4'[t] - v5'[t]) == vb / r4 - v4[t] / r4 - α v23[t];
node5 = c3(v4'[t] - v5'[t]) == v5[t] / r5;
ics = {v1[0] == 0, v2[0] == vb / 2, v3[0] == vb / 2, v4[0] == vb, v5[0] == 0};

微分代数方程式系を解き，可視化する：

Wolfram Language code:

sol = First@NDSolve[{node1, node2, node3, node4, node5, ics} //. params, {v1, v2, v3, v4, v5}, {t, 0, 0.2}, Method -> {"EquationSimplification" -> "Residual"}, MaxSteps -> 100000] ;

	
Plot[Evaluate[{v3[t], v4[t]} /. sol], {t, 0, 0.2}, PlotRange -> All]

Akzo-Nobelの化学反応

次のシステムは，2つの物質FLBとZHUが混合され，その間常にがシステムに加えられるという科学プロセスを記述している．この化学物質名は架空のものである．このシステムとそれに関連するダイナミクスは[N93]に記載されている．反応方程式は下に示す．それぞれの反応と関連して反応速度があり，指定されている．代数的制約は，FLBとZHU，およびその混合体が一定の関係を保つよう要求している最後の反応から派生している．

化学反応の反応速度とパラメータを定義する：

Wolfram Language code:

Fin = klA (pCO2 / H - CO2[t]);
r1 = k1 FLB[t]^2 CO2[t]^1 / 2;
r2 = k2 FLBT[t] ZHU[t];
r3 = (k2 / KK) FLB[t] ZLA[t];
r4 = k3 FLB[t] ZHU[t]^2 CO2[t];
r5 = k4 FLB[t] ZHU[t];
	
params = {k1 -> 18.7, k2 -> 0.58, k3 -> 0.09, k4 -> 0.42, KK -> 34.4, klA -> 3.3, Ks -> 115.83, pCO2 -> 0.9, H -> 737};

系の速度方程式を生成する：

Wolfram Language code:

eqns = {FLB'[t] == -2 r1 + r2 - r3 - r4, CO2'[t] == -0.5 r1 - r4 - 0.5r5 + Fin, FLBT'[t] == r1 - r2 + r3, ZHU'[t] == -r2 + r3 - 2 r4, ZLA'[t] == r2 - r3 + r5};
eqEqn = Ks FLB[t] ZHU[t] == FLBZHU[t];
	
ic = {FLB[0] == 0.444, CO2[0] == 0.00123, FLBT[0] == 0, ZHU[0] == 0.007, ZLA[0] == 0};

系を解き，可視化する：

Wolfram Language code:

sol = NDSolve[{eqns, eqEqn, ic} /. params, {FLB, ZHU, CO2, ZLA}, {t, 0, 200}];
Plot[Evaluate[#[t] /. sol], {t, 0, 200}, PlotRange -> All, PlotLabel -> #[t]]& /@ {FLB, ZHU, CO2, ZLA}

車軸問題

車軸モデルは，でこぼこ道での車軸の動作をシミュレーションするために設計された簡単な多体系である．この系の概略を以下に示す．このモデルは2つの車輪からできている．左の車輪は平坦な面の上を動き，右の車輪はでこぼこの表面上を動くものとされている．でこぼこの表面は正弦曲線としてモデル化される．両方の車輪は，以下に示すように質量の無視できる2つのバネを通して，その動きを車のシャーシに伝達する．車のシャーシは質量Mのバーとしてモデル化される．このモデルの詳細はMazzia et. al. [MM08]で見ることができる．これは指数3の系である．

Wolfram Language code: Clear[yB, xL, yL, xB, L, m, Lleft, Lright, ϵ, eqns, ics, ω, carParams, params, vars, carAxis, M];

左側のタイヤは常に平面にあり，右側のタイヤは定期的にでこぼこの上を走る．

右側のタイヤがでこぼこの上を走るときの動きを定義する：

Wolfram Language code: yB[t_] = h Sin[ω t];

タイヤとタイヤの間の距離は，常に固定されたままでなければならず，車のシャーシは常に固定された長さでなければならない．

Wolfram Language code: constraints = {xL[t] xB[t] + yL[t] yB[t] == 0, (xL[t] - xR[t])^2 + (yL[t] - yR[t])^2 == L^2};

ニュートンの法則を使って，運動の方程式を導く：

Wolfram Language code:

m = ((ϵ^2M)/2);
xB[t_] = Sqrt[L^2 - yB[t]^2];
Lleft = Sqrt[xL[t]^2 + yL[t]^2];
Lright = Sqrt[(xR[t] - xB[t])^2 + ( yR[t] - yB[t])^2];
eqns = (⁠|                                                                                     |
| ----------------------------------------------------------------------------------- |
| m xL''[t] == ((L0 - Lleft)xL[t]/Lleft) + λ1[t] xB[t] + 2 λ2[t] (xL[t] - xR[t])       |
| m yL''[t] == ((L0 - Lleft) yL[t]/Lleft ) + λ1[t] yB[t] + 2 λ2[t] (yL[t] - yR[t]) - m |
| m xR''[t] == (L0 - Lright)((xR[t] - xB[t])/Lright) - 2 λ2[t] (xL[t] - xR[t])         |
| m yR''[t] == (L0 - Lright) ((yR[t] - yB[t])/Lright) - 2 λ2[t] (yL[t] - yR[t]) - m    |⁠);
ics = {xL[0] == 0, yL[0] == 1 / 2, yR[0] == 1 / 2, yL'[0] == 0, yR'[0] == 0};

系に関連付けられたパラメータを定義する：

Wolfram Language code:

carParams = {L -> 1, L0 -> (1/2), ϵ -> 10^-2, M -> 10, h -> 10^-1, ω -> 10};
vars = {xL, yL, xR, yR, λ1, λ2};

解いて，解を可視化する：

Wolfram Language code:

carAxis = NDSolve[{eqns, constraints, ics} /. carParams, vars, {t, 0, 3}, Method -> {"IndexReduction" -> Automatic}, AccuracyGoal -> 7] ;

Wolfram Language code:

Animate[
	With[{t = tt}, 
	Show[
	(* ground *)
	Plot[Evaluate[{h * Exp[-5 * (x - 1.2) ^ 2] * Sin[ω t]} - 0.14 /. carParams], {x, -0.9, 2}, Filling -> -0.5, FillingStyle -> Hue[.1, .2, .4], Prolog -> {Polygon[{Scaled[{0, 0}], Scaled[{1, 0}], Scaled[{1, 1}], Scaled[{0, 1}]}, VertexColors -> {Hue[.6, .2, 1], Hue[.6, .2, 1], Hue[.6, .7, .5], Hue[.6, .7, .5]}]}, PlotStyle -> Hue[.1, .4, .2], PlotRange -> {{-0.9, 1.9}, {-0.5, 1.5}}, Frame -> True, Axes -> False], 
	Graphics[{
	(* Axle *)
	{Thickness[0.03], Hue[.6, 1, .3], Line[{{xL[t] - 0.2, yL[t]} / 2, {(xB[t]) + 0.2, (yR[t] - 0 * yB[t]) / 2}}]}, 
	(* tires *)
	{Thickness[0.09], GrayLevel[.2], Line[{{Mean[{0 - 0.22, xL[t] - 0.22}], 0}, {Mean[{0 - 0.22, xL[t] - 0.22}], yL[t]}}]}, 
	{Thickness[0.09], GrayLevel[.2], Line[{{Mean[{xB[t] + 0.22, xR[t] + 0.22}], yB[t]}, {Mean[{xB[t] + 0.22, xR[t] + 0.22}], yR[t]}}]}, 
	(* springs *)
	{Thickness[0.025], Hue[.1, 1, .9], Line[{{0 + 0.01, yL[0] / 2 + 0.03}, {xL[t] + 0.01, yL[t] + 0.2}}]}, 
	{Thickness[0.025], Hue[.1, 1, .9], Line[{{xB[t], yB[t] / 2 + 0.28}, {xR[t] - 0.01, yR[t] + 0.2}}]}, 
	(* chassis *)
	{Thickness[0.01], Hue[.6, 1, .3], Line[{{xL[t] - 0.03, yL[t] + 0.2}, {xR[t] + 0.03, yR[t] + 0.2}}]}, 
	(* body *)
	{GrayLevel[1, .3], EdgeForm[{GrayLevel[1, .5]}], GeometricTransformation[
	{Polygon[.5 + .0045{{-185., -105.}, {-195., -100.}, {-200., -95.}, {-215., 0}, {-215., 45.}, {-210., 50.}, {-205., 50.}, {-200., 65.}, {-190., 75.}, {-220., 75.}, {-230., 80.}, {-230., 95.}, {-225., 100.}, {-205., 105.}, {-190., 105.}, {-185., 100.}, {-180., 90.}, {-175., 90.}, {-160., 125.}, {-145., 155.}, {-135., 165.}, {-120., 172.}, {-75., 177.}, {0., 180.}, {75., 177.}, {120., 172.}, {135., 165.}, {145., 155.}, {160., 125.}, {175., 90.}, {180., 90.}, {185., 100.}, {190., 105.}, {205., 105.}, {225., 100.}, {230., 95.}, {230., 80.}, {220., 75.}, {190., 75.}, {200., 65.}, {205., 50.}, {210., 50.}, {215., 45.}, {215., 0}, {200., -95.}, {195., -100.}, {185., -105.}}], Polygon[.5 + .0045{{-145., 80.}, {-155., 90.}, {-145., 125.}, {-130., 155.}, {-120., 160.}, {-55., 165.}, {0., 165.}, {55., 165.}, {120., 160.}, {130., 155.}, {145., 125.}, {155., 90.}, {145., 80.}}], Polygon[.5 + .0045{{-190., 20.}, {-185., 25.}, {-150., 25.}, {-110., 20.}, {-100., 15.}, {-100., 10.}, {-110., 0.}, {-140., -0.08875739644975056}, {-160., -10.}, {-175., -10.}, {-185., -5.}, {-190., 5.}}], Polygon[.5 + .0045{{110., -0.08875739644975056}, {140., -0.08875739644975056}, {160., -10.}, {175., -10.}, {185., -5.}, {190., 5.}, {190., 20.}, {185., 25.}, {150., 25.}, {110., 20.}, {100., 15.}, {100., 10.}}], Polygon[.5 + .0045{{-85., 15.}, {-75., 20.}, {-25., 25.}, {0., 25.}, {25., 25.}, {75., 20.}, {85., 15.}, {90., 5.}, {90., -5.}, {85., -30.}, {75., -60.}, {65., -70.}, {55., -75.}, {-55., -75.}, {-65., -70.}, {-75., -60.}, {-85., -30.}, {-90., -5.}, {-90., 5.}}], Polygon[.5 + .0045{{-85., -50.}, {-90., -40.}, {-190., -40.}, {-205., -25.}, {-195., -50.}}], Polygon[.5 + .0045{{85., -50.}, {195., -50.}, {205., -25.}, {190., -40.}, {90., -40.}}]}, RotationTransform[ArcTan[xR[t] - xL[t], yR[t] - yL[t]]]]}
	} /. carParams /. carAxis], ImageSize -> 400
	]
	]
	, {tt, 0, .2Pi}, SaveDefinitions -> True, AnimationRunning -> False]

楕円型偏微分方程式と放物型偏微分方程式の1Dでの結合

次は，非線形の形式で結合された放物型偏微分方程式と楕円型偏微分方程式を解くための合成例題である．偏微分方程式の空間的離散性により，かなり多数の代数的変数を含む系になる．方程式は，系に対して解析的解が存在するように設定される．結合系の詳細は以下で与えられる．

線の方法を使って偏微分方程式の空間的要素を離散化し，支配方程式から常微分方程式系を構築する．これらが微分方程式を形成する：

Wolfram Language code:

Clear[u, v, nx, eqn1, eqn2, ic1, ic2, vars, eqn1, eqn2, ic1, ic2, vars, sol, PDEsol, order, p, npts];
npts = 100; order = 12; p = 10;
nx = Range[0, 1, 1 / (npts - 1)];
{ddx, d2dx2} = Map[NDSolve`FiniteDifferenceDerivative[Derivative[#], nx, "DifferenceOrder" -> order]["DifferentiationMatrix"]&, {1, 2}] ;
u = Array[u$[#][t]&, npts]; v = Array[v$[#][t]&, npts];
	
eqn1 = Thread[D[u, t] - D[v, t] == d2dx2.u + u (ddx.v) - v (ddx.u) + (2 π t)^2v];
eqn2 = Thread[d2dx2.v + (2 π t)^2 u == 0];

系の境界条件を指定する．これらが系の代数方程式を形成する：

Wolfram Language code:

eqn1[[{1, -1}]] = Thread[u[[{1, -1}]] == v[[{1, -1}]]] ;
eqn2[[1]] = v[[1]] == 0;
eqn2[[-1]]  =   ddx[[-1]].v == 2π t Cos[2 π t] ;

系の初期条件を指定する：

Wolfram Language code:

ic1 = Thread[u == 0] /. t -> 0;
ic2 = Thread[v == 0] /. t -> 0;

偏微分方程式を微分代数方程式系として解く：

Wolfram Language code:

vars = Flatten[{u, v}] /. x__[t] :> x ; 
PDEsol = First[NDSolve[{eqn1, eqn2, ic1, ic2}, vars, {t, 0, 3}]];

解を可視化する：

Wolfram Language code:

time = Flatten@First[vars[[1]] /. PDEsol];
{usol, vsol} = Flatten[Table[Transpose[{nx, ConstantArray[ti, npts], # /. PDEsol}] /. t -> ti, {ti, time[[1]], time[[2]], 0.1}], 1]& /@ {u, v};
GraphicsRow[ListPlot3D[#, PlotRange -> All]& /@ {usol, vsol}, ImageSize -> 500]

多種の食物網モデル（楕円型偏微分方程式と放物型偏微分方程式の2Dでの結合）

この例は，多種の食物網のモデルである[B86]．典型的な捕食被食問題とは異なり，この問題には空間的要素が存在する．捕食被食関係は一連の拡散パラメータを使って2Dの空間領域上で定義される．一般的な結合偏微分方程式は以下で与えられる．

上の系では，最初の種が被食者，次の種が捕食者である．ここではの特殊な場合が考慮される．つまり，1種の捕食者と1種の被食者がいるということである．空間領域は20×20のメッシュで離散化される．最終的な微分代数方程式系には2×400個の変数があることになる．捕食者の一様分布は最初に指定される．これは矛盾した初期状態になるが，組込みのアルゴリズムで修正される．

モデルに必要なパラメータを定義する：

Wolfram Language code:

Clear[p, α, β, a, bv, di, bi, ri, nx, grid, nL, order, x, y, dfdx, dfdy, d2fdx2, d2fdy2, Lapf, species, R, eqns, ics, vars, speciesSol, c];
p = 1; α = 50; β = 1000;
	
bv = Table[If[i ≤ p, 1, -1], {i, 2p}];
di = Table[If[i ≤ p, 1, 0.5], {i, 2 p}] ; 
a = SparseArray[{{i_, i_} -> -1, {i_, j_} /; (i ≤ p  && j > p) -> -0.5 10 ^ -6, {i_, j_} /; (i > p  && j ≤ p) -> 10 ^ 4}, {2 p, 2 p}] ;
	
bi[x_, y_] := (1 + α x y   + β Sin[4 Pi x] Sin[4 Pi y]) * bv;
ri[{x_, y_}, c_List] :=   c * (bi[x, y] + a.c);
	
nx = N@Range[0, 1, 1 / 19]; 
grid = Flatten[Outer[List, nx, nx], 1];
nL = Length[grid];
order = 2;
	
{left, right, bot, top} = Flatten[Position[grid, #]]& /@ {{0., y_}, {1., y_}, {x_, 0.}, {x_, 1.}};

最初と2つ目の導関数に対して有限差分微分行列を生成する．それからラプラス行列を生成する：

Wolfram Language code:

{dfdx, dfdy, d2fdx2, d2fdy2} = Map[ NDSolve`FiniteDifferenceDerivative[#, {nx, nx}, "DifferenceOrder" -> order]["DifferentiationMatrix"]&, {{1, 0}, {0, 1}, {2, 0}, {0, 2}}];
Lapf = d2fdx2 + d2fdy2;

食物網の種を定義する：

Wolfram Language code: species = Transpose[Table[Array[c[i][#][t]&, nL], {i, 2 p}]];

このモデルの方程式を生成する：

Wolfram Language code:

R = Apply[ri, MapThread[List, {grid, species}], 1] ; 
eqns  = Table[Thread[If[i ≤ p, 1, 0] * D[species[[All, i]], t] == R[[All, i]] + di[[i]] * Lapf.species[[All, i]]], {i, 2 p}];
Do[
	eqns[[i]][[left]] = Thread[dfdx[[left]].species[[All, i]] == 0];
	eqns[[i]][[right]] = Thread[dfdx[[right]].species[[All, i]] == 0];
	eqns[[i]][[bot]] = Thread[dfdy[[bot]].species[[All, i]] == 0];
	eqns[[i]][[top]] = Thread[dfdy[[top]].species[[All, i]] == 0];
	, {i, 2 p}];

モデルの初期条件を定義する：

Wolfram Language code:

ics = Flatten[Table[Thread[(species[[All, i]]  /. t -> 0) == If[i ≤ p, 10 + i * (grid /. {x_, y_} :> (16x (1 - x) y (1 - y)) ^ 2), 10 ^ 5]], {i, 2 p}]];

食物網モデルをシミュレーションする：

Wolfram Language code:

vars = Flatten[species] /. x___[t] :> x;
speciesSol = First[NDSolve[{eqns, ics}, vars, {t, 0, 1}]];

修正された初期条件をプロットする．捕食者の分布が変化している：

Wolfram Language code:

csol = (species /. t -> 0) /. speciesSol ; 
Row[ListPlot3D[Join[grid, csol[[All, {#}]], 2], ImageSize -> 200]& /@ Range[2 p]]

および

における解を可視化すると，指定されたパラメータについて定常状態に速く達することが分かる：

Wolfram Language code: speciesRange = {{0, 750}, {0, 7.5 10 ^ 6}};

Wolfram Language code:

csol = (species /. t -> 0.0001) /. speciesSol ;
Row[Table[ ListPlot3D[Join[grid, csol[[All, {i}]], 2], ImageSize -> 200], {i, 2p}]]

Wolfram Language code:

csol = (species /. t -> 0.1) /. speciesSol ;
Row[Table[ ListPlot3D[Join[grid, csol[[All, {i}]], 2], ImageSize -> 200, PlotRange -> {{0, 1}, {0, 1}, speciesRange[[i]]}], {i, 2p}]]

Wolfram Language code:

csol = (species /. t -> 1) /. speciesSol ;
Row[Table[ ListPlot3D[Join[grid, csol[[All, {i}]], 2], ImageSize -> 200, PlotRange -> {{0, 1}, {0, 1}, speciesRange[[i]]}], {i, 2p}]]

Andrewsの圧搾器メカニズム

Andrewsの圧搾器メカニズムは7つの異なるリンクを含む多体系である．さまざまなリンクの次元およびその指定は以下にある通りである．この系は7つの角度を使って完全に定義することができる．角度は次の系で q₁–q₇としてラベル付けされている．また，リンクの次元と接続に関連付けられた6つの代数的制約がある．これにより14次元の指数3の系になる．このシステムの詳細はHairerとWanner[HW96]の論文を参照されたい．

さまざまなリンクに関連付けられた次元は以下に示す通りである．7つのリンクがK₁–K₇でラベル付けされている．

モデルで使われるパラメータの表記とその説明は，HairerとWannerの論文[HW96]に記載されている．

パラメータのいくつかをクリアする：

Wolfram Language code: Clear[params0, params1, params2, params3, varQ, varL, vars, m, L, force, f, g, r, u, q, i, a, z, Subscript];

パラメータを定義する：

Wolfram Language code:

params0 = {Subscript[x, a] -> -0.06934, Subscript[y, a] -> -0.00227, Subscript[x, b] -> -0.03635, Subscript[y, b] -> 0.03273, Subscript[x, c] -> 0.014, Subscript[y, c] -> 0.072};
params1 = {mom -> 0.033, d -> 0.028, Subscript[e, a] -> 0.01421, Subscript[r, r] -> 0.007, Subscript[s, a] -> 0.01874, Subscript[s, d] -> 0.02, Subscript[t, b] -> 0.00916, Subscript[u, b] -> 0.00449, Subscript[d, a] -> 0.0115, Subscript[z, f] -> 0.02, Subscript[r, a] -> 0.00092, Subscript[s, b] -> 0.01043, Subscript[z, t] -> 0.04, u -> 0.04, Subscript[c, 0] -> 4530, e -> 0.02, Subscript[f, a] -> 0.01421, Subscript[s, s] -> 0.035, Subscript[s, c] -> 0.018, Subscript[t, a] -> 0.02308, Subscript[u, a] -> 0.01228, Subscript[L, 0] -> 0.07785};
params2 = {Subscript[m, 1] -> 0.04325, Subscript[m, 2] -> 0.00365, Subscript[m, 3] -> 0.02373, Subscript[m, 4] -> 0.00706, Subscript[m, 5] -> 0.07050, Subscript[m, 6] -> 0.00706, Subscript[m, 7] -> 0.05498};
params3 = {Subscript[i, 1] -> 2.194 10^-6, Subscript[i, 2] -> 4.410 10^-7, Subscript[i, 3] -> 5.255 10^-6, Subscript[i, 4] -> 5.667 10^-7, Subscript[i, 5] -> 1.169 10^-5, Subscript[i, 6] -> 5.667 10^-7, Subscript[i, 7] -> 1.912 10^-5};
andrewsParams = Join[params0, params1, params2, params3];

系の運動は7つの角度を使って定義することができる．この角度はここでと定義されている．運動の方程式は，形式の方程式系になるハミルトニアン原理を使って導かれる．ここでは質量行列，は力のベクトル，は制約条件のベクトル，λ はラグランジュの乗数である．

角度とラグランジュの乗数に対する変数を定義する：

Wolfram Language code:

vQ = {Subscript[q, 1][t], Subscript[q, 2][t], Subscript[q, 3][t], Subscript[q, 4][t], Subscript[q, 5][t], Subscript[q, 6][t], Subscript[q, 7][t]};
vL = {Subscript[λ, 1][t], Subscript[λ, 2][t], Subscript[λ, 3][t], Subscript[λ, 4][t], Subscript[λ, 5][t], Subscript[λ, 6][t]};
vars = Head /@ Join[vQ, vL];

系に関連付けられた質量行列

を定義する：

Wolfram Language code:

M = SparseArray[{
	{1, 1} -> Subscript[m, 1] Subscript[r, a]^2 + Subscript[m, 2] (Subscript[r, r]^2 - 2 Subscript[d, a] Subscript[r, r] Cos[Subscript[q, 2][t]] + Subscript[d, a]^2) + Subscript[i, 1] + Subscript[i, 2], 
	{1, 2} -> Subscript[m, 2](Subscript[d, a]^2 - Subscript[d, a] Subscript[r, r] Cos[Subscript[q, 2][t]]) + Subscript[i, 2], 
	{2, 1} -> Subscript[m, 2] (Subscript[d, a]^2 - Subscript[d, a] Subscript[r, r] Cos[Subscript[q, 2][t]]) + Subscript[i, 2], 
	{2, 2} -> Subscript[m, 2] Subscript[d, a]^2 + Subscript[i, 2], {3, 3} -> Subscript[m, 3] (Subscript[s, a]^2 + Subscript[s, b]^2) + Subscript[i, 3], {4, 4} -> Subscript[m, 4] (e - Subscript[e, a])^2 + Subscript[i, 4], 
	{5, 4} -> Subscript[m, 4] ((e - Subscript[e, a])^2 + Subscript[z, t] (e - Subscript[e, a]) Sin[Subscript[q, 4][t]]) + Subscript[i, 4], 
	{4, 5} -> Subscript[m, 4] ((e - Subscript[e, a])^2 + Subscript[z, t] (e - Subscript[e, a]) Sin[Subscript[q, 4][t]]) + Subscript[i, 4], 
	{5, 5} -> Subscript[m, 4](Subscript[z, t]^2 + 2 Subscript[z, t] (e - Subscript[e, a]) Sin[Subscript[q, 4][t]] + (e - Subscript[e, a])^2) + Subscript[m, 5] (Subscript[t, a]^2 + Subscript[t, b]^2) + Subscript[i, 4] + Subscript[i, 5], 
	{6, 6} -> Subscript[m, 6] (Subscript[z, f] - Subscript[f, a])^2 + Subscript[i, 6], {6, 7} -> Subscript[m, 6] ((Subscript[z, f] - Subscript[f, a])^2 - u (Subscript[z, f] - Subscript[f, a]) Sin[Subscript[q, 6][t]]) + Subscript[i, 6], 
	{6, 7} -> Subscript[m, 6]((Subscript[z, f] - Subscript[f, a])^2 - u (Subscript[z, f] - Subscript[f, a]) Sin[Subscript[q, 6][t]]) + Subscript[i, 6], 
	{7, 6} -> Subscript[m, 6] ((Subscript[z, f] - Subscript[f, a])^2 - u (Subscript[z, f] - Subscript[f, a]) Sin[Subscript[q, 6][t]]) + Subscript[i, 6], 
	{7, 7} -> Subscript[m, 6]((Subscript[z, f] - Subscript[f, a])^2 - 2 u (Subscript[z, f] - Subscript[f, a]) Sin[Subscript[q, 6][t]] + u^2) + Subscript[m, 7] (Subscript[u, a]^2 + Subscript[u, b]^2) + Subscript[i, 6] + Subscript[i, 7]}
	];

系に関連付けられたさまざまな力

を定義する：

Wolfram Language code:

Subscript[x, d] = Subscript[s, d] Cos[Subscript[q, 3][t]] + Subscript[s, c] Sin[Subscript[q, 3][t]] + Subscript[x, b];
Subscript[y, d] = Subscript[s, d] Sin[Subscript[q, 3][t]] - Subscript[s, c] Cos[Subscript[q, 3][t]] + Subscript[y, b];
L = Sqrt[(Subscript[x, d] - Subscript[x, c])^2 + (Subscript[y, d] - Subscript[y, c])^2];
force = -Subscript[c, 0] (L - Subscript[L, 0]) / L;
Subscript[f, x] = force (Subscript[x, d] - Subscript[x, c]);
Subscript[f, y] = force (Subscript[y, d] - Subscript[y, c]);

Wolfram Language code:

F = (⁠|                                                                         |
| ----------------------------------------------------------------------- |
| mom - m2 da rr q2'[t] (q2'[t] + 2 q1'[t]) Sin[q2[t]]                    |
| m2 da rr q1'[t]^2Sin[q2[t]]                                             |
| fx (sc Cos[q3[t]] - sd Sin[q3[t]]) + fy (sd Cos[q3[t]] + sc Sin[q3[t]]) |
| m4 zt (e - ea) q5'[t]^2Cos[q4[t]]                                       |
| -m4 zt (e - ea) q4'[t] (q4'[t] + 2 q5'[t]) Cos[q4[t]]                   |
| -m6 u (zf - fa) q7'[t]^2Cos[q6[t]]                                      |
| m6 u (zf - fa) q6'[t] (q6'[t] + 2 q7'[t]) Cos[q7[t]]                    |⁠)  //Flatten;

さまざまなリンクに関連付けられた制約条件ベクトル

を定義する：

Wolfram Language code:

g = (⁠|                                                                                  |
| -------------------------------------------------------------------------------- |
| rr Cos[q1[t]] - d Cos[q1[t] + q2[t]] - ss Sin[q3[t]] - xb                        |
| rr Sin[q1[t]] - d Sin[q1[t] + q2[t]] + ss Cos[q3[t]] - yb                        |
| rr Cos[q1[t]] - d Cos[q1[t] + q2[t]] - e Sin[q4[t] + q5[t]] - zt Cos[q5[t]] - xa |
| rr Sin[q1[t]] - d Sin[q1[t] + q2[t]] + e Cos[q4[t] + q5[t]] - zt Sin[q5[t]] - ya |
| rr Cos[q1[t]] - d Cos[q1[t] + q2[t]] - zf Cos[q6[t] + q7[t]] - u Sin[q7[t]] - xa |
| rr Sin[q1[t]] - d Sin[q1[t] + q2[t]] - zf Sin[q6[t] + q7[t]] + u Cos[q7[t]] - ya |⁠)//Flatten;

方程式

とパラメータの特定の値を持つ制約方程式を構築する：

Wolfram Language code: eqns = Thread /@ {M.D[vQ, {t, 2}] == F - Transpose[D[g, {vQ}]].vL , g == 0} /. andrewsParams;

特定の角度の集合で休止する，系の初期状態を定義する：

Wolfram Language code:

ics = Thread /@ {vQ == {-0.06171389001427, 0, 0.4552798191630704, 0.22266839016588588, 0.4873649795438425, -0.222668390165, 1.2305474445498212}, D[vQ, t] == 0} /. t -> 0;

系を解き，角度と力をプロットする．リンクの運動は時間の経過に伴い速くなる．直線はリンクの色に基づいてカラーコードされている：

Wolfram Language code:

solAndrews = First@NDSolve[{eqns, ics}, vars, {t, 0, 0.05}, Method -> {"IndexReduction" -> {Automatic, "ConstraintMethod" -> None}, "EquationSimplification" -> "Residual"}, MaxSteps -> ∞];

Wolfram Language code:

Plot[Evaluate[Mod[vQ /. solAndrews, 2 Pi, -Pi]], {t, 0, 0.05}, ImageSize -> Medium, Frame -> True, FrameLabel -> {"Time (t)", "Angles (q)"}, 
	PlotRange -> {Automatic, {-0.7, 1.5}}, PlotStyle -> {{Dotted, StandardBrown}, {Dotted, StandardGreen}, {Thick, StandardBlue}, {Thick, StandardBrown}, {Thick, LightRed}, {Thick, StandardMagenta}, {Thick, StandardOrange}}, PlotLegends -> (Style[#, 14]& /@ (vQ /. x___[t_] :> x))]

Wolfram Language code:

Plot[Evaluate[F  /. andrewsParams /. solAndrews], {t, 0, 0.05}, ImageSize -> Medium, Frame -> True, FrameLabel -> {"Time (t)", "Forces (SubscriptBox[F, i])"}, 
	PlotRange -> {Automatic, {-0.2, 0.2}}, 
	PlotStyle -> {{Dotted, StandardBrown}, {Dotted, StandardGreen}, {Thick, StandardBlue}, {Thick, StandardBrown}, {Thick, LightRed}, {Thick, StandardMagenta}, {Thick, StandardOrange}}, PlotLegends -> (Style[StringForm["SubscriptBox[K, `1`]", #], 14]& /@ Range[7])]

このメカニズムの動きを可視化する：

Wolfram Language code:

Animate[
	p1 = {Subscript[r, r] Cos[Subscript[q, 1][tt]], Subscript[r, r] Sin[Subscript[q, 1][tt]]};
	p2 = {Subscript[r, r] Cos[Subscript[q, 1][tt]] - d Cos[Subscript[q, 1][tt] + Subscript[q, 2][tt]], Subscript[r, r]Sin[Subscript[q, 1][tt]] - d Sin[Subscript[q, 1][tt] + Subscript[q, 2][tt]]};
	p3 = {Subscript[x, b] + Subscript[s, c] Sin[Subscript[q, 3][tt]] + Subscript[s, d] Cos[Subscript[q, 3][tt]], Subscript[y, b] - Subscript[s, c] Cos[Subscript[q, 3][tt]] + Subscript[s, d] Sin[Subscript[q, 3][tt]]};
	p3c = {Subscript[x, b] + Subscript[s, a] Sin[Subscript[q, 3][tt]] + Subscript[s, b] Cos[Subscript[q, 3][tt]], Subscript[y, b] - Subscript[s, a] Cos[Subscript[q, 3][tt]] + Subscript[s, b] Sin[Subscript[q, 3][tt]]};
	p4 = {Subscript[x, a] + Subscript[z, t]Cos[Subscript[q, 5][tt]], Subscript[y, a]  + Subscript[z, t] Sin[Subscript[q, 5][tt]]};
	p5 = {Subscript[x, a] + Subscript[t, a] Cos[Subscript[q, 5][tt]] - Subscript[t, b] Sin[Subscript[q, 5][tt]], Subscript[y, a] + Subscript[t, a] Sin[Subscript[q, 5][tt]] + Subscript[t, b] Cos[Subscript[q, 5][tt]]};
	p6 = {Subscript[x, a] + u Sin[Subscript[q, 7][tt]], Subscript[y, a] - u Cos[Subscript[q, 7][tt]]};
	p7 = {Subscript[x, a] + Subscript[u, a] Sin[Subscript[q, 7][tt]] - Subscript[u, b] Cos[Subscript[q, 7][tt]], Subscript[y, a] - Subscript[u, a] Cos[Subscript[q, 7][tt]] - Subscript[u, b] Sin[Subscript[q, 7][tt]]};
	
	link1 = {Thickness[0.05], LightDarkSwitched[Black, White], CapForm["Round"], Line[{{0, 0}, p1}]} ;
	link2 = {Thickness[0.05], Green, EdgeForm[LightDarkSwitched[Black, White]], CapForm["Round"], Line[{p1, p2}]} ;
	link3 = {Thickness[0.05], Blue, CapForm["Round"], Line[{p2, p3c}], Line[{{Subscript[x, b], Subscript[y, b]}, p3c}], Line[{p3c, p3}]} ;
	link4 = {Thickness[0.05], Brown, CapForm["Round"], Line[{p2, p4}]};
	link5 = {Thickness[0.05], LightRed, CapForm["Round"], Line[{p2, p6}]};
	link6 = {Thickness[0.05], Magenta, CapForm["Round"], Line[{{Subscript[x, a], Subscript[y, a]}, p5}], Line[{p5, p4}]} ;
	link7 = {Thickness[0.05], Orange, CapForm["Round"], Line[{{Subscript[x, a], Subscript[y, a]}, p7}], Line[{p7, p6}]} ;
	spring = {Thickness[0.02], Purple, CapForm["Round"], Line[{p3, {Subscript[x, c], Subscript[y, c]}}]};
	
	disk0 = {White, EdgeForm[{LightDarkSwitched[Black, White], Thick}], Disk[{0, 0}, 0.0017]};
	disk1 = {White, EdgeForm[{LightDarkSwitched[Black, White], Thick}], Disk[p1, 0.0017]};
	disk2 = {White, EdgeForm[{LightDarkSwitched[Black, White], Thick}], Disk[{Subscript[x, b], Subscript[y, b]}, 0.0017]} ;
	disk3 = {White, EdgeForm[{LightDarkSwitched[Black, White], Thick}], Disk[p2, 0.0017]} ;
	disk4 = {White, EdgeForm[{LightDarkSwitched[Black, White], Thick}], Disk[p4, 0.0017]};
	disk5 = {White, EdgeForm[{LightDarkSwitched[Black, White], Thick}], Disk[p3, 0.001]};
	disk6 = {White, EdgeForm[{LightDarkSwitched[Black, White], Thick}], Disk[p6, 0.0017]};
	disk7 = {White, EdgeForm[{LightDarkSwitched[Black, White], Thick}], Disk[{Subscript[x, a], Subscript[y, a]}, 0.0017]};
	
	Graphics[{link1, disk0, link2, link3, link4, link5, link6, link7, spring, disk1, disk2, disk3, disk4, disk6, disk7, disk5} /. andrewsParams /. solAndrews, 
	PlotRange -> {{-0.09, 0.02}, {-0.03, 0.08}}, Frame -> True, ImageSize -> 400]
	, {tt, 0.008, 0.0276}, AnimationRunning -> False, SaveDefinitions -> True]

Top

微分代数方程式の数値解法

はじめに

微分代数方程式の指数

微分方程式の取扱い

微分代数方程式の解法

不変量のある常微分方程式

質量行列を持つ系を解く

微分代数方程式の指数減少

指数減少法

構造的接続行列

Pantelides法

構造行列法

制約法

ダミー導関数

投影法

偏微分代数方程式の指数減少

微分代数方程式の矛盾のない初期化

選点法

選点の方向

選点

選点の範囲

初期条件の指定

QR法

BLT法

矛盾した初期条件

イベントを使った微分代数方程式の再初期化

微分代数方程式の例題

トーラスの問題

タンクからの圧力の解放

非伸縮性スライダークランク

接続軸にバネ制振装置を装備したスライダクランク

2リンク平面ロボットアーム

n 振り子問題

トランジスタ増幅回路

Akzo-Nobelの化学反応

車軸問題

楕円型偏微分方程式と放物型偏微分方程式の1Dでの結合

多種の食物網モデル（楕円型偏微分方程式と放物型偏微分方程式の2Dでの結合）

Andrewsの圧搾器メカニズム

関連するテクニカルノート