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

odex2a() void odex2a ( int  n, void(*)(int, double, double *, double *)  f2, double *  t, double  y[], double  yp[], double  tout, double  tend, double *  rtol, double *  atol, int  itol, int  mode, double  work[], int  lwork, int  iwork[], int  liwork, int *  info  ) 2階常微分方程式の初期値問題 (補外法) 目的 本プログラムは2階常微分方程式の初期値問題 d2y/dt2 = f(t, y), ただし t = t0 において y = y0, y' = y'0 の解を求める. ただし, y および y' は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. また, t0, y0 および y'0 はそれぞれ t, y および y' の既知の初期値である. 本プログラムは, 補外法 (シュテルマー公式) プログラム ODEX2 (文献 (1)) を書き直したものである. 引数 [in] n 微分方程式の数. (n >= 1) [in] f2 微分方程式の関数値を求めるユーザー定義サブルーチンで, 次のように定義すること. void f2(int n, double t, const double y[], double ypp[]) { d2y/dt2 を求め ypp[] に入れる. } ただし, n は方程式の数, d2y/dt2 (= f2(t, y[])) は与えられた t および y[] における2次微分の計算値である. [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,out] yp[] 配列 yp[lyp] (lyp >= n) [in] t の初期値 t0 における微分係数 y' の初期値 y'0. [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] rtol スカラー (itol = 0 の場合) または 配列 rtol[lrtol] (itol = 1 の場合) (lrtol >= 2*n) (rtol または rtol[i] >= 0) 求める解の精度を指定する相対誤差許容値を表す. 本パラメータは itol の指定によりスカラーまたは配列を選択できる. 許容値は atol と共に各ステップにおける局所誤差テストに用いられ, y[] および yp[] の各要素が次式を満たすように自動的に選ばれたステップ幅を用いて積分が行われる. (局所誤差) <= rtol*abs(y[i]) + atol および rtol*abs(yp[i]) + atol (itol = 0 の場合) (局所誤差) <= rtol[i]*abs(y[i]) + atol[i] および rtol[i + n]*abs(yp[i]) + atol[i + n] (itol = 1 の場合(i = 0 〜 n - 1)) rtol と atol (または rtol[i] と atol[i]) が同時に 0 であってはならない. [in] atol スカラー (itol = 0 の場合) または 配列 atol[latol] (itol = 1 の場合) (latol >= 2n) (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]: ステップ幅の初期値. (デフォルト値 = 0.0001) work[2]: ステップ幅の最大値 (デフォルト値 = abs(tend - t)) work[4] および work[5]: ステップ幅選択パラメータ. facmin/work[5] <= j番目のhnew/hold <= 1/facmin となるようにステップ幅が選ばれる. ただし, facmin = work[4]^(1/(2*j - 1)). (デフォルト値: work[4] = 0.02, work[5] = 4) work[6] および work[7]: 次数選択パラメータ. w[k-2] <= w[k-1]*work[6] であれば次数を減らし, w[k-1] <= w[k-2]*work[7] であれば次数を増やす. (デフォルト値: work[6] = 0.8, work[7] = 0.9) work[8] および work[9]: ステップ幅制御の安全係数. hnew = h*work[9]*(work[8]*tol/err)^(1/(j-1)) とする. (デフォルト値: work[8] = 0.65, work[9] = 0.94) work[10]: ステップ幅を減らすときの係数. (デフォルト値 = 0.5) [out] work[1]: 最後に使用されたステップ幅. [in] lwork 配列 work[]のサイズ. (abs(lwork) >= n*(7*km + 14) + km*(2*n*km + 5) + 40 (km = iwork[2])) lwork < 0 であればその絶対値を採用し, work[0]〜work[19] をすべて 0 に設定する. [in,out] iwork[] 配列 iwork[liwork] 作業領域. ただし, iwork[0]〜iwork[19] はプログラムのためのパラメータ領域として使用される. info = 0 の呼び出し時に 0 に設定されている入力パラメータはそのデフォルト値が使用される. [in] iwork[1]: 許される最大ステップ数 (デフォルト値 = 10000) (iwork[1] < 0 の場合, デフォルト値を使用する) iwork[2]: 補外表の最大列数 km (km >= 3). (デフォルト値 = 9) (km < 3 の場合, デフォルト値を使用する.) iwork[3]: ステップ数列の選択. (デフォルト値 = 1 (mode != 3 の場合), 4 (mode = 3 の場合)) = 1: 2, 4, 6, 8, 10, 12, 14, 16, ... = 2: 2, 4, 8, 12, 16, 20, 24, 28, ... = 3: 2, 4, 6, 8, 12, 16, 24, 32, ... = 4: 2, 6, 10, 14, 18, 22, 26, 30, ... mode = 3 の場合, iwork[3] >= 4 でなければならない. iwork[3] の値が範囲外の場合, デフォルト値が使用される. iwork[6]: 補外式の次数を定めるパラメータ (1 <= iwork[6] <= 8). (デフォルト値 = 6) iwork[6] の値が範囲外の場合, デフォルト値が使用される) iwork[8]: 密出力式の誤差推定. mode = 3 のときに有効. (デフォルト値 = 0) = 0: 誤差推定を行う. = 1: 誤差推定を行わない. 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) >= km + 60 (km = iwork[2])) 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: (エラー) 最大ステップ数を超えた. 文献 (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))