Maple で常微分方程式 (ODE)の数値解を扱う際の効率に関しては、様々な考察があります。このページは、解の値やプロットの計算時の dsolve/numeric の効率的な使用、および、すべての dsolve/numeric プロシージャ形式の解についての全般的な考察に関する情報の中央保管場所としての役割を果たします。

環境変数 Digits の設定は、dsolve/numeric の効率に極めて重要です。Digits が evalhf(Digits) の値(現在、すべての対応プラットフォームで 15)以下に設定されている場合、解の計算はハードウェアの浮動小数点演算を使用して実行し、53 binaryの精度の Digits を使用ます。

Digits が15以下の場合、計算がハードウェアモードで実行されていると言い、可能な場合には、ハードウェア浮動小数で高速な評価を行うevalhf が使用されます。こうすると効率が大幅に促進され、場合によっては、 程度の因数によって計算を加速することができます。

さらに、これをサポートする方法としては、ハードウェアモードでの演算中には compile オプションを使用することができます。このオプションにより、 に加えて  の因数により、計算を加速することができます。

ハードウェアモードでの演算中には、値が 15 以下である限り Digits の設定は関与しません。従って、 で実行される計算では、たとえば  や  で実行した計算と同一の結果が出ます。

Digits := 16:

dsn16 := dsolve({diff(y(t),t,t)+y(t),y(0)=0,D(y)(0)=1}, numeric, abserr=1e-12, relerr=1e-12, maxfun=0);

tt := time():

dsn16(50*Pi);

t_sw := time()-tt;

Digits := 10:

dsn10 := dsolve({diff(y(t),t,t)+y(t),y(0)=0,D(y)(0)=1}, numeric, abserr=1e-12, relerr=1e-12, maxfun=0);

tt := time():

dsn10(50*Pi);

t_hw := time()-tt;

t_sw/t_hw;

dsn10 := dsolve({diff(y(t),t,t)+y(t),y(0)=0,D(y)(0)=1}, numeric, abserr=1e-12, relerr=1e-12, maxfun=0, compile=true);

tt := time():

dsn10(50*Pi);

t_hc := time()-tt;

t_sw/t_hc, t_hw/t_hc;

dsolve/numeric で出力される大部分のプロシージャは、compute-as-needed となります。すなわち、その系は解の値が要求されたときに、初期点または最後に計算された解の点のいずれから積分されます。

ここで「大部分」と表現したのは、特定の問題またはオプションは 解の保存を強制したり、暗示したりするためです。これには、BVP 問題の解、 解を事前計算するための range オプションを使用した IVP の解、または、piecewise 出力形式を使用した解が含まれます。

compute-as-needed 手法で深刻な問題となるのは、IVP 問題の安定性が明らかでないため、積分の方向が反転できないという点でです。

たとえば、IVP 問題の初期点が  であり、 で値が要求され、その後、次の要求が  で行なわれる場合、IVP 系は、両要求に対して  からすべてに積分を実行する必要があります。逆に、 での値が最初に要求され、次に  での値が要求される場合、2 番目の要求では方向の反転が不要であるため、2 番目の要求については、積分は  から実行すれば良いのです。 

(とりわけ)こうした理由により、特殊な IVP の解をプロットするルーチンである plots[odeplot] が準備されており、これを使用して ODE の解のプロットを効率的に作成します。plot のような標準のプロットルーチンを使用すると、処理される点の次数に対する制限に関して把握する手段がないため、非効率になる場合があります。さらに、特殊なルーチンを使用すると、ODE 積分用の特殊出力モデルが可能となり、解法とプロット作成との間の強化された(さらには効率化された)結び付きサポートします。

dsn := dsolve({diff(y(t),t,t)+y(t),y(0)=0,D(y)(0)=1}, numeric, abserr=1e-12, relerr=1e-12, maxfun=0, output=listprocedure);

tt := time():

plots[odeplot](dsn,[t,y(t)],0..50*Pi,numpoints=500);

t_spec := time()-tt;

tt := time():

plot(op([2,2],dsn),0..50*Pi);

t_gen := time()-tt;

t_gen/t_spec;

IVP 問題では、入力系は dsolve/numeric/IVP で説明されているように、既に解かれた 1 階形式に変換する必要があります。一部の系においては、この記号解のステップは非常に大きな負担をかける場合があります。これは特に、自動生成される濃密系やホモトピー問題からの系で発生します。

幸いにも導関数において系が線形の場合 (1 階に変換後)には、implicit オプションを使用することができます。implicit オプションは記号の backsolve (解のプロシージャの作成時に 1 度実行される)を遅延させ、純粋な数値の backsolves (積分中の導関数の各評価において実行される)にします。

記号式の膨らみや、大きな式であることから発生する関連丸め誤差の発生の防止に加え、この方法には、未知の記号係数に代わって系の係数の局所数値でピボットが行われるため、削除プロセスにおいて変数を解く際により安定した選択が行われるという利点があります。

eq1 := cos(t)*diff(x(t),t)+(1/2-sin(t))*diff(y(t),t)=1/2-y(t);

eq2 := sin(t)*diff(x(t),t)+cos(t)*diff(y(t),t)=x(t);

dsn := dsolve({eq1,eq2,x(0)=0,y(0)=0}, numeric, implicit=true);

plots[odeplot](dsn,[x(t),y(t)],0..6*Pi);

別の dsolve/numeric の解に依存する ODE 問題の dsolve/numeric の解を計算する必要があることがあります。たとえば、dsolve/numeric 問題から  の解を得るとします。その解は次に、強制関数といった別の問題の一部として使用されるとします。

時間の尺度が同等な場合、2 つの系を単純に組み合わせることにより、遥かに効率よく解を求めることができます。これは、dsolve/numeric の解のプロシージャに依存する系を evalhf (第 1 セクション参照)のもとで積分することができないために、積分の速度が遥かに遅くなるためです。

sys1 := {diff(x(t),t,t)+x(t)=0, x(0)=0, D(x)(0)=1};

sol1 := dsolve(sys1, numeric, output=listprocedure);

xf := op([2,2],sol1);

sys2 := {diff(y(t),t,t)=xf(t), y(0)=0, D(y)(0)=1};

sol2 := dsolve(sys2, numeric, known={xf});

tt := time():

sol2(100);

t_separate := time()-tt;

sysc := sys1 union eval(sys2,xf(t)=x(t));

solc := dsolve(sysc, numeric);

tt := time():

solc(100);

t_combined := time()-tt;

t_separate/t_combined;

時間の尺度が異なる場合 (たとえば、第 2 の系が  に依存し、このとき  が別の ODE 問題の解である場合)、デフォルトの stiff および non-stiff の方法では、所定の範囲では区分形式の出力が可能です。区分解の計算は、第 2 の問題の望ましい積分の制限を十分に超えて実行することが必要となる可能性があることに注意してください。これは、ソルバが積分中に現在の範囲を超えて 進むためです。

sol1p := dsolve(sys1, numeric, output=piecewise, range=0..110):

xp := op([2,2],sol1p):

type(xp,specfunc(anything,piecewise)), nops(xp);

sys2p := {diff(y(t),t,t)=xf, y(0)=0, D(y)(0)=1}:

sol2p := dsolve(sys2p, numeric);

tt := time():

sol2p(100);

t_piecewise := time()-tt;

t_separate/t_piecewise;

最後に、必要であれば、(区分的でない)ストレートな range の解を使用することができます。この方法では、少なくとも解の再計算をせずに済みますが、evalhf モードでは実行されないため、選択肢の中では最も非効率な方法になります。