dasrt
ゼロ交差するDAE
呼び出し手順
[r,nn,[,hd]] = dasrt(x0,t0,t [,atol,[rtol]],res [,jac],ng, surf [,info] [,hd])
引数
- x0
y0(ydot0は 0 を初期推定値としてdasslで推定されたもの) または行列[y0 ydot0]のどちらかです.g(t,y0,ydot0)は 0 に等しい必要があります.ydot0の推定値のみが既知の場合,info(7)=1と設定してください.- y0
初期条件の実数列ベクトル.
- ydot0
t0におけるyの時間微分の 実数列ベクトル (推定値の場合もある).
- t0
初期時間を表す実数.
- t
実数のスカラーおよびベクトル. 解を求める時間を指定します.
info(2)=1を設定することにより, dasslの各ステップ点での解を得ることができることに注意してください.- nn
2個のエントリ
[times num]を有するベクトル.timesは面が交差した時刻の値,numは交差した面の数です- atol, rtol
実数スカラーまたは
yと同じ大きさの列ベクトル.atol, rtolはそれぞれ絶対および相対許容誤差を指定します. ベクトルの場合,yの各要素毎に許容誤差が指定されます.- res
外部 (関数またはリストまたは文字列).
g(t,y,ydot)の値を計算します.次のようになります :Scilab関数.
この場合, その呼出し手順を
[r,ires]=res(t,x,xdot)とする 必要があり,resは 残差r=g(t,x,xdot)とエラーフラグiresを返す必要があります.resがrの計算に 成功した場合にはires = 0,(t,x,xdot)のローカルな残差が定義されない 場合には=-1, パラメータが許容範囲外の場合は=-2となります.リスト.
外部ルーチンのこの形式は 関数にt,y,ydot以外のパラメータを指定する際に使用されます. 以下のような形式とします:
list(res,x1,x2,...)
ただし,ここで関数
resの呼び出し手順は 以下のようになりますr = res(t,y,ydot,x1,x2,...)
resはこの場合も(t,y,ydot,x1,x2,...)の関数としてr=g(t,y,ydot)を返します.警告: 関数にオプションの引数を指定しない場合, この形式を使用する必要はありません
文字列.
ScilabにリンクされるCまたはFortranサブルーチンの名前を指します.
Cの場合, 呼出し手順は以下のようにします:
void res(double *t, double y[], double yd[], double r[], int *ires, double rpar[], int ipar[])Fortranの場合, 呼出し手順は以下のようにします:
subroutine res(t,y,yd,r,ires,rpar,ipar) double precision t, y(*),yd(*),r(*),rpar(*) integer ires,ipar(*)
rparおよびipar配列は存在しますが使用できません.
- jac
外部ルーチン (関数またはリストまたは文字列). パラメータ
cjの指定した値に関してdg/dy+cj*dg/dydotの値を計算します.Scilab関数.
その呼び出し手順は
r=jac(t,y,ydot,cj)とする必要があり,jac関数はr=dg(t,y,ydot)/dy+cj*dg(t,y,ydot)/dydotを 返す必要があります. ただし,cjは実数スカラーですリスト.
以下のようにする必要があります
list(jac,x1,x2,...)
ただし,この場合の関数
jacの呼び出し手順は 以下のようになりますr = jac(t,y,ydot,cj,x1,x2,...)
この場合も
jacは(t,y,ydot,cj,x1,x2,...)の関数としてdg/dy+cj*dg/dydotを返します.文字列.
ScilabにリンクされるFortranサブルーチンの名前を指す 必要があります
Cの場合,呼び出し手順は次のようになります:
void jac(double *t, double y[], double yd[], double pd[], double *cj, double rpar[], int ipar[])Fortranの場合, 呼び出し手順は次のようになります:
subroutine jac(t,y,yd,pd,cj,rpar,ipar) double precision t, y(*),yd(*),pd(*),cj,rpar(*) integer ipar(*)
- surf
外部ルーチン (関数またはリストまたは文字列).
ng個の要素を有する列ベクトルsurf(t,y)の値を計算します. 各要素は一つの面を定義します. 以下のように定義されます:Scilab関数.
呼び出し手順は
surf(t,y)となりますリスト.
以下のようになります
list(surf,x1,x2,...)
ただし,この場合の関数
surfの 呼び出し手順は以下のようになりますr = surf(t,y,x1,x2,...)
文字列.
ScilabにリンクされるCまたはFortranサブルーチンの名前を指します
Cの場合,呼び出し手順は次のようになります:
void surf(int *ny, double *t, double y[], int *ng, double gout[])
Fortranの場合,呼び出し手順は次のようになります:
subroutine surf(ny,t,y,ng,gout) double precision t, y(*),gout(*) integer ny,ng
- info
7個の要素を含むオプションのリスト. デフォルト値は list([],0,[],[],[],0,0);- info(1)
gを評価できる時間の最大値 を指定する実数スカラー. または,時間に制限を課さない場合には,空の行列[]- info(2)
dasrtが その中間的な計算値を返す(flag=1), もしくはユーザが指定した時間においてのみ値を返す (flag=0) かどうかを指定するフラグ.- info(3)
2個の要素のベクトルで,jac;r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)"により計算された 帯行列の定義[ml,mu]を指定します.jacが完全な行列を返す場合,info(3)=[]を指定します.- info(4)
実数スカラーで, ステップの大きさの最大値を指定します. 制限をしない場合には,
info(4)=[]を指定してください.- info(5)
実数スカラーで,ステップの大きさの初期値を指定します. 指定しない場合には,
info(5)=[]を指定してください.- info(6)
解が非負であることが既知の場合に
info(6)=1に設定, そうでない場合はinfo(6)=0に設定します.- info(7)
ydot0が単なる推定値の場合はinfo(7)=1,g(t0,y0,ydot0)=0の場合はinfo(7)=0を指定します.
- hd
実数ベクトルで,
dasrtのコンテキストを保存し, 積分を再開することを可能にします.- r
実数行列. 各列はベクトル
[t;x(t);xdot(t)]です. ただしtは解を計算した時間です.
説明
陰的な微分方程式を解きます.
g(t,y,ydot) = 0 y(t0) = y0 and ydot(t0) = ydot0
面を横切る時間および
到達した面の数をnnに返します.
詳細な例は SCIDIR/tests/unit_tests/dassldasrt.tst で見つけることができます
例
//dy/dt = ((2*log(y)+8)/t -5)*y, y(1) = 1, 1<=t<=6 //g1 = ((2*log(y)+8)/t - 5)*y //g2 = log(y) - 2.2491 y0 = 1; t = 2:6; t0 = 1; y0d = 3; atol = 1.d-6; rtol = 0; ng = 2; deff('[delta, ires] = res1(t,y,ydot)', 'ires=0; delta=ydot-((2*log(y)+8)/t-5)*y') deff('rts = gr1(t,y)', 'rts = [((2*log(y)+8)/t-5)*y;log(y)-2.2491]') [yy, nn] = dasrt([y0,y0d], t0, t, atol, rtol, res1, ng, gr1); //(Should return nn=[2.4698972 2])
参照
| Report an issue | ||
| << daskr | Differential Equations, Integration | dassl >> |