dasrt

パラメータ

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])

参照

• ode — 常微分方程式ソルバ
• dassl — 微分代数方程式
• impl — 微分代数方程式
• fort — Fortran or C user routines call
• link — dynamic linker
• external — Scilabオブジェクト, 外部関数またはルーチン