daskr

引数

x0

は y0 (ydot0 は ゼロを初期推定値として dassl により推定) または行列 [y0 ydot0]のどちらかです. g(t,y0,ydot0) はゼロに等しい必要があります. ydot0 の推定値のみ既知の場合, info(7)=1 を指定してください.

y0

初期条件の実数列ベクトル.

ydot0

y at t0 (推定値でも可)におけるyの時間微分の実数列ベクトル.

t0

実数, 初期時間.

t

実数スカラーまたはベクトル. 解を得る時刻を指定します. info(2)=1と指定することにより, dasslの各ステップ点で解を得ることができます.

nn

2個のエントリ [times num]を有するベクトル. times は面が交差する点における時間の値, num は交差する面の数です.

atol, rtol

yと同じ大きさの実数スカラーまたは列ベクトル. atol, rtol はそれぞれ解の絶対及び相対許容誤差を 表します. この許容誤差は yの各要素について指定されます.

res

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

• Scilab関数.

  呼び出し手順は以下とします [r,ires]=res(t,y,ydot) かつ res は残差 r=g(t,y,ydot) かつ エラーフラグ iresを返すこと. res は, rの計算に成功した場合はires = 0, 残差が(t,y,ydot) において局所的に定義されないには =-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

external (関数またはリストまたは文字列). 指定した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にリンクされたCまたは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

external (関数またはリストまたは文字列). k 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

14 個の要素を有するリスト. デフォルト値は list([],0,[],[],[],0,[],0,[],0,[],[],[],1)です.

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]を指定します. info(3)=[]の場合, jac は通常の行列 を返します. info(8)=1の場合, ダミーとして扱われます.

info(4)

最大ステップサイズを指定する実数スカラー値. 制限しない場合は info(4)=[]を指定します.

info(5)

初期ステップサイズを指定する実数スカラー値. 指定しない場合は info(5)=[]を指定します.

info(6)

解が非負であることが既知な場合は info(6)=1, それ以外は info(6)=0 を指定します.

info(7)

ydot0 がg(t0,y0,ydot0)=0のように指定された場合, info(7)=[] と指定します. それ以外の場合は info(7)=[+-1, ..., +-1]とします. ここで,y(i)が微分変数の場合はinfo(7)(i)=1, y(i)が代数変数(関数g(t,y,ydot)に微分値が陽に現れない)の場合は info(7)(i)=-1とします.

info(8)

直接法 / Krylov. ソルバーでKrylov反復法を使用したい場合,info(8)=1と指定し, psolにルーチンを指定します. それ以外の場合(dasslの直接法)はinfo(8)=0を指定します.

info(9)

Krylovパラメータ. info(8)=0とした場合はダミー引数として処理されます. それ以外の場合は info(9)=[maxl kmp nrmax epli]を指定します, ただし:

- maxl = mSPIGMRアルゴリズムの最大反復回数 (デフォルト min(5,neq)),

- kmp = SPIGMRアルゴリズムで直交化を行うベクトルの数 (デフォルト maxl),

- nrmax = 各非線形反復におけるSPIGMRアルゴリズムのリスタート最大数 (デフォルト 5),

- epli = SPIGMRアルゴリズムの収束安定定数 (デフォルト 0.05).

info(10)

初期条件. info(7)=[]の場合はダミーとして処理されます. 初期条件の計算直後にソルバーを停止させる場合は info(10)=1と指定し,それ以外は info(10)=0と指定します.

info(11)

psol用のプリコンディショナ計算およびLU分解ルーチン. info(8)=0の場合はダミーとして扱われます. external psolが特定の ルーチンを使用する必要がある場合は,info(11)=1を指定し, pjacルーチンを指定します. それ以外はinfo(11)=0とします.

info(12)

全変数についてエラー制御をローカルに行いたい場合, info(12)=[]を指定します. それ以外は, info(12)=[+-1, ..., +-1]を指定します. ただし, y(i)が微分変数の場合は info(12)(i)=1, y(i)が代数変数(微係数が関数g(t,y,ydot)に陽に現れない) の場合は info(12)(i)=-1とします.

info(13)

経験的なパラメータ. info(7)=[]の場合はダミーとして扱われます. それ以外の場合, info(13)=[mxnit mxnj mxnh lsoff stptol epinit]を指定します, ただし:

