オプションに numeric および method=rkf45 を指定した dsolve コマンドは、Taylor 級数法 (Taylor series method) を用いて数値解を求めます。この計算法は高精度な解を得るために使用することができます。

注意 : この計算法で低い精度の解を得る場合は、常に他の計算法に比べてより多くの時間を使います。したがって、この計算法は高精度な解を求めたい場合にのみ使用することを薦めます。

この解法では副解法について二つの選択をします。lazyseries choice を指定すると、dsolve は全ての積分の前に雑な級数展開を用いた Taylor 級数係数を生成させるルーチンを構成します。これはデフォルトであり、一般に最も効果的な計算法となります。series choice を指定すると、各積分ステップ dsolve,series は局所的な級数展開の計算を呼び出します。

taylorseries による計算法では以下の options が使用可能となります。 :

'output'

dsolve の出力形式を指定します。

'known'

ユーザ定義の既知関数を指定します。基本的な用法に関しては dsolve[numeric] をご参照ください。ここに書かれているように、numeric を指定した評価関数で定義された手続きに用件を加えます。この計算法では、常に `diff/` の規則を必要とします。

`diff/` の規則は独立変数によって表された関数の全微分を与えます。注意 : ここでのキーワードは全微分 (total) であり、その例は以下で見ることができます。これは rosenbrock による計算法のもので、diff の規則は偏微分を与えなければなりません。

任意の全次数が必要となり、これは既知 (known) 関数の導関数が、別のユーザ定義の関数で与えられていることを意味します。このとき、独立変数によって表された全微分の計算には、関数に `diff/` の規則がなければなりません。再び、この新しいユーザ定義の関数は、別のユーザ定義の関数を用いて全微分を与えます。そしてさらに再び、独立変数によって表された全微分の計算には、関数にうんざりするほどの `diff/` の規則がなければなりません。

最後の要件として、'known' オプションはデフォルトの taylorseries に関して choice だけを使用することができます。これはすなわち lazyseries を使用することになります。

'startinit'

あらゆる要求を含む解の計算のために、以前に計算された結果から計算を連続して続けるべきか、あるいは新たに初期条件から計算するかを表示するように指定します。

注意 : 続ける場合は、初期条件に比例して、積分の支持は決して逆にされることはありません。

'abserr'

個々のステップのための絶対誤差を制御する正の浮動小数点数を指定します。デフォルトの値は Float(1,-Digits) となっています。

'minstep'

最小のステップ幅を与える浮動小数点数を指定します。デフォルトの値は 0.00001 となっています。

'order'

積分で扱う Taylor 級数展開のオーダー指定します。デフォルトの値は max(22,trunc(3*Digits/2)) で与えられています。

'range'

前計算によって記憶される解の範囲を指定します。このオプションを指定した場合、全ての返された解の値は記憶された Taylor 級数展開の一つの形式 - 計算法の各ステップの記憶された一つ - で計算されます。このオプションの使用は storage を on にする効果があります。

'storage'

全てのステップで解を記憶するための指定となります。これは、既に解が計算されて記憶された領域では Taylor 多項式を使用します。デフォルトでは false となっています。

'steppast'

解を積分する場合に、自然なステップ幅 ( つまり、誤差の制御において禁止されるようなステップ幅  ) を、独立変数の値が (内部的に) 要求された点の先に置かれている場合でも使用するための指定となります。この解は、要求された点を含んでいるような領域にわたって Taylor 多項式の形式で計算されます。これは密 (dense) な主力にはより効果的であり、要求された点の先の解を続けて計算することができます。デフォルトでは true となっています。

結果は plots パッケージの odeplot 関数を使用してプロットすることができます。

Digits := 20:
dsys1 := {diff(x(t),t)=y(t), diff(y(t),t)=x(t)+y(t),
          x(0)=2, y(0)=1};

dsol1 := dsolve(dsys1, numeric, method=taylorseries,
                output=listprocedure, abserr=Float(1,-20)):
