deabm deabm_r derkfa derkfa_r dop853a dop853a_r dopn1210 dopn1210_r dopn43 dopn43_r dopn64 dopn64_r dopn86 dopn86_r dopri5a dopri5a_r dverka dverka_r odex2a odex2a_r odexa odexa_r retarda retarda_r ylaga

retarda_r() void retarda_r ( int  n, double *  t, double  y[], double  tout, double  tend, double *  rtol, double *  atol, int  itol, int  mode, double  work[], int  lwork, int  iwork[], int  liwork, int *  info, double *  tt, double  yy[], double  yyp[], int *  irev  ) 遅延微分方程式の数値解 (5(4)次 ドルマン・プリンス法) (リバースコミュニケーション版) 目的 本プログラムは1階の遅延微分方程式の初期値問題 dy/dt = f(t, y(t), y(t - τ), y(t - τ2), ...), ただし t = t0 において y = y0 の解を求める. ただし, y は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. また, t0 および y0 はそれぞれ t および y の既知の初期値である. f は過去の値 y(t - τ), y(t - τ2), ... にも依存する. 本プログラムは, 遅延微分方程式用の 5(4)次 ドルマン・プリンス法プログラム RETARD (文献 (1)) を書き直したものである. 本プログラムは retarda のリバースコミュニケーション版である. 引数 [in] n 微分方程式の数. (n >= 1) [in,out] t 独立変数 t を表す. 本プログラムは t の初期値から tend までの積分を行う. mode の設定により, 途中結果を返すために tout またはステップごとに戻ることができる. [in] t の初期値 t0. [out] 積分終了時には t = tend, 途中結果を返すために戻ったときには mode の設定により t = tout または t = 直前のステップの終点の値 を返す. [in,out] y[] 配列 y[ly] (ly >= n) 従属変数 y を表す. [in] t の初期値 t0 における y の初期値 y0. [out] t における y の値 (数値解). [in] tout mode = 2 または 3 の場合に, 途中結果の確認／出力を行う t を表す. mode = 0 または 1 では tout は参照されない. mode = 2 または 3 では t = tout における y を求め, それぞれ info = 2 または 3 として戻る. 続いて新たな tout における解を求めるために積分を継続したい場合には, 新たな tout に変更して(info を含む)他の変数を変更せずに再度呼び出しを行うことができる. t < tout <= tend でなければならない. 先に tend に達した場合にはその時点で積分を終了し t = tend, info = 0 として戻る. [in] tend 積分を終了する点を表す. tend に達すると積分を終了し, t = tend, info = 0 として戻る. 積分を行うのは前進方向 (tend > t) でなければならない. [in] rtol スカラー (itol = 0 の場合) または 配列 rtol[lrtol] (itol = 1 の場合) (lrtol >= n) (rtol または rtol[i] >= 0) 求める解の精度を指定する相対誤差許容値を表す. 本パラメータは itol の指定によりスカラーまたは配列を選択できる. 許容値は atol と共に各ステップにおける局所誤差テストに用いられ, y[] の各要素が次式を満たすように自動的に選ばれたステップ幅を用いて積分が行われる. (局所誤差) <= rtol*abs(y[i]) + atol (itol = 0 の場合) (局所誤差) <= rtol[i]*abs(y[i]) + atol[i] (itol = 1 の場合(i = 0 〜 n - 1)) rtol と atol (または rtol[i] と atol[i]) が同時に 0 であってはならない. [in] atol スカラー (itol = 0 の場合) または 配列 atol[latol] (itol = 1 の場合) (latol >= n) (atol または atol[i] >= 0) 求める解の精度を指定する絶対誤差許容値を表す. 本パラメータは itol の指定によりスカラーまたは配列を選択できる. 許容値は rtol と共に各ステップにおける局所誤差テストに用いらる (上記 rtol 参照). [in] itol rtol および atol がスカラーか配列かを指定する. = 0: rtol および atol はスカラー (または要素数 1 の配列). = 1: rtol および atol は配列. [in] mode 動作モード. 本プログラムは t0 から tend までの積分を行うが, 中間結果の確認／出力の方法により4つの動作モードが提供される. どのモードでも t = tend に達したときには積分を終了し, info = 0 を返す. = 0: tend まで戻らない. tout は無視される. = 1: 成功したステップごとに戻る (info = 1 を返す). tout は無視される. t にはそのステップの終了点を返す. 最終ステップでは t = tend となるようにステップ幅が調整される. = 2: 積分途中で tout における中間結果を返すために戻る (t = tout, info = 2 を返す). 続けて次の tout における解を求めるために積分を継続するためは, tout を新たな値に変えて再度呼び出しを行うことができる. tout 直前のステップでは t = tout となるようにステップ幅が調節される. = 3: 積分途中で tout における中間結果を返すために戻る (t = tout, info = 3 を返す). 続けて次の tout における解を求めるために積分を継続するためは, tout を新たな値に変えて再度呼び出しを行うことができる. mode = 2 と異なり tout における y の値は補間により求められ, tout 直前のステップでのステップ幅の調節は行わず mode = 1 と同じステップを踏む. [in,out] work[] 配列 work[lwork] 作業領域. ただし, work[0]〜work[19] はプログラムのためのパラメータ領域として使用される. info = 0 の呼び出し時に 0 に設定されている入力パラメータはそのデフォルト値が使用される. [in] work[0]: ステップ幅の初期値. (デフォルト値 = プログラムが自動推定する) work[2]: ステップ幅の最大値 (デフォルト値 = abs(tend - t)) work[4] および work[5]: ステップ幅選択パラメータ. work[4] <= hnew/hold <= work[5] となるようにステップ幅が選ばれる. (デフォルト値: work[4] = 0.2, work[5] = 10) work[8]: ステップ幅選択時の安全係数. 0.0001 < work[8] < 1 でなければならない. この範囲外であればデフォルト値が使用される. (デフォルト値 = 0.9) work[11]: ステップ幅制御の安定化のためのパラメータ. 0 < work[11] <= 0.2 でなければならない. この範囲外であればデフォルト値が使用される. (デフォルト値 = 0.04 (work[11] = 0 または work[11] > 0.2 の場合), 0 (work[11] < 0 の場合)) [out] work[1]: 最後に使用されたステップ幅 [in] lwork 配列 work[]のサイズ. (abs(lwork) >= 8*n + ngrid + mxst*(5*n + 1) + 41 (ngrid = iwork[5], mxst = iwork[6])) lwork < 0 であればその絶対値を採用し, work[0]〜work[19] をすべて 0 に設定する. [in,out] iwork[] 配列 iwork[liwork] 整数作業領域. ただし, iwork[0]〜iwork[19] はプログラムのためのパラメータ領域として使用される. info = 0 の呼び出し時に 0 に設定されている入力パラメータはそのデフォルト値が使用される. [in] iwork[1]: 許される最大ステップ数. (デフォルト値 = 100000) iwork[3]: iwork[3]回のステップごとにスティフ性のテストを行う. iwork[3]が負の値であればテストは行わない. (デフォルト値 = 1000) iwork[5]: グリッド点の数 ngrid. (デフォルト値 = 0) 積分区間内に解の導関数の不連続点(グリッド点)がある場合, それを指定して計算精度および効率を改善することができる. グリッド点を, work[41], ..., work[40 + ngrid] に設定しておかなければならない(積分区間の始点と終点は除く). iwork[6]: 過去何ステップ分の情報を保持しておくかを指定する. ylaga はこれを使用して過去の値 y(t - τ) を求める. (デフォルト値 = 1000) iwork[12]: カウンタ iwork[13]〜iwork[17] のリセット方法を指定する. (デフォルト値 = 0) = 0: info = 0 のときにリセットする. = 1: info = 0 であってもリセットしない. [out] 統計情報(カウンタ). iwork[13]: 関数評価回数. iwork[15]: 全ステップ数. iwork[16]: 誤差評価により採用されたステップ数. iwork[17]: 誤差評価により不採用だったステップ数. [in] liwork 配列 iwork[]のサイズ. (abs(liwork) >= 50) liwork < 0 であればその絶対値を採用し, iwork[0]〜iwork[19] をすべて 0 に設定する. [in,out] info [in] 制御コード. = 0: 最初の呼び出し時(新たに問題を開始する場合)には info = 0 と設定する. 全ての変数の初期化を行ってから計算を開始する. = 1, 2, 3: info = 1, 2 または 3 で戻った場合, 計算を継続するために tout だけを変更し, info の値をそのままにして再呼び出しすることができる. [out] リターンコード. = 0: 正常終了. tend までの積分が完了した. < 0: (-info)番目の入力パラメータの誤り. = 1: mode = 1 の中間結果出力のために戻った. 再呼び出しすることより次のステップまで進むことができる. = 2: mode = 2 の中間結果出力のために t = tout において戻った. tout を再設定して再呼び出しすることができる. = 3: mode = 3 の中間結果出力のために t = tout において戻った. tout を再設定して再呼び出しすることができる. = 11: (エラー) 最大ステップ数を超えた. = 12: (エラー) ステップ幅が小さくなりすぎたため計算が継続できない. = 13: (エラー) スティフな問題の可能性がある. [out] tt irev = 1〜9 の場合: yyp[]を参照のこと. [out] yy[] 配列 yy[lyy] (lyy >= n) irev = 1〜9 の場合: yyp[]を参照のこと. [in] yyp[] 配列 yyp[lyyp] (lyyp >= n) irev = 1〜9 の場合: t = tt および y = yy[] における微分値 dy/dt = f(t, y(t), y(t - τ), y(t - τ2), ...) を計算して, 再呼び出し時に yyp[] に設定しておくこと. y(t - τ), y(t - τ2), ... は ylaga を使用して求めることができる. ylaga の引数 cont[] および icont[] には work[] および iwork[] を渡すこと. [in,out] irev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行い irev を変更せずに再び本プログラムを呼び出すこと. = 0: 処理終了. 正常終了かどうかは info をチェックすること. = 1〜9: t = tt および y = yy[] における微分の計算値を yyp[] に設定する. yyp[] 以外の変数を変更してはならない. 文献 (1) 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))