- mxnit = ヤコビアンまたはプリコンディショナ評価毎のニュートン反復の最大回数 (デフォルト:info(8)=0の場合は 5, それ以外は 15),

- mxnj = ヤコビアンまたはプリコンディショナ評価の最大回数 (デフォルト:info(8)=0の場合は 6, それ以外は 2),

- mxnh = info(7) ≠ []　の場合に試行する人工的なステップサイズ パラメータHの最大値 (デフォルト 5),

- lsoff = ライン探索アルゴリズムをオフにするフラグ (lsoff = 0 はライン探索有効を意味します, lsoff = 1 はオフを意味します) (デフォルト: 0),

- stptol = ライン探索アルゴリズムの最小スケール刻み (デフォルト: (丸め誤差単位)^(2/3)),

- epinit = ニュートン反復収束判定の振幅係数 (デフォルト 0.01).

info(14)

冗長性. 出力を最小化する場合は info(14)=1, 最大化する場合は info(14)=2, それ以外は info(14)=0を指定します.

psol

external (関数またはリストまたは文字列). 線形システムP*x = bを解きます. ただし, P はプリコンディショナで, ルーチン pjac が事前に計算して, wp 及び iwpに保存します.

• Scilab関数.

  呼び出し手順は [r, ier] = psol(wp, iwp, b) とします. psol 関数はシステムの解rと エラーフラグierを返す必要があります.

• リスト.

  以下のような形式とします

  list(psol,x1,x2,...)
  

  ただし, 関数 psolの呼び出し手順は以下とします

  psol(wp,iwp,b,x1,x2,...)
  

  psol は rの解を返します.

• 文字列.

  ScilabにリンクされたCまたはFortranサブルーチンの名前を指定します

  Cの場合,呼び出し手順は以下とします:

  void psol (int*neq, double*t, double*y, double*ydot, double*savr,
    double*wk, double*cj, double*wght, double*wp, int*iwp, double*b, double*eplin, int*ier, double*rpar, int*ipar)
  

  ただし,配列 wp および iwp はプリコンディショナ Pの行列要素を行圧縮形式で有します.

  Fortranの場合は以下とします:

  subroutine psol (neq, t, y, ydot, savr, wk, cj, wght,
    wp, iwp, b, eplin, ier, rpar, ipar)
  double precision t,y(*),ydot(*),savr(*),wk(*),cj,wght(*),wp(*),
    b(*),eplin,rpar(*)
  integer neq,iwp(*),ier,ipar(*)
  

pjac

external (関数またはリストまたは文字列). 指定したパラメータcjについて dg/dy + cj*dg/dydotを計算し, doubleおよびint配列にLU分解を行います.

• Scilab関数.

  呼び出し手順は [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr) となり, 結果として, 配列 wp および iwp には全てのプリコンディショナ情報が 圧縮スパース行形式で返されます.

• リスト.

  以下の形式とします

  list(pjac,x1,x2,...)
  

  ただし,関数 pjacの呼び出し手順は以下とします

  pjac(neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)
  

  pjac は, (neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)の関数として 分解した dg/dy + cj*dg/dydot を返します.

• 文字列.

  ScilabにリンクされたCまたはFortran形式のサブルーチン名を指定します

  Cの場合,呼び出し手順は以下とします:

  void pjac (double*res, int*ires, int*neq, double*t, double*y, double*ydot, double*rewt, double*savr,
  double*wk, double*h, double*cj, double*wp, int*iwp, int*ier, double*rpar, int*ipar)
  

  Fortranの場合, 以下のとします:

  subroutine pjac (res, ires, neq, t, y, ydot, rewt, savr,
    wk, h, cj, wp, iwp, ier, rpar, ipar)
  double precision res(*), t, y(*), ydot(*), rewt(*), savr(*),
    wk(*), h, cj, wp(*), rpar(*)
  integer ires, neq, iwp(*), ier, ipar(*)
  

hd

dasslコンテキストを保存し,積分を再開できるようにする実数ベクトル.

r

実数行列 . 各列はベクトル[t;x(t);xdot(t)]で, t は解が計算される添字の時刻です.

説明

陰的微分方程式の解.

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

局面交差した時刻と面交差した数をnnに返します.

詳細な例は SCI/modules/differential_equations/tests/unit_tests/daskr.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] = daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1);
nn
// nn = [2.4698972 2] を返します

参照

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

履歴