|
|
◆ _dverk_r()
| void _dverk_r |
( |
int |
n, |
|
|
double * |
t, |
|
|
double |
y[], |
|
|
double |
tout, |
|
|
double |
tol, |
|
|
double |
c[], |
|
|
int |
lc, |
|
|
double |
work[], |
|
|
int |
lwork, |
|
|
int * |
info, |
|
|
double * |
tt, |
|
|
double |
yy[], |
|
|
double |
yyp[], |
|
|
int * |
irev |
|
) |
| |
常微分方程式の初期値問題 (6(5)次 ルンゲ・クッタ・ヴァーナー法) (リバースコミュニケーション版)
- 目的
- 本ルーチンは1階の常微分方程式の初期値問題
dy/dt = f(t, y), ただし t = t0 において y = y0
本ルーチンは6(5)次ルンゲ・クッタ・ヴァーナー法のプログラムである. 微分係数の計算に時間がかからない非スティフな問題に適している.
dverk_rはdverkのリバースコミュニケーション版である.
- 引数
-
| [in] | n | 微分方程式の数. (n >= 1) |
| [in,out] | t | 本ルーチンはtからtoutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される.
[in] 独立変数tの初期値.
[out] 独立変数tの最終ステップの最後の点の値(通常toutに等しい). 解がこの点まで正常に求められたことを示す. toutの値を変えて本ルーチンを再呼び出しすることにより新たな点まで積分を継続することができる. |
| [in,out] | y[] | 配列 y[ly] (ly >= n)
[in] tの初期値における従属変数y[]の初期値.
[out] 最終のt(通常toutに等しい)において求められた解(の近似値). |
| [in] | tout | 解を求めたい点を設定する. 積分を行うのはtについて前進方向(tout > t)でも後退方向(tout < t)でもよい.
要求精度を満たすように自動的に選ばれたステップ幅を用いてtからtoutに向かって解が求められる. |
| [in] | tol | 誤差許容値.
本ルーチンは, 大域誤差がtolに比例するように, 局所誤差のノルムを制御する. ノルムは, c[] パラメータで指定される誤差制御方式に依存する重み付きの, 最大ノルムである. デフォルトの重みは, k番目の要素について, 1/max(1, abs(y(k))) である. これは絶対および相対誤差制御を混合したものになる. |
| [in,out] | c[] | 配列 c[lc]
Array c[lc]
[in]
入力パラメータ (デフォルト動作をさせるにはすべて 0 に設定せよ).
c[0]: エラー制御. 局所誤差のノルムは重み付き推定誤差ベクトルの最大ノルムを使う. c[0]はその重みを指定する.
= 0: 重み = 1/max(1, abs(y(k))) (絶対および相対誤差制御を混合したもの)
= 1: 重み = 1 (絶対誤差制御)
= 2: 重み = 1/abs(y(k)) (相対誤差制御)
= 3: 重み = 1/max(abs(c[1]), abs(y(k))) (abs(y(k))が下限値abs(c[1])より小さくない限り相対誤差制御)
= 4(*): 重み = 1/max(abs(c[k+30]), abs(y[k])) (個別の下限値が使われる)
= 5(*): 重み = 1/abs(c[k+30])
(他の値であれば c[0] = 0 とみなす.)
(*) 4または5の場合, 配列c[]の大きさをn+30以上とし, c[30]〜c[n+29]を適切な値に設定しておかなければならない.
c[1]: 下限値. c[0] = 3 の場合に使われる.
c[2]: 最小ステップ幅 (hmin). c[2] = 0の場合, デフォルト値 = 10*max(dwarf, rreb*max((yの重み付きノルム)/tol, abs(x))) を使う. ただし, dwarfは非常に小さな正数 (= 1.0e-50), rrebは相対丸め誤差限界値である.
c[3]: ステップ幅の初期値 (hstart) = abs(c[3]). ただし, hminおよびhmaxの制限がある. c[3] = 0の場合, デフォルト値 = hmax*tol^(1/6) を使う.
c[4]: 問題のスケール尺度 (scale). c[4] = 0の場合, デフォルト値 = 1 を使う. 本パラメータは, hmaxの決定および判定基準の変更のために使われる. scaleの値が大きければより高精度の結果を求めようとする.
c[5] 最大ステップ幅 (hmax).
c[5] = 0 の場合:
c[4] = 0 の場合: hmax = 2
c[4] != 0 の場合: hmax = 2/abs(c[4])
c[5] != 0 の場合:
c[4] = 0 の場合: hmax = abs(c[5])
c[4] != 0 の場合: hmax = min(abs(c[5]), 2/abs(c[4]))
c[6]: 関数評価回数の最大値. c[6] = 0の場合, 評価回数のチェックを行わない. c[6] < 0の場合, abs(c[6])を使う.
c[7]: 割り込み 1.
= 0: 割り込みを行わない.
!= 0: ステップ幅の暫定値を選んだ後, 次のステップのためのhtrialおよびttrialを選ぶ直前に割り込みを行う. info = 2で戻り, 次にinfo = 2で再呼び出しすると割り込み時点から計算を再開する.
c[8]: 割り込み 2.
= 0: 割り込みを行わない.
!= 0: 最新の試行ステップの結果を採用するかどうかを決めた直後に割り込みを行う. 採用する場合にはinfo = 5, 採用しない場合にはinfo = 6で戻る. info = 5あるいは6で再呼び出しすると割り込み時点から計算を再開する. この場合, 必要ならばinfoを変更してもよい. 例えば, 本来採用されないステップを強制的に採用させる, あるいはその逆のことができる.
= 2: 上記に加えて密出力を可能にする. すなわち, ind = 5 (最後のステップでは ind = 3) のときには dverk_int ルーチンを使って直近のステップの区間内で補間値を求めることができる. 採用されたステップあたり1回の追加関数評価が必要になる.
c[30], c[31], ... c[n+29]: 下限値 (c[0] = 4 または 5 の場合).
[out]
出力パラメータ.
c[9]: 相対丸め誤差限界 (rreb)
c[10]: 非常に小さな正数 (dwarf)
c[11]: yの重み付きノルム
c[12]: 最小ステップ幅 (hmin)
c[13]: ステップ幅 (hmag)
c[14]: スケール尺度 (scale)
c[15]: 最大ステップ幅 (hmax)
c[16]: 試行ステップのtの値 (ttrial)
c[17]: 試行ステップのhの値 (htrial)
c[18]: 推定誤差
c[19]: 前回のtout
c[20]: toutに到達したことを示すフラグ
c[21]: 成功したステップ数
c[22]: 連続して失敗したステップ数
c[23]: 関数評価回数 |
| [in] | lc | 配列 c[]のサイズ. (lc >= 24, ただし c[0] = 4 または 5 のときは lc >= n + 30) |
| [out] | work[] | 配列 work[lwork]
作業領域. |
| [in] | lwork | 配列 work[]のサイズ. (lwork >= 9*n (ただし, c[8] = 2 の場合は lwork >= 12*n)) |
| [in,out] | info | [in] 制御コード.
最初の呼び出し時には info = 1 または 2 と設定すること.
= 1: すべてのパラメータを自動的にデフォルト値に初期化してから積分を開始する.
= 2: パラメータの初期化を行わずに積分を開始する. ユーザーは呼び出し前にあらかじめ入力パラメータ(c[0]〜c[8])を設定しておかなけらばならない.
= 3: 正常終了後の新toutによる再呼び出し.
= 4: 割り込み 1 からの再開 (c[7]を参照せよ).
= 5: 割り込み 2 (ステップ採用時) からの再開 (c[8]を参照せよ).
= 6: 割り込み 2 (ステップ不採用時) からの再開 (c[8]を参照せよ).
[out] リターンコード. ユーザーはこの値をチェックし, info = 3〜6 であれば必要により次のアクションとして本ルーチンを再度呼び出すことができる.
= -1: 入力パラメータ n の誤り (n < 1)
= -2: 入力パラメータ t の誤り (継続呼び出しで t = tout または t != 前回のtout)
= -5: 入力パラメータ tol の誤り (tol <= 0)
= -7: 入力パラメータ lc の誤り (lcが小さすぎる)
= -9: 入力パラメータ lwork の誤り (lworkが小さすぎる)
= -10: 入力パラメータ info の誤り (info != 1〜6)
= 3: 正常終了 (t = tout). toutを再設定して再呼び出しすることができる. 他の変数は変更してはならない.
= 4: 割り込み 1
= 5: 割り込み 2 (ステップを採用)
= 6: 割り込み 2 (ステップを不採用)
= 11: 関数評価回数の最大値を超えた
= 12: hminの値がhmaxより大きい (tolが小さすぎる可能性がある)
= 13: 特定ステップサイズがhminに等しいか小さく, 要求精度を満たすことができなかった (tolが小さすぎる可能性がある) |
| [out] | tt | irev = 1〜10: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のtの値を返す. |
| [out] | yy[] | 配列 yy[lyy] (lyy >= n)
irev = 1〜10: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のyの値を返す. |
| [in] | yyp[] | 配列 yyp[lyyp] (lyyp >= n)
irev = 1〜10: t(= tt)およびy(= yy[])における微分係数, すなわち, yyp[i] = dyi/dt = fi(tt, yy[0], ..., yy[n-1]) (i = 0 〜 n-1) を計算して, 再呼び出し時に返す. |
| [in,out] | irev | リバースコミュニケーションの制御変数.
[in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない.
[out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと.
= 0: 処理終了. 正常終了かどうかはinfoをチェックすること.
= 1〜10: ttおよびyy[]における微分係数の計算値をyyp[]に設定すること. yyp[]以外の変数を変更してはならない. |
- 出典
- netlib/ode
|
1.9.6