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

dopn86() void dopn86 ( 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  ) 常微分方程式の初期値問題 (8(6)次ルンゲ・クッタ・ニュストレム法) (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' の既知の初期値である. 本プログラムは 8(6)次ルンゲ・クッタ・ニュストレム法プログラムである. 係数 (RKN8(6)9FM) は文献(1)を, 密出力のための補間式は文献(2)を参照せよ. 引数 [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' の値(数値解). 注 - 本プログラムでは mode = 3 で y' の値(補間値 yp[])を返さない. y だけでなく y' の値も必要な場合には mode = 2 を使用すること. [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 >= n) (rtol または rtol[i] >= 0) 求める解の精度を指定する相対誤差許容値を表す. 本パラメータは itol の指定によりスカラーまたは配列を選択できる. 許容値は atol と共に各ステップにおける局所誤差テストに用いられ, y[] および yp[] の各要素が次式を満たすように自動的に選ばれたステップ幅を用いて積分が行われる. itol = 0 の場合 (i = 0 〜 n - 1): (局所誤差) <= max(rtol*abs(y[i]) + atol, rtol*abs(yp[i]) + atol) (iwork[11] = 0 の場合) (局所誤差) <= rtol*abs(y[i]) + atol (iwork[11] = 1 の場合) itol = 1 の場合 (i = 0 〜 n - 1): (局所誤差) <= max(rtol[i]*abs(y[i]) + atol[i], rtol[i]*abs(yp[i]) + atol[i]) (iwork[11] = 0 の場合) (局所誤差) <= rtol[i]*abs(y[i]) + atol[i] (iwork[11] = 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 と同じステップを踏む. 注 - 本プログラムでは mode = 3 で y' の値(補間値 yp[])を返さない. y だけでなく y' の値も必要な場合には mode = 2 を使用すること. [in,out] work[] 配列 work[lwork] 作業領域. ただし, work[0]〜work[19] はプログラムのためのパラメータ領域として使用される. info = 0 の呼び出し時に 0 に設定されている入力パラメータはそのデフォルト値が使用される. [in] work[0]: ステップ幅の初期値 (デフォルト値 = プログラムにより自動推定される) work[2]: ステップ幅の最大値 (デフォルト値 = abs(tout - 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) [out] work[0]: 最後に使用されたステップ幅. [in] lwork 配列 work[]のサイズ. (lwork >= 15*n + 40) lwork < 0 であれば abs(lwork) を使用し, 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: y および y' の両方の推定誤差を使用してステップの採用を判定する. = 1: y の推定誤差のみを使用してステップの採用を判定する. 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[]のサイズ (liwork >= 30) liwork < 0 であれば abs(liwork) を使用し, 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: (エラー) ステップ幅が小さくなりすぎたため計算が継続できない. 参考文献 (1) J. R. Dormand, M. E. A. El-Mikkawy, P. J. Prince, "High Order Embedded Runge-Kutta-Nystrom Formulae", IMA J. of Numerical Analysis, 7, 423-430 (1987). (2) Ch. Tsitouras, G. Papageorgiou "Interpolating Runge-Kutta-Nystrom Methods of High Order", Intern. J. Computer Math., Vol. 41, 209-217 (1993).