dsolve/numeric/IVP - ODE 初期値問題の数値解を求める
使い方
dsolve(odesys, numeric, vars, options)
dsolve(procopts, numeric, vars, options)
パラメータ
odesys - 集合またはリスト; (連立)常微分方程式および初期条件
numeric - 名前; dsolve に数値解を求めるように指示
vars - (オプション) odesys に関する従属変数、あるいは従属変数の集合またはリスト
procopts - (odesys が設定されていない場合に必要) 系を定義する手続きを指定するためのオプション
options - (オプション) keyword = value の形をした等式
|
説明
|
|
•
|
numeric または type=numeric のオプションと初期値問題(IVP)を与えて dsolve コマンドを呼び出すと、dsolve は ODE または連立 ODE の IVP に関する数値解を求めます。オプションの等式 method=numericmethod (ここで numericmethod は rkf45, rosenbrock, dverk78, lsode, gear, taylorseries, classical のいずれか)が与えられると、dsolve は数値解を得るためにその手法を使用します。
|
'output' = keyword or array
'range' = numeric..numeric
'abserr' = numeric
'relerr' = numeric
'stiff' = boolean
'known' = name or list of names
•
|
一階の変換では、上記にある abserr のためのインプット形式のリストは例外となります。
|
|
最もな場合、全ての解の計算のための abserr に対する値は一つで十分です。この例外は、方程式系の他の変数より差のスケールの大きい値を持つ変数の導関数のためのいくつかの従変数がある場合となります。この場合、abserr のためのリストの入力形式は方程式系のそれぞれの従変数の絶対誤差許容量の差、あるいは方程式系の解の要素を用いて指定することが可能です。これらの差は、方程式系の従変数の許容の指定が、常に各変数の解の計算が独立と考えられている間の変数の導関数のためのものです。
|
|
abserr に依存、あるいは要素となるように指定するため、出力内の従変数のオーダーは既知である必要があります。これは従変数 (dsolve の vars パラメータ )を指定することで、リストとしてコントロールすることが可能です。
|
|
これらのことから従変数 abserr は、指定した変数と指定した abserr に 1-1 対応の関係があることが容易に理解できます。単位要素の abserr は、基本的な考え方が同じであるため少し扱いにくい部分もありますが、各変数のオーダーは一つの abserr の値よりも依存している必要かもしれません。
|
|
例えば、 x'',y''',z' を含む問題について、もし変数が [x(t),y(t),z(t)] として指定されたならば、単位要素 abserr のために 6 つの値が必要となり、 [x(t),x'(t),y(t),y'(t),y''(t),z(t)] に対応していることになります。
|
•
|
デフォルトの (non-stiff と stiff 問題に対する ) IVP による計算法のみが、abserr に対するリスト形式をサポートしています。
|
•
|
taylorseries を除く全ての IVP 法は、 実数の独立変数を持つ複素数値の IVP を扱う機能を持っています。
|
•
|
以下のオプションは IVP で指定できるもので、主に手続きの形での入力 (procopts) 系の指定
|
'number' = integer
'procedure' = procedure
'start' = numeric
'initial' = array
'procvars' = list
|
さらに、解法の制御や刻み幅の制御に関連するものです。 :
|
'startinit' = boolean
'minstep' = numeric
'maxstep' = numeric
'initstep' = numeric
'stop_cond' = list
'implicit' = boolean
•
|
手続きに基づくオプションを議論する前に、連立 ODE の IVP を数値的に扱う Maple の処理について説明しておきます。最初に、全ての連立微分方程式は、新たな従属変数を導入して1階の連立微分方程式に変換されます。例えば、単独の ODE
|
3 y''' = 2 y' - x y = 0
|
を数値的に解くには、新たな従属変数 u=y', v=y'' が導入され、単独の3階 ODE は、3つの従属変数を含んだ3つの1階連立 ODE になります。さらに、連立1階微分方程式の導関数は分離しています。そうすると
|
y' = u
u' = v
v' = ( 2 u - x y ) / 3
|
この時点で、ODE は1階の連立方程式に変換されているので、手続きに基づいたオプションをより厳密に記述することが可能です。これ以降はすべて、連立方程式や従属変数についての記述は、1階の連立方程式に関するものに限ります。
|
|
オプション
|
|
|
従属変数の導関数を評価する手続きを指定する整数です。これは proc(N,x,Y,YP) という形でなくてはなりません。ここで N は従属変数の個数、 x は導関数を評価する独立変数の値、Y は x での従属変数の値の配列、YP は手続きが導関数の計算値を入れるための出力配列を表します。
|
|
上の例では、変数の対応は Y[1]=y, Y[2]=u および Y[3]=v で与えられ、以下の手続きによって 連立 ODEの導関数が正しく評価されます。
|
solproc := proc(N, x, Y, YP)
YP[1] := Y[2];
YP[2] := Y[3];
YP[3] := (2*Y[2] - x*Y[1])/3;
end;
|
1階の連立方程式における従属変数の初期値を指定する配列です。この配列の値は、引数 procedure で与えられる手続きで使用されるものと同じ次数で表されなくてはなりません。
|
|
解の procedure における添え字付きの値に対応する従属変数の名前と、initial 値の配列を指定するリストです。上の例では、 procvars=[y(x),u(x),v(x)] が使われます。もとの ODE の変数で解を表すには、procvars=[y(x),diff(y(x),x),diff(y(x),x,x)] が使われます。
|
•
|
単独の ODE あるいは連立の ODE が直接入力された場合(連立問題では 手続きのオプションを使用しない)、 新しい従属変数として導関数に変数名を付け直すことで、Maple は内部的に連立方程式を1階に変換します。この時、number, procedure, start, initial, procvars の値は、入力から得ます。 連立方程式と手続きの形の入力が両方与えられた場合は、連立方程式の形は無視され、警告が発せられます。
|
|
注意: 連立の ODE が主導関数(連立方程式の各従属変数で最も高次の導関数)について解いた形として与えられる場合、連立方程式の代数的取り扱いは、(新しい従属変数として導関数の変数名を付け直すのではなく)、手続きの形を得る中で行われます。これにより、連立方程式内にある導関数の実際の計算をよりよく制御することが可能になります。連立 ODE の方程式のいずれかが解いた形でないと、解いた形を得るために solve が呼び出されます。これは、ある種の好ましくない副次効果を起こすことがあります(例えば、solve(y-(x+1)^20-1,y) を考えてみて下さい)。
|
|
主導関数について連立の ODE が線形である場合にのみ、true の値をとる引数です。またこの引数は、連立方程式により定義された問題で次の手法と共に使用されます: rkf45, dverk78, gear, lsode, classical 。デフォルトでは、このオプションは false となっています。しかし、true に設定することにより、 dsolve[numeric] 方程式を解かない形のままにして、導関数値が必要となる度に1次方程式を解いて、導関数値の計算を行います。これは、 n 個の関数の評価をステップ毎に必要とする手法が、ステップ毎に n 回1次方程式を解くことが必要であることを意味します。
|
|
このオプションを使用することで、ほとんどの場合では計算が遅くなりますが、主導関数の分離により式が膨大になり過ぎたり、導関数の値を計算する際に不良な条件を持つ関数が現れるような複雑な係数を持つ連立 ODE に対しては、必要なことです。この挙動は、多項式系の解を微分方程式系の解にホモトピー法を使用してマッピングする際にしばしば起こります。
|
|
連続した計算に関する数値積分の挙動を制御するブール値です。startinit=true が指定されると、それぞれの計算は問題の初期点から始められます。 startinit=false (デフォルトの設定)では、積分の方向を逆転しないならば、前の計算が終わった点からそれぞれの計算が始められます。積分方向を逆転する場合には、計算は初期点から始められます。
|
|
デフォルトの解法 rkf45 と rosenbrock では、初期点と前に解を求められた他の点との間に新しい独立変数の値がない場合にのみ、解が再計算されます。これは積分方向を決して逆転させない効果があり、また既に計算した区間での解の評価に全く負担をかけません(解の値は補間によって与えられます)。パラメータ startinit はまた、これらの解法に対して、解の値が求められる時に再計算するよう、これらの解法に伝えます。
|
•
|
与えられた入力された問題の解の値の計算に加え、手続きの出力 (つまり、output=procedurelist, listprocedure, operator) は、追加の会話形式で提供します(例の最後の箇所を参照してください)。
|
|
手続き dsn に対して、dsn('initial')と使用すると、計算に使用した初期値を返します。
|
|
最後に計算された点での解の値を返します (たとえば、特異値問題など、ある点で積分できない問題に有用です)。
|
|
問題に対する新しい初期データとして、手続きが value 内で値を使用するようにします。procedurelist の出力に対して、value は、n あるいは n+1 の数値のリストになり、 従属変数値と導関数、あるいは、独立変数値と導関数にそれぞれ対応する必要があります。これらは、出力に与えられるものと同じ次数である必要があります。 listprocedure と operator の出力に対して、value は、 procedurelistと同じか、あるいは、呼び出された手続きに相当する変数の初期値に対応する1つの数値になります。注意: 手続き listprocedure あるいは operator のいずれか1つの変更は、dsolve の呼び出しで返される他のすべての手続きに影響します。それらは、独立ではありません。
|
|
|
|
例
|
|
y(x)=sin(x) の初期値問題: 両方向から積分します:
>
|
dsys1 := {diff(y(x),x,x)+y(x)=0, y(1)=sin(1), D(y)(1)=cos(1)};
|
| (2.1) |
>
|
dsol1 := dsolve(dsys1, numeric, method=gear);
|
| (2.2) |
| (2.3) |
| (2.4) |
連立 IVP :
>
|
deq2 := {diff(x(t),t)=y(t),diff(y(t),t)=x(t)+y(t)};
|
| (2.5) |
>
|
ic2 := {x(0)=2,y(0)=1}:
dsol2 := dsolve(deq2 union ic2, numeric, output=listprocedure,
range=0..1);
|
| (2.6) |
>
|
dsol2x := subs(dsol2,x(t)): dsol2y := subs(dsol2,y(t)):
dsol2x(0), dsol2y(0);
|
| (2.7) |
>
|
dsol2x(0.4), dsol2y(0.4);
|
| (2.8) |
>
|
dsol2x(1.0), dsol2y(1.0);
|
| (2.9) |
手続きの指定を用いた Lorenz 系:
>
|
dproc3 := proc(N,t,Y,YP)
YP[1] := sigma * (Y[2] - Y[1]);
YP[2] := r*Y[1] - Y[1]*Y[3] - Y[2];
YP[3] := Y[1]*Y[2] - b*Y[3];
end:
sigma := 10: r := 28: b := 8/3:
ic3 := array([5,5,5]):
dvars3 := [x(t),y(t),z(t)]:
dsol3 := dsolve(numeric, number=3, procedure=dproc3,
start=0, initial=ic3, procvars=dvars3);
|
| (2.10) |
| (2.11) |
| (2.12) |
| (2.13) |
次の例は、( (exp(x)*y')'+sin(x)*y(x)=0 の形の)問題で、ここでは変数の自然な名前付けは望むものではありません。1番目の従属変数を y(x), 2番目を v(x)=exp(x)*y'(x) として選択し、問題の手続き指定を使用します。
>
|
dproc4 := proc(N,x,Y,YP)
YP[1] := exp(-x)*Y[2];
YP[2] := -sin(x)*Y[1];
end;
|
| (2.14) |
>
|
ic4 := array([1,0]):
dvars4 := [y(x),exp(x)*diff(y(x),x)]:
dsol4 := dsolve(numeric, number=2, procedure=dproc4,
start=0, initial=ic4, procvars=dvars4);
|
| (2.15) |
| (2.16) |
| (2.17) |
| (2.18) |
複素数値の初期値問題:
>
|
dsys1 := {diff(y(x),x)+I*y(x)=0,y(0)=1+2*I};
|
| (2.19) |
>
|
dsol5 := dsolve(dsys1, numeric, method=gear);
|
| (2.20) |
| (2.21) |
| (2.22) |
| (2.23) |
会話形式の例:
>
|
dsys := {diff(y(x),x)=y(x)^2, y(0)=1};
|
| (2.24) |
>
|
dsoli := dsolve(dsys, numeric);
|
| (2.25) |
計算を実行します。
| (2.26) |
初期値、最後に計算した値、計算方法をたずねます。
| (2.27) |
| (2.28) |
| (2.29) |
特異点に至るまで計算します。
Error, (in dsoli) cannot evaluate the solution further right
of .99999999, probably a singularity
| |
特異点に最も近い、この方法で値を得ることができる点をたずねます。
| (2.30) |
初期条件を変更して、再び計算します。
>
|
dsoli('initial'=[1,1/2]);
|
| (2.31) |
Error, (in dsoli) cannot evaluate the solution further right
of 2.9999999, probably a singularity
| |
| (2.32) |
listprocedure の出力:
>
|
dsoli2 := dsolve(dsys, numeric, output=listprocedure);
|
| (2.33) |
| (2.34) |
| (2.35) |
計算と質問
| (2.36) |
| (2.37) |
初期の x 点のみ変更し、y(1.5) を計算します。
| (2.38) |
| (2.39) |
両方を変更します。
>
|
dsy('initial'=[1,1/2]);
|
| (2.40) |
Error, (in dsy) cannot evaluate the solution further right
of 2.9999999, probably a singularity
| |
>
|
dsx('last'),dsy('last');
|
| (2.41) |
|
|
参照
|
|
dsolve[classical], dsolve[dverk78], dsolve[Error_Control], dsolve[gear], dsolve[lsode], dsolve[maxfun], dsolve[numeric], dsolve[numeric,IVP], dsolve[rkf45], dsolve[rosenbrock], dsolve[Stiffness], dsolve[Stop_Conditions], dsolve[taylorseries], plots[odeplot]
|
|