_deabm _deabm_r _derkf _derkf_int _derkf_r _dop853 _dop853_r _dopri5 _dopri5_r _doprin _doprin_r _dverk _dverk_int _dverk_r _odex _odex2 _odex2_r _odex_r _retard _retard_r

_retard_r() void _retard_r ( int  n, double *  t, double  y[], double  tout, double *  rtol, double *  atol, int  itol, int  iout, double  work[], int  lwork, int  iwork[], int  liwork, double  rcont[], int  lrcont, int  icont[], int  licont, int *  info, double *  tt, double  yy[], double  yyp[], int *  irtrn, int *  irev  ) 遅延微分方程式の数値解 (5(4)次 ドルマン・プリンス法) (リバースコミュニケーション版) 目的 本ルーチンは1階の遅延微分方程式の初期値問題 dy/dt = f(t, y(t), y(t-τ), ...), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. retardは5(4)次ドルマン・プリンス法に基づく陽的ルンゲ・クッタ法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている. 詳細については下記参考文献を参照せよ. retard_rはretardのリバースコミュニケーション版である. 引数 [in] n 微分方程式の数. (n >= 1) [in,out] t 本ルーチンはtからtoutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. [in] 独立変数 t の初期値. [out] ステップ終了時の最終値(通常toutに等しい). 解がこの点まで正常に求められたことを示す. 新たなtoutで再呼び出しすることにより, 新たな点まで積分を継続することができる. [in,out] y[] 配列 y[ly] (ly >= n) [in] tの初期値における従属変数y[]の初期値. [out] 最終のt(通常toutに等しい)において求められた解(の近似値). [in] tout 解を求めたい点を設定する. 積分を行うのはtについて前進方向(tout > t)でなければならない. 要求精度を満たすように自動的に選ばれたステップ幅を用いてtからtoutに向かって解が求められる. [in] rtol スカラー(itol = 0 の場合) または 配列 rtol[lrtol] (itol = 1 の場合) (lrtol >= n) (rtolまたはrtol[i] >= 0) 相対誤差許容値. y[i]の局所誤差がおおむね次式以下になるようにする. rtol*abs(y[i]) + atol (i = 0〜n-1) (itol = 0 の場合) または rtol[i]*abs(y[i]) + atol[i] (i = 0〜n-1) (itol = 1 の場合). rtolとatol, または, rtol[i]とatol[i] (i = 0〜n-1) が同時に0になってはいけない. [in] atol スカラー(itol = 0 の場合) または 配列 atol[latol] (itol = 1 の場合) (latol >= n) (atolまたはatol[i] >= 0) 絶対誤差許容値. y[i]の局所誤差がおおむね次式以下になるようにする. rtol*abs(y[i]) + atol (i = 0〜n-1) (itol = 0 の場合) または rtol[i]*abs(y[i]) + atol[i] (i = 0〜n-1) (itol = 1 の場合). rtolとatol, または, rtol[i]とatol[i] (i = 0〜n-1) が同時に0になってはいけない. [in] itol rtolおよびatolがスカラーか配列かを指定する. = 0: rtol および atolはスカラー. = 1: rtol および atolは配列. (他の値であれば itol = 0 とみなす.) [in] iout 中間結果を出力するために戻るかどうか指定 (irev = 50, 51を参照). = 0: 戻らない = 1: 密出力付きの中間出力のために戻る. この場合, iwork[4]に密出力を行う要素数 (nrdens)を設定しなければならない. (他の値であれば iout = 0 とみなす.) [in,out] work[] 配列 work[lwork] 作業領域. work[0]〜work[19]はプログラムのためのパラメータ領域として使用される. 入力パラメータを0に設定すると, パラメータ値としてはそれぞれのデフォルト値を使用する. [in] work[0]: ステップ幅の初期値 (デフォルト値 = プログラムにより自動推定される) work[1]: ステップ幅の最大値 (デフォルト値 = tout - t) work[3] および work[4]: ステップ幅選択パラメータ. work[3] <= hnew/hold <= work[4] となるようにステップ幅が選ばれる. (デフォルト値: work[3] = 0.2, work[4] = 10) work[7]: ステップ幅推定時の安全係数 (0.0001 < work[7] < 1). work[7] <= 0.0001 または work[7] >= 1 であればデフォルト値を使用する. (デフォルト値 = 0.9) work[8]: ステップ幅制御の安定化のためのパラメータβ (work[8] <= 0.2). work[8] < 0であればβは0に設定される. (デフォルト値 = 0.04) work[20], ..., work[20 + ngrid - 1]: 積分の際にグリッド点とみなすべき点を指定する (t < work[20] < ... < work[20 + ngrid - 1] = tout). [out] work[0]: 最後に使用されたステップ幅 [in] lwork 配列 work[]のサイズ. (lwork >= 8*n + ngrid + 20 (ngrid: iwork[5]を参照)) lwork < 0 であれば abs(lwork) を使用し, work[0]〜work[19]を 0 に初期化する. [in,out] iwork[] 配列 iwork[liwork] 作業領域. iwork[0]〜iwork[19]はプログラムのためのパラメータ領域として使用される. 入力パラメータを0に設定すると, パラメータ値としてはそれぞれのデフォルト値を使用する. [in] iwork[1]: 許される最大ステップ数. (デフォルト値 = 100000) (iwork[1] < 0 の場合, デフォルト値を使用する) iwork[3]: iwork[3]回のステップごとにスティフ性のテストを行う. iwork[3]が負の値であればテストは行わない. (デフォルト値 = 1000) iwork[5]: ngrid = グリッド点の数. (デフォルト値 = 1) 積分区間中に解の導関数の不連続点(グリッド点)がある場合, それを指定しておくと計算精度および効率を改善することができる. iwork[5] <= 0 であれば, ngrid = 1 とみなす. iwork[5] > 0 であれば, work[20], ..., work[20 + ngrid - 1] を設定しておかなければならない. 最後の点 work[20 + ngrid - 1] = tout でなければならない(iwork[5] = 0 のときには自動的に設定される). iwork[7]: nrdens = 密出力を行う要素数(fおよびsoloutで使用). iwork[7] < 0 であれば, nrdens = abs(iwork[7]) とする. (デフォルト値 = n) 0 < nrdens < n の場合, icont[]を設定しなければならない(下記参照). iwork[12]: iwork[13]〜iwork[17]をリセットするタイミングを指定する. = 0: 本ルーチンが呼び出されるたびにリセットする. != 0: info = 0 のときにリセットする. [out] iwork[13]: nfcn = 関数評価回数 iwork[15]: nstep = ステップの計算回数 iwork[16]: naccept = 採用されたステップ数 iwork[17]: nreject = 不採用だったステップ数 (誤差評価の結果による) (最初のステップの不採用はカウントされない) [in] liwork 配列 iwork[]のサイズ. (liwork >= 20) liwork < 0 であれば abs(liwork) を使用し, iwork[0]〜iwork[19]を 0 に初期化する. [out] rcont[] 配列 rcont[lrcont] y(t-τ)の計算および密出力の制御情報のための領域. [in] lrcont 配列 rcont[]のサイズ. (lrcont >= mxst*(5*nrdens + 2) (nrdens: iwork[7]を参照, mxst: 記録する(ylagによりさかのぼることができる)ステップ数)) [in,out] icont[] 配列 icont[licont] y(t-τ)の計算または密出力を行う要素番号のリスト. nrdens = n または iwork[7] < 0 の場合, 自動的に 1, ..., nrdens と設定される. 0 < nrdens < n の場合, 密出力を行う要素番号を設定しておかなければならない. nrdensについてはiwork[7]を参照せよ. [in] licont 配列 icont[]のサイズ. (licont >= nrdens (nrdens: iwork[7]を参照)) [in,out] info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). ngrid = 1 のときにのみ使用可能. [out] = -1: 入力パラメータ n の誤り (n < 1) = -5: 入力パラメータ rtol の誤り (rtol < 0 または rtol[i] < 0) = -5: 入力パラメータ rtol または atol の誤り (rtol = 0 かつ atol = 0, または rtol[i] = 0 かつ atol[i] = 0) = -6: 入力パラメータ atol の誤り (atol < 0 または atol[i] < 0) = -9: 入力パラメータ work[20 + i - 1] の誤り (昇順でない) = -10: 入力パラメータ lwork の誤り (lworkが小さすぎる) = -12: 入力パラメータ liwork の誤り (liworkが小さすぎる) = -14: 入力パラメータ lrcont の誤り (lrcontが小さすぎる) = -16: 入力パラメータ licont の誤り (licontが小さすぎる) = -17: 入力パラメータ info の誤り (info != 0 かつ info != 1) = -22: 入力パラメータ irev の誤り (irev != 1 〜 9, 50, 51) = 1: 正常終了 = 2: irtrn による中断 (正常終了) = 11: 最大ステップ数を超えた = 12: ステップサイズが小さくなり過ぎた = 13: スティフな問題の可能性がある (中断された) = 14: ylagにより中断された (mxstが小さすぎるなどの原因が考えられる) [out] tt irev = 1〜9: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のtの値を返す. [out] yy[] 配列 yy[lyy] (lyy >= n) irev = 1〜9: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のyの値を返す. [in] yyp[] 配列 yyp[lyyp] (lyyp >= n) irev = 1〜9: t(= tt)およびy(= yy[])における微分係数, すなわち, yyp[i] = dyi/dt = fi(t, y(t), y(t-τ), ...) (i = 0 〜 n-1) を計算して, 再呼び出し時に設定すること. y(t - τ) は ylag_r 使って求めることができる. ylag_r(i, x, rcont, icont, &yy, &irev); ただし, i は要素番号(0 〜 n - 1), x は t - τ である. rcont および icont は変更せずに渡さなければならない. irev = 0 のときに yy に結果が入る. t0 - τ <= x <= t0 (ただし, t0 は t の初期値) であれば irev = 1 または 2 で戻るので, x における y[i] の初期値を yy に与えること. [in,out] irtrn [in] irev = 50, 51: 積分を中断したい場合以外変更してはならない. irtrnを負の値に設定すると, 積分を中断して info = 2 で終了する. [out] irev = 50, 51: irev = 50 または 51で最初に戻ったときに0, 最後に戻ったときに2, それ以外のときに1を返す. [in,out] irev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはinfoをチェックすること. = 1〜9: ttおよびyy[]における微分係数の計算値をyyp[]に設定すること. yyp[]以外の変数を変更してはならない. = 50, 51: iout = 1の場合, 中間結果を出力するためにステップが成功するごとにこのコードで戻る. tにおけるy[]の値が渡される. ttは前回のtの値である. iout = 1 の場合には, 制御情報rcont[]およびicont[]により密出力を行うことができる. t以前の任意の点での解y[i] (0 <= i <= n-1)をylag_r()を呼び出すことにより求めることができる. yyp[]の説明を参照せよ. 出典 E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations. Nonstiff Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1993) 邦訳: 「常微分方程式の数値解法Ⅰ 基礎編」スプリンガージャパン (2007)