dsolve の DAE extension method (DAE 拡張法) は、DAE 問題を解く目的で拡張された ODE IVP の標準的な数値法で構成されています。

具体的には、標準的な ODE および (希望の従属変数に対して分離可能な) 指数が 1 に制限された DAE が解けるようにソルバが拡張されています。方程式システムは、index reduction (指数減少法) によって、システムの指数が 1 になるまで簡約されます。

注意：指数は、DAE システムがどれだけ ODE システムから離れているかを示す指標として考えることができます。

DAE 拡張法では、各ステップで改良されたソルバが適用され、各ステップの最後には、問題に対する追加制約が満たされているかを確認するために、解が逆投影されます。これらの追加制約は、指数減少プロセスの中で現れるもので、指数減少を必要とする方程式システムにのみ存在します。

dsolve には、rkf45 法の拡張である rkf45_dae、ck45 法の拡張である ck45_dae、および  rosenbrock 法の拡張である rosenbrock_dae の 3 つの拡張法があります。

拡張であるこれらの数値法は、たとえば、イベントや範囲解など元のソルバの機能の多くを持ち合わせています。そのため、このページでは主に両者の相違点および拡張法固有のオプションについてのみ説明します。

まず異なる点は、関連後処理後のシステムは通常かなり複雑であるため、拡張法では optimization がデフォルトで有効になっている点です。

問題を解くには、初期条件に一貫性が求められます。つまり、初期条件は問題に存在する隠れた制約事項をすべて満たしている必要があります。制約事項が十分に満たされていない場合は、エラーが発生し、満たされていない条件の情報が表示されますので、それらを調整します。

さらに、differential が false (偽) の場合は、問題で指数が 1 であるいずれの変数に対しても自動的に数値計算を行うことができます。これらの変数については、初期条件を指定する必要がありません。

機能面では、元のソルバと異なり、拡張ソルバは実数値問題にのみ使用可能であるという点と、問題を解くための初期値が対話的に変更できないという点の 2 つの大きな違いがあります。

DAE 拡張ソルバに固有のオプションは 3 つです：

'projection'= boolean

このオプションは、各ステップの後に、ソルバで直接使用されなかった制約事項に対して解を逆投影するかどうかを指定します。このオプションのデフォルト値は true です。DAE IVP の安定性が高く、速度が焦点の場合は、このオプションを false に設定することで、より短時間で妥当な解を導き出すことができます。

'differential'= boolean

このオプションは、システムが完全に微分可能になるまで指数減少を継続するよう (つまり、改良型ソルバの指数 1 機能の不適用を) 指定します。このオプションのデフォルト値は false です。ソルバの指数 1 機能を使用すると、大抵指数が完全に簡約されたシステムより安定性が高く、硬さ (stiffness) がより小さく、通常は、より早く解けるシステムを導き出すことができます。

'implicit'= boolean

このオプションは、評価が実行されるまで拡張システムを部分的に未解決 (暗黙) のまま残してもかまわないことを指定します。DAE 問題の数値解を得るために解決する必要があるシステムの部分と暗黙的に残せるシステムの部分はソルバによって判断されます。このオプションは、基本的には IVP 問題 (dsolve[numeric,IVP] を参照)  で使用される implicit オプションの高機能版で、有用なオプションです。数値積分を計算するためのシステム導量をより単純なものにしますが、式の評価には線形解決しか使用されません。そのため、主な導関数に対してほぼ解決されているシステムでは、問題の解決に直接法より若干時間がかかることになります。主な導関数に対して難易度が高い問題については、早い解決が期待できます。

他の留意事項については、dsolve[numeric,DAE] のページでも説明されています。

dsys := {
   diff(x(t),t,t) = -2*lambda(t)*x(t),
   diff(y(t),t,t) = -2*lambda(t)*y(t)-Pi^2,
   x(t)^2+y(t)^2 = 1,
   x(0)=0, D(x)(0)=1/10, y(0)=-1, D(y)(0)=0,
   lambda(0) = 1/200+Pi^2/2
};

dsol1 := dsolve(dsys, numeric, method=rkf45_dae, maxfun=0);

t1 := time():

r1 := dsol1(2000);

t1 := time()-t1;

dsol2 := dsolve(dsys, numeric, method=rkf45_dae,
                projection=false, maxfun=0);

t2 := time():

r2 := dsol2(2000);

t2 := time()-t2:

map(rhs,r1)-map(rhs,r2);

t1,t2;

dsol3 := dsolve(dsys, numeric, method=rkf45_dae,
                differential=true, maxfun=0);

t3 := time():

r3 := dsol3(2000);

t3 := time()-t3:

map(rhs,r1)-map(rhs,r3);

t1,t3;