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

dverka() void dverka ( int  n, void(*)(int, double, double *, double *)  f, double *  t, double  y[], double  tout, double  tend, double  tol, int  mode, double  work[], int  lwork, int  iwork[], int  liwork, int *  info  ) 常微分方程式の初期値問題 (6(5)次 ルンゲ・クッタ・ヴァーナー法) 目的 本プログラムは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, y は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. また, t0 および y0 はそれぞれ t および y の既知の初期値である. 本プログラムは, 6(5)次 ルンゲ・クッタ・ヴァーナー法プログラム DVERK (文献 (1)) に Enrightら (文献 (2)) の補間式を使用した密出力機能を追加して書き直したものである. 引数 [in] n 微分方程式の数. (n >= 1) [in] f 微分方程式の関数値を求めるユーザー定義サブルーチンで, 次のように定義すること. void f(int n, double t, const double y[], double yp[]) { dy/dt を求め yp[] に入れる. } ただし, n は方程式の数, dy/dt (= f(t, y[])) は与えられた t および y[] における微分の計算値である. [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 <= tout < t でなければならない). 先に tend に達した場合にはその時点で積分を終了し t = tend, info = 0 として戻る. [in] tend 積分を終了する点を表す. tend に達すると積分を終了し, t = tend, info = 0 として戻る. 積分を行うのは前進方向 (tend > t) でも後退方向 (tend < t) でもよい. [in] tol 誤差許容値を表す. (tol > 0) 大域誤差が tol に比例するように局所誤差のノルムが制御される. ノルムは, iwork[11] パラメータで指定される誤差制御方式による重み付き最大ノルムが使用される. デフォルトの重みは 1/max(1, abs(y(i))) で, 絶対誤差制御および相対誤差制御を混合したものである. [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]: ステップ幅の初期値. (デフォルト値 = hmax*tol^(1/6)) work[2] 最大ステップ幅 (hmax) を決めるパラメータ. hmax = min(abs(work[2]), 2/abs(scale)) (ただし, scale = 0 のときは hmax = abs(work[2])). (デフォルト値: hmax = 2/abs(scale) (ただし, scale = 0 のときは hmax = 2)) work[11]: 最小ステップ幅 (hmin). (デフォルト値 = 10*max(dwarf, rreb*max((y の重み付きノルム)/tol, abs(t))). ただし, dwarf は非常に小さな正数(マシン精度) (= 1.0e-50), rreb は相対丸め誤差限界値である.) work[12]: 問題のスケール尺度 (scale). (デフォルト値 = 1)n 本パラメータは, hmax の決定および判定基準の変更のために使われる. この値が大きければより高精度の結果を求めようとする. work[13]: iwork[11] = 3 の場合の下限値. [in] lwork 配列 work[] のサイズ. (abs(lwork) >= 12*n + 40) lwork < 0 であればその絶対値を採用し, work[0]〜work[19] をすべて 0 に設定する. [in,out] iwork[] 配列 iwork[liwork] 整数作業領域. ただし, iwork[0]〜iwork[19] はプログラムのためのパラメータ領域として使用される. info = 0 の呼び出し時に 0 に設定されている入力パラメータはそのデフォルト値が使用される. [in] iwork[1]: 許される最大ステップ数. (デフォルト値 = 10000) iwork[11]: エラー制御方法. 誤差制御のための推定誤差ベクトルの重みを指定する. (デフォルト値 = 0) (値が下記以外の場合, デフォルト値が使用される) = 0: 重み = 1/max(1, abs(y(i))) (絶対および相対誤差制御を混合したもの) = 1: 重み = 1 (絶対誤差制御) = 2: 重み = 1/abs(y(i)) (相対誤差制御) = 3: 重み = 1/max(abs(work[13]), abs(y(i))) (abs(y(i)) が下限値 abs(work[13]) より小さくない限り相対誤差制御) 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) >= 40) 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: (エラー) ステップ幅の選択に失敗した (tol が小さすぎる可能性がある). = 13: (エラー) ステップ幅が小さくなりすぎた (tol が小さすぎる可能性がある). 文献 (1) netlib/ode (2) W. H. Enright, et al.: "Interpolants for Runge-Kutta Formulas", ACM Transactions on Mathematocal Software, Vol.12, No.3. (1986)