|
|
◆ _odex2()
| void _odex2 |
( |
int |
n, |
|
|
void(*)(int, double, double *, double *) |
f, |
|
|
double * |
t, |
|
|
double |
y[], |
|
|
double |
yp[], |
|
|
double |
tout, |
|
|
double * |
rtol, |
|
|
double * |
atol, |
|
|
int |
itol, |
|
|
void(*)(int, double, double, double *, double *, int, double *, int *, int *) |
solout, |
|
|
int |
iout, |
|
|
double |
work[], |
|
|
int |
lwork, |
|
|
int |
iwork[], |
|
|
int |
liwork, |
|
|
double |
rcont[], |
|
|
int |
lrcont, |
|
|
int |
icont[], |
|
|
int |
licont, |
|
|
int * |
info |
|
) |
| |
常微分方程式の初期値問題 (補外法) (2階微分方程式用)
- 目的
- 本ルーチンは2階の常微分方程式の初期値問題
d2y/dt2 = f(t, y), ただし t = t0 において y = y0, y' = y'0
odex2はシュテルマー公式にもとづく補外法のプログラムである. ステップ制御, 次数選択および密出力機能を持っている.
詳細については下記参考文献を参照せよ.
- 引数
-
| [in] | n | 微分方程式の数. (n >= 1) |
| [in] | f | 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. void f(int n, double t, const double y[], double ypp[])
{
ypp[i] = tおよびy[]における2次微分係数の計算値 (i = 0 〜 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 >= 2n) (rtolまたはrtol[i] >= 0)
相対誤差許容値. y[i]およびyp[i]の局所誤差がおおむね次式以下になるようにする.
rtol*abs(y[i])+atol および rtol*abs(yp[i])+atol (i = 0〜n-1) (itol = 0 の場合)
または
rtol[i]*abs(y[i])+atol[i] および rtol[i+n]*abs(yp[i])+atol[i+n] (i = 0〜n-1) (itol = 1 の場合).
rtolとatol, または, rtol[i]とatol[i] (i = 0〜2n-1) が同時に0になってはいけない. |
| [in] | atol | スカラー(itol = 0 の場合) または 配列 atol[latol] (itol = 1 の場合) (latol >= 2n) (atolまたはatol[i] >= 0)
絶対誤差許容値. y[i]およびyp[i]の局所誤差がおおむね次式以下になるようにする.
rtol*abs(y[i])+atol および rtol*abs(yp[i])+atol(i = 0〜n-1) (itol = 0 の場合)
または
rtol[i]*abs(y[i])+atol[i] および rtol[i+n]*abs(yp[i])+atol[i+n](i = 0〜n-1) (itol = 1 の場合).
rtolとatol, または, rtol[i]とatol[i] (i = 0〜2n-1) が同時に0になってはいけない. |
| [in] | itol | rtolおよびatolがスカラーか配列かを指定する.
= 0: rtol および atolはスカラー.
= 1: rtol および atolは配列.
(他の値であれば itol = 0 とみなす.) |
| [in] | solout | 中間結果を出力するユーザー定義サブルーチンで, ステップが成功するごとに呼び出される. 次のように定義すること. void solout(int nr, double told, double t, double y[], double yp[], int n, double rcont[], int icont[], int *irtrn)
{
nr番目のステップtでのy[]およびyp[]の値が渡されるのでこれを出力する. toldは前回のtの値, nは方程式の次数である.
irtrnの入力値は, 初回, 中間, または, 最終呼び出しのとき, それぞれ 0, 1 または 2である. irtrnを使用して積分を中断することもできる. soloutでirtrnを負の値に設定すると積分を中断して info = 2 で終了する.
iout = 2 の場合には, 制御情報 rcont[] および icont[] により密出力を行うことができる. 区間[told, t]内の任意の点t2での解y[i] (0 <= i <= n-1)を次の関数呼び出しにより求めることができる.
y[i] = contx2(i, t2, rcont, icont);
}
|
| [in] | iout | soloutの呼び出し指定.
= 0: soloutの呼び出しを行わない
= 1: soloutの呼び出しを行う
= 2: 密出力のためにsoloutの呼び出しを行う. この場合, iwork[7]に密出力を行う要素数(nrdens)を設定しなければならない.
(他の値であれば iout = 0 とみなす.) |
| [in] | work[] | 配列 work[lwork]
作業領域.
work[0]〜work[19]はプログラムのためのパラメータ領域として使用される. 入力パラメータを0に設定すると, パラメータ値としてはそれぞれのデフォルト値を使用する.
[in]
work[0]: ステップ幅の初期値.
= 1/(f'のノルム) であり, 通常0.1〜0.001が適切. この値の選択はシビアではない. (デフォルト値 = 0.0001)
work[1]: ステップ幅の最大値 (デフォルト値 = tout - t)
work[3] および work[4]: ステップ幅選択パラメータ. j番目の対角エントリーについての新しいステップ幅は次式の範囲で選ばれる.
facmin/work[4] <= j番目のhnew/hold <= 1/facmin
ただし, facmin = work[3]^(1/(2*j - 1))
(デフォルト値: work[3] = 0.02, work[4] = 4)
work[5], work[6]: 次数選択パラメータ.
w[k-2] <= w[k-1]*work[5] : ステップ幅を減らす
w[k-1] <= w[k-2]*work[6] : ステップ幅を増やす
(デフォルト値: work[5] = 0.8, work[6] = 0.9)
work[7], work[8]: ステップ幅制御の安全係数.
hnew = h*work[8]*(work[7]*tol/err)^(1/(j - 1))
(デフォルト値: work[7] = 0.65, work[8] = 0.94)
work[9]: 安定性のチェック結果がNGの場合にステップ幅を減らすが, そのときに使う係数 (デフォルト値 = 0.5)
[out]
work[0]: 最後に採用されたステップのステップ幅 |
| [in] | lwork | 配列 work[]のサイズ. (lwork >= n*(2*km + 8) + 5*km + km*(2*km + 3)*nrdens + 20 (km: iwork[2]を参照, nrdens: iwork[7]を参照))
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[2]: km = 補外表の最大列数 (iwork[2] >= 3) (デフォルト値 = 9) (iwork[2] < 3 の場合, デフォルト値を使用する)
iwork[3]: ステップ数列の選択 (iout = 2 の場合 iwork[3] >= 4) (iwork[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, ...
(デフォルト値 = 1 (iout = 0または1), 4 (iout = 2))
iwork[6]: 補外式の次数を定める.
mu = 2*kappa - iwork[6] + 1 (1 <= iwork[6] <= 8)
(デフォルト値 = 6) (iwork[6]の値が範囲外の場合, デフォルト値を使用する)
iwork[7]: nrdens = 密出力を行う要素数. iwork[7] < 0 であれば, nrdens = abs(iwork[7]) とする. (デフォルト値 = n)
nrdensは iout = 2 の場合のみ有効である. その他の場合, nrdens = 0 とみなす.
0 < nrdens < n の場合, icont[] を設定しなければならない(下記参照).
iwork[8]: 密出力の式における誤差推定
= 0: 行う
= 1: 行わない
(デフォルト値 = 0) (iout = 2 の時に有効)
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 >= km + 20 (km: iwork[2]を参照))
liwork < 0 であれば abs(liwork) を使用し, iwork[0]〜iwork[19]を 0 に初期化する. |
| [out] | rcont[] | 配列 rcont[lrcont]
密出力の制御情報のための領域. iout = 2 の場合にのみ使用される. |
| [in] | lrcont | 配列 rcont[]のサイズ. (lrcont >= (2*km + 6)*nrdens (km: iwork[2]を参照, nrdens: iwork[7]を参照)) |
| [in,out] | icont[] | 配列 icont[licont]
密出力を行う要素番号のリスト. iout = 2 の場合にのみ使用される.
nrdens = n または iwork[7] < 0 の場合, 自動的に 1, ..., nrdens と設定される.
0 < nrdens < n の場合, 密出力を行う要素番号を設定しておかなければならない.
nrdensについてはiwork[7]を参照せよ. |
| [in] | licont | 配列 icont[]のサイズ. (licont >= nrdens (nrdens: iwork[7]を参照)) |
| [in,out] | info | [in]
= 0: 初期化して計算を開始する(新たな問題を解く).
= 1: toutだけを変えて計算を続ける(前回呼び出しの計算を再開する).
[out]
= -1: 入力パラメータ n の誤り (n < 1)
= -7: 入力パラメータ rtol の誤り (rtol < 0 または rtol[i] < 0)
= -7: 入力パラメータ rtol または atol の誤り (rtol = 0 かつ atol = 0, または rtol[i] = 0 かつ atol[i] = 0)
= -8: 入力パラメータ atol の誤り (atol < 0 または atol[i] < 0)
= -13: 入力パラメータ lwork の誤り (lworkが小さすぎる)
= -15: 入力パラメータ liwork の誤り (liworkが小さすぎる)
= -17: 入力パラメータ lrcont の誤り (lrcontが小さすぎる)
= -19: 入力パラメータ licont の誤り (licontが小さすぎる)
= -20: 入力パラメータ info の誤り (info != 0 かつ info != 1)
= 1: 正常終了
= 2: soloutによる中断 (正常終了)
= 11: 最大ステップ数を超えた |
- 出典
- 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)
- 目的
- 本ルーチンは2階の常微分方程式の初期値問題
d2y/dt2 = f(t, y), ただし t = t0 において y = y0, y' = y'0
odex2はシュテルマー公式にもとづく補外法のプログラムである. ステップ制御, 次数選択および密出力機能を持っている.
詳細については下記参考文献を参照せよ.
- 引数
-
| [in] | n | 微分方程式の数. (n >= 1) |
| [in] | f | 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. void f(int n, double t, const double y[], double ypp[])
{
ypp[i] = tおよびy[]における2次微分係数の計算値 (i = 0 〜 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 >= 2n) (rtolまたはrtol[i] >= 0)
相対誤差許容値. y[i]およびyp[i]の局所誤差がおおむね次式以下になるようにする.
rtol*abs(y[i])+atol および rtol*abs(yp[i])+atol (i = 0〜n-1) (itol = 0 の場合)
または
rtol[i]*abs(y[i])+atol[i] および rtol[i+n]*abs(yp[i])+atol[i+n] (i = 0〜n-1) (itol = 1 の場合).
rtolとatol, または, rtol[i]とatol[i] (i = 0〜2n-1) が同時に0になってはいけない. |
| [in] | atol | スカラー(itol = 0 の場合) または 配列 atol[latol] (itol = 1 の場合) (latol >= 2n) (atolまたはatol[i] >= 0)
絶対誤差許容値. y[i]およびyp[i]の局所誤差がおおむね次式以下になるようにする.
rtol*abs(y[i])+atol および rtol*abs(yp[i])+atol(i = 0〜n-1) (itol = 0 の場合)
または
rtol[i]*abs(y[i])+atol[i] および rtol[i+n]*abs(yp[i])+atol[i+n](i = 0〜n-1) (itol = 1 の場合).
rtolとatol, または, rtol[i]とatol[i] (i = 0〜2n-1) が同時に0になってはいけない. |
| [in] | itol | rtolおよびatolがスカラーか配列かを指定する.
= 0: rtol および atolはスカラー.
= 1: rtol および atolは配列.
(他の値であれば itol = 0 とみなす.) |
| [in] | solout | 中間結果を出力するユーザー定義サブルーチンで, ステップが成功するごとに呼び出される. 次のように定義すること. void solout(int nr, double told, double t, double y[], double yp[], int n, double rcont[], int icont[], int *irtrn)
{
nr番目のステップtでのy[]およびyp[]の値が渡されるのでこれを出力する. toldは前回のtの値, nは方程式の次数である.
irtrnの入力値は, 初回, 中間, または, 最終呼び出しのとき, それぞれ 0, 1 または 2である. irtrnを使用して積分を中断することもできる. soloutでirtrnを負の値に設定すると積分を中断して info = 2 で終了する.
iout = 2 の場合には, 制御情報 rcont[] および icont[] により密出力を行うことができる. 区間[told, t]内の任意の点t2での解y[i] (0 <= i <= n-1)を次の関数呼び出しにより求めることができる.
y[i] = contx2(i, t2, rcont, icont);
}
|
| [in] | iout | soloutの呼び出し指定.
= 0: soloutの呼び出しを行わない
= 1: soloutの呼び出しを行う
= 2: 密出力のためにsoloutの呼び出しを行う. この場合, iwork[7]に密出力を行う要素数(nrdens)を設定しなければならない.
(他の値であれば iout = 0 とみなす.) |
| [in] | work[] | 配列 work[lwork]
作業領域.
work[0]〜work[19]はプログラムのためのパラメータ領域として使用される. 入力パラメータを0に設定すると, パラメータ値としてはそれぞれのデフォルト値を使用する.
[in]
work[0]: ステップ幅の初期値.
= 1/(f'のノルム) であり, 通常0.1〜0.001が適切. この値の選択はシビアではない. (デフォルト値 = 0.0001)
work[1]: ステップ幅の最大値 (デフォルト値 = tout - t)
work[3] および work[4]: ステップ幅選択パラメータ. j番目の対角エントリーについての新しいステップ幅は次式の範囲で選ばれる.
facmin/work[4] <= j番目のhnew/hold <= 1/facmin
ただし, facmin = work[3]^(1/(2*j - 1))
(デフォルト値: work[3] = 0.02, work[4] = 4)
work[5], work[6]: 次数選択パラメータ.
w[k-2] <= w[k-1]*work[5] : ステップ幅を減らす
w[k-1] <= w[k-2]*work[6] : ステップ幅を増やす
(デフォルト値: work[5] = 0.8, work[6] = 0.9)
work[7], work[8]: ステップ幅制御の安全係数.
hnew = h*work[8]*(work[7]*tol/err)^(1/(j - 1))
(デフォルト値: work[7] = 0.65, work[8] = 0.94)
work[9]: 安定性のチェック結果がNGの場合にステップ幅を減らすが, そのときに使う係数 (デフォルト値 = 0.5)
[out]
work[0]: 最後に採用されたステップのステップ幅 |
| [in] | lwork | 配列 work[]のサイズ. (lwork >= n*(2*km + 8) + 5*km + km*(2*km + 3)*nrdens + 20 (km: iwork[2]を参照, nrdens: iwork[7]を参照))
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[2]: km = 補外表の最大列数 (iwork[2] >= 3) (デフォルト値 = 9) (iwork[2] < 3 の場合, デフォルト値を使用する)
iwork[3]: ステップ数列の選択 (iout = 2 の場合 iwork[3] >= 4) (iwork[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, ...
(デフォルト値 = 1 (iout = 0または1), 4 (iout = 2))
iwork[6]: 補外式の次数を定める.
mu = 2*kappa - iwork[6] + 1 (1 <= iwork[6] <= 8)
(デフォルト値 = 6) (iwork[6]の値が範囲外の場合, デフォルト値を使用する)
iwork[7]: nrdens = 密出力を行う要素数. iwork[7] < 0 であれば, nrdens = abs(iwork[7]) とする. (デフォルト値 = n)
nrdensは iout = 2 の場合のみ有効である. その他の場合, nrdens = 0 とみなす.
0 < nrdens < n の場合, icont[] を設定しなければならない(下記参照).
iwork[8]: 密出力の式における誤差推定
= 0: 行う
= 1: 行わない
(デフォルト値 = 0) (iout = 2 の時に有効)
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 >= km + 20 (km: iwork[2]を参照))
liwork < 0 であれば abs(liwork) を使用し, iwork[0]〜iwork[19]を 0 に初期化する. |
| [out] | rcont[] | 配列 rcont[lrcont]
密出力の制御情報のための領域. iout = 2 の場合にのみ使用される. |
| [in] | lrcont | 配列 rcont[]のサイズ. (lrcont >= (2*km + 6)*nrdens (km: iwork[2]を参照, nrdens: iwork[7]を参照)) |
| [in,out] | icont[] | 配列 icont[licont]
密出力を行う要素番号のリスト. iout = 2 の場合にのみ使用される.
nrdens = n または iwork[7] < 0 の場合, 自動的に 1, ..., nrdens と設定される.
0 < nrdens < n の場合, 密出力を行う要素番号を設定しておかなければならない.
nrdensについてはiwork[7]を参照せよ. |
| [in] | licont | 配列 icont[]のサイズ. (licont >= nrdens (nrdens: iwork[7]を参照)) |
| [in,out] | info | [in]
= 0: 初期化して計算を開始する(新たな問題を解く).
= 1: toutだけを変えて計算を続ける(前回呼び出しの計算を再開する).
[out]
= -1: 入力パラメータ n の誤り (n < 1)
= -7: 入力パラメータ rtol の誤り (rtol < 0 または rtol[i] < 0)
= -7: 入力パラメータ rtol または atol の誤り (rtol = 0 かつ atol = 0, または rtol[i] = 0 かつ atol[i] = 0)
= -8: 入力パラメータ atol の誤り (atol < 0 または atol[i] < 0)
= -13: 入力パラメータ lwork の誤り (lworkが小さすぎる)
= -15: 入力パラメータ liwork の誤り (liworkが小さすぎる)
= -17: 入力パラメータ lrcont の誤り (lrcontが小さすぎる)
= -19: 入力パラメータ licont の誤り (licontが小さすぎる)
= -20: 入力パラメータ info の誤り (info != 0 かつ info != 1)
= 1: 正常終了
= 2: soloutによる中断 (正常終了)
= 11: 最大ステップ数を超えた |
- 出典
- 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)
|
1.9.6