|
|
◆ _doprin_r()
| void _doprin_r |
( |
int |
n, |
|
|
double * |
t, |
|
|
double |
y[], |
|
|
double |
yp[], |
|
|
double |
tout, |
|
|
double * |
rtol, |
|
|
double * |
atol, |
|
|
int |
itol, |
|
|
int |
iout, |
|
|
double |
work[], |
|
|
int |
lwork, |
|
|
int |
iwork[], |
|
|
int |
liwork, |
|
|
int * |
info, |
|
|
double * |
tt, |
|
|
double |
yy[], |
|
|
double |
yypp[], |
|
|
int * |
irtrn, |
|
|
int * |
irev |
|
) |
| |
常微分方程式の初期値問題 (7(6)次ルンゲ・クッタ・ニュストレム法) (2階微分方程式用) (リバースコミュニケーション版)
- 目的
- 本ルーチンは2階の常微分方程式の初期値問題
d2y/dt2 = f(t, y), ただし t = t0 において y = y0, y' = y'0
本ルーチンはドルマンとプリンスによる7(6)次ルンゲ・クッタ・ニュストレム法のプログラムである. ステップ制御アルゴリズムを持っている.
詳細については下記参考文献を参照せよ.
doprin_rはdoprinのリバースコミュニケーション版である.
- 引数
-
| [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,out] | yp[] | 配列 yp[lyp] (lyp >= n)
[in] tの初期値における微分値y'1〜y'nの初期値.
[out] 最終のt(通常toutに等しい)において求められた微分値(の近似値). |
| [in] | tout | 解を求めたい点を設定する. 積分を行うのはtについて前進方向(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 | soloutの呼び出し指定.
= 0: soloutの呼び出しを行わない
= 1: soloutの呼び出しを行う
(他の値であれば 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)
[out]
work[0]: 最後に使用されたステップ幅 |
| [in] | lwork | 配列 work[]のサイズ. (lwork >= 6*n + 20)
lwork < 0 であれば abs(lwork) を使用し, work[0]〜work[19]を 0 に初期化する. |
| [in,out] | iwork[] | 配列 iwork[liwork]
作業領域.
iwork[0]〜iwork[19]はプログラムのためのパラメータ領域として使用される. 入力パラメータを0に設定すると, パラメータ値としてはそれぞれのデフォルト値を使用する.
[in]
iwork[1]: 許される最大ステップ数. (デフォルト値 = 10000) (iwork[1] < 0 の場合, デフォルト値を使用する)
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 に初期化する. |
| [in,out] | info | [in]
= 0: 初期化して計算を開始する(新たな問題を解く).
= 1: toutだけを変えて計算を続ける(前回呼び出しの計算を再開する).
[out]
= -1: 入力パラメータ n の誤り. (n < 1)
= -6: 入力パラメータ rtol の誤り (rtol < 0 または rtol[i] < 0)
= -6: 入力パラメータ rtol または atol の誤り (rtol = 0 かつ atol = 0, または rtol[i] = 0 かつ atol[i] = 0)
= -7: 入力パラメータ atol の誤り (atol < 0 または atol[i] < 0)
= -11: 入力パラメータ lwork の誤り. (lwork が小さ過ぎる)
= -13: 入力パラメータ liwork の誤り. (lwork が小さい)
= -14: 入力パラメータ info の誤り. (info != 0, 1)
= -19: 入力パラメータ irev の誤り. (irev != 1 〜 10, 50, 51)
= 1: 正常終了.
= 2: irtrn による中断 (正常終了)
= 11: 収束しなかった (ステップ数の最大値を超えたなど). |
| [out] | tt | irev = 1〜10: 再呼び出し時に2次微分係数値を求めてyypp[]に入れるべき点のtの値を返す.
irev = 50, 51: 現在のtの値. |
| [out] | yy[] | 配列 yy[lyy] (lyy >= n)
irev = 1〜10: 再呼び出し時に2次微分係数値を求めてyypp[]に入れるべき点のyの値を返す.
irev = 50, 51: 現在のtにおけるyの値. |
| [in,out] | yypp[] | 配列 yypp[lyypp] (lyypp >= n)
[in] irev = 1〜10: t(= tt)およびy(= yy[])における2次微分係数, すなわち, yypp[i] = d2yi/dt2 = fi(tt, yy[0], ..., yy[n-1]) (i = 0 〜 n-1) を計算して, 再呼び出し時に設定すること.
[out] irev = 50, 51: 現在のtにおけるypの値. |
| [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〜10: ttおよびyy[]における2次微分係数の計算値をyypp[]に設定すること. yypp[]以外の変数を変更してはならない.
= 50, 51: iout = 1 の場合, 中間結果を出力するためにステップが成功するごとにこのコードで戻る. t, y, yp の値が tt, yy[], yypp[] に返される. |
- 参考文献
- E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations I", Springer Series in Computational Mathematics, Springer-Verlag (1987)
- J. R. Dormand and P. J. Prince, "New Runge-Kutta algorithms for numerical simulation in dynamical astoronomy", Celestial Mechanics 18, 223?232 (1978)
|
1.9.6