dsol1x := eval(x(t),dsol1): dsol1y := eval(y(t),dsol1):
dsol1x(0), dsol1y(0);

dsol1x(0.4), dsol1y(0.4);

dsol1x(1.0), dsol1y(1.0);

dsola1:= dsolve(dsys1, numeric, method=taylorseries,
                output=array([0,0.2,0.4,0.6,0.8,1]),
                abserr=Float(1,-20));

dsys2 := {(D@@2)(x)(t)=-y(t), (D@@2)(y)(t)=D(x)(t)+y(t)};

init2 := {x(0)=1, D(x)(0)=0, y(0)=0, D(y)(0)=1};

dsol2 := dsolve(dsys2 union init2, numeric, method=taylorseries,
                output=array([0,.6,1.1,1.5,2.3,2.5]),
                abserr=Float(1,-20) ):
dsol2[1,1][1]=dsol2[2,1][6,1],dsol2[1,1][5]=dsol2[2,1][6,5];

Digits := 30:
deqn3 := { diff(y(t),t$3) - 2*diff(y(t),t$2) + 2*y(t) };

init3 := { y(0)=1, D(y)(0)=1, (D@@2)(y)(0)=1 };

dsol3 := dsolve(deqn3 union init3, numeric, method=taylorseries,
                abserr=Float(1,-30)):
dsol3(0.5);

Digits := 20:
dsys4 := {diff(y(t),t,t)+y(t), y(0)=0, D(y)(0)=1};

dsol4 := dsolve(dsys4, numeric, method=taylorseries,
                abserr=Float(1,-20), range=0..10*Pi);

fsolve(t->rhs(dsol4(t)[2]),3..4);

%-evalf(Pi);

fsolve(t->rhs(dsol4(t)[2]),6..7);

fsolve(t->rhs(dsol4(t)[2]),9..10);

fsolve(t->rhs(dsol4(t)[2]),25..30);

dsol5a := dsolve(dsys4, numeric, method=taylorseries,
                 abserr=Float(1,-20));

tt := time():
for i to 100 do dsol5a(Pi*i/100): end do:
time()-tt;

dsol5b := dsolve(dsys4, numeric, method=taylorseries,
                 abserr=Float(1,-20), steppast=false);

tt := time():
for i to 100 do dsol5b(Pi*i/100): end do:
time()-tt;

dsys6 := {(t+1)^2*diff(y(t),t,t) + (t+1)*diff(y(t),t)
          + ((t+1)^2-0.25)*y(t) = 0,
          y(0) = 0.6713967071418030,
          D(y)(0) = 0.09540051444747446};

ds1 := dsolve(dsys6, numeric, method=taylorseries['series'],
          abserr=Float(1,-20));

tt := time():
ds1(3.0);

time()-tt;

ds1 := dsolve(dsys6, numeric, method=taylorseries[lazyseries],
          abserr=Float(1,-20));

tt := time():
ds1(3.0);

time()-tt;

f(x,y(x)) = Int(exp(-t*2/100)/10,t=0..x)-y(x);

f := proc(x,y) local t;
   if not type(evalf(x),'numeric') or
      not type(evalf(y),'numeric') then
      'procname'(x,y);
   else
      evalf(Int(exp(-t*2/100)/10,t=0..x))-y;
   end if;
end proc;

`diff/f` := proc(x,y);
   if args[3]=x then
      diff(x,args[3])*exp(-x*2/100)/10-f(x,y);
   elif args[3]=y then
      error "should not happen";
   else
      error "unable to differentiate";
   end if;
end proc;

dsol := dsolve({diff(y(x),x)=f(x,y(x)), y(0)=0}, numeric,
                known=f,method=taylorseries);

dsol(10);

Barton, D.; Willers, I.M.; and Zahar, R.V.M. Mathematical Software, pp. 369-389. Edited by J.R. Rice. New York: Academic Press, 1972 .