dassl

引数

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の各ステップ点での解を得ることができることに注意してください.

atol, rtol

実数スカラーまたはyと同じ大きさの列ベクトル. atol, rtol はそれぞれ絶対および相対許容誤差を指定します. ベクトルの場合,yの各要素毎に許容誤差が指定されます.

res

外部 (関数またはリストまたは文字列). g(t,y,ydot)の値を計算. 以下のようにします :

• Scilab関数.

  呼び出し手順は [r,ires]=res(t,y,ydot) とする必要があり, res は残差 r=g(t,y,ydot) と誤差フラグiresを 返す必要があります. resがrの計算に 成功した場合にはires = 0, (t,x,xdot)のローカルな残差が定義されない 場合には =-1, パラメータが許容範囲外の場合は =-2 となります.

• リスト.

  外部ルーチンのこの形式は 関数にパラメータを指定する際に使用されます. 以下のような形式とします:

  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(*)

info

7 個の要素を含むオプションのリスト. デフォルト値は list([],0,[],[],[],0,0);

info(1)

gを評価できる時間の最大値 を指定する実数スカラー. または,時間に制限を課さない場合には,空の行列 []

info(2)

dassl が その中間的な計算値を返す(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

実数ベクトルで,dasslのコンテキストを保存し, 積分を再開することを可能にします.

r

実数行列. 各列はベクトル[t;x(t);xdot(t)]です. ただし t は解を計算した時間です.

説明

dassl関数は代数微分方程式を積分し, 指定した時間における y の 値を返します.

g(t,y,ydot)=0
y(t0)=y0  and   ydot(t0)=ydot0

例

function [r, ires]=chemres(t, y, yd)
   r=[-0.04*y(1)+1d4*y(2)*y(3)-yd(1)
       0.04*y(1)-1d4*y(2)*y(3)-3d7*y(2)*y(2)-yd(2)
       y(1)+y(2)+y(3)-1];
   ires=0
endfunction
function pd=chemjac(x, y, yd, cj)
    pd=[-0.04-cj , 1d4*y(3)               , 1d4*y(2);
         0.04    ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);
         1       , 1                      , 1       ]
endfunction
y0=[1;0;0];
yd0=[-0.04;0.04;0];
t=[1.d-5:0.02:.4,0.41:.1:4,40,400,4000,40000,4d5,4d6,4d7,4d8,4d9,4d10];
y=dassl([y0,yd0],0,t,chemres);
info=list([],0,[],[],[],0,0);
info(2)=1;
y1=dassl([y0,yd0],0,4d10,chemres,info);
y2=dassl([y0,yd0],0,4d10,chemres,chemjac,info);
//パラメータとしてオプション引数を使用
//-----------------------------------
function [r, ires]=chemres(t, y, yd, a, b, c)
   r=[-a*y(1)+b*y(2)*y(3)-yd(1)
       a*y(1)-b*y(2)*y(3)-c*y(2)*y(2)-yd(2)
       y(1)+y(2)+y(3)-1];
   ires=0
endfunction
function pd=chemjac(x, y, yd, cj, a, b, c)
    pd=[-a-cj , b*y(3)             , b*y(2);
         a    ,-b*y(3)-2*c*y(2)-cj ,-b*y(2);
         1    , 1                  , 1       ]
endfunction
y3=dassl([y0,yd0],0,t,list(chemres,0.04,1d4,3d7),list(chemjac,0.04,1d4,3d7));
// Cコードを使用
//------------
// - Cコードを作成
rescode=['void chemres(double *t, double y[], double yd[], double r[], int *ires, double rpar[], int ipar[])'
         ' {'
         '   r[0] = -0.04*y[0]+1.0e4*y[1]*y[2]                -yd[0];'
         '   r[1] =  0.04*y[0]-1.0e4*y[1]*y[2]-3.0e7*y[1]*y[1]-yd[1];'
         '   r[2] =       y[0]+y[1]+y[2]-1;'
         '   *ires = 0;'
         ' }'];
jaccode=['void chemjac(double *t, double y[], double yd[], double pd[], double *cj, double rpar[], int ipar[])'
         ' {'
         '   /* first column*/'
         '   pd[0] = -0.04-*cj;'
         '   pd[1] =  0.04;'
         '   pd[2] =  1.0;'
         '    /* second column*/'
         '   pd[3] =  1.0e4*y[2];'
         '   pd[4] = -1.0e4*y[2]-2*3.0e7*y[1]-*cj;'
         '   pd[5] =  1.0;'
         '    /* third column*/'
         '   pd[6] =  1.0e4*y[1];'
         '   pd[7] = -1.0e4*y[1];'
         '   pd[8] =  1.0;'
         ' }'];
mputl([rescode;jaccode],fullfile(TMPDIR,'mycode.c')) //Cファイルを作成
// - コンパイル
ilib_for_link(['chemres','chemjac'],fullfile(TMPDIR,'mycode.c'),[],'c','',fullfile(TMPDIR,'loader.sce'));//compile
// - Scilabにリンク
exec(fullfile(TMPDIR,'loader.sce')) //incremental linking
// - dasslをコール
y4=dassl([y0,yd0],0,t,'chemres','chemjac');

参照

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