debdf debdf_r radaua radaua_r rodasa rodasa_r seulexa seulexa_r

debdf_r() void debdf_r ( int  n, double *  t, double  y[], double  tout, double *  rtol, double *  atol, int  itol, int  mode, int  itstop, int  idjac, double  work[], int  lwork, int  iwork[], int  liwork, int *  info, double *  tt, double  yy[], double  yyp[], int  ldyypd, double  yypd[], int *  irev  ) 常微分方程式の初期値問題 (1〜5次 後退微分公式 (BDF)) (リバースコミュニケーション版) 目的 本ルーチンは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. 本ルーチンは, 微分方程式ソルバー・パッケージDEPACに収録されているプログラムのひとつである(DEPACは, derkf, deabmおよびdebdfからなる). debdfは可変次数(1〜5次)後退微分公式のプログラムである. 3つの中で最も複雑である. 主に, スティフな微分方程式の低〜中精度な解を求めるために設計されている. 問題が強いスティフ性を持つ場合, derkfおよびdeabmはdebdfに比べて極めて非効率である. しかし, 非スティフな問題に対しては, debdfはderkfおよびdeabmに比べて非効率である. これは, より多くのメモリを必要とし, よりオーバーヘッドが大きく, 次数が低いため効率的に精度を上げることができないためである. debdfは A. C. Hindmarsh によるプログラム LSODE の修正版のドライバー・ルーチンである. debdf_rはdebdfのリバースコミュニケーション版である. 引数 [in] n 微分方程式の数. (n >= 1) [in,out] t 本ルーチンはtからtoutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. 新たなtoutにおける解を求めるために積分を継続することができる(インターバル・モードとよぶ). その場合, 2回目以降の呼び出しでは, tは前回のtoutに等しくなければならない. また, toutまでの各中間ステップにおいて戻ることもできる(中間結果出力モードとよぶ). 本モードは解の挙動をみたい場合に使うとよい. モードはパラメータmodeで指定できる. [in] 独立変数tの初期値. [out] 独立変数tの最終ステップの最後の点の値 (インターバル・モードでは通常toutに等しい). 解がこの点まで正常に求められたことを示す. [in,out] y[] 配列 y[ly] (ly >= n) [in] tの初期値における従属変数y[]の初期値. [out] 最終のt(インターバル・モードにおいてはtoutに等しい)において求められた解(の近似値). [in] tout 解を求めたい点を設定する. tout = t とすることができ, その場合はtにおける微分係数を計算して戻る. 積分を行うのはtについて前進方向(tout > t)でも後退方向(tout < t)でもよい. ただし, 初期化時でなければ積分方向を変更することはできない. 要求精度を満たすように自動的に選ばれたステップ幅を用いてtからtoutに向かって解が求められる. 必要であれば, 各中間ステップにおける解とその微分係数を確認するために戻ることができる(中間結果出力モード). ただし, その場合であっても本来のとおりtoutを指定しなければならない. [in,out] rtol スカラー(itol = 0 の場合) または 配列 rtol[lrtol] (itol = 1 の場合) (lrtol >= n) (rtolまたはrtol[i] >= 0) 求める解の精度を指定する相対誤差許容値. 本パラメータはitolに従ってスカラーまたは配列のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる. 許容値は各ステップにおける局所誤差テストに使われ, y[]の各要素がおおむね次式を満たすようにする (i = 0 〜 n-1).   abs(局所誤差) <= rtol*abs(y[i]) + atol (itol = 0 の場合)     または   abs(局所誤差) <= rtol[i]*abs(y[i]) + atol[i] (itol = 1 の場合) rtol = 0 と設定するとその要素については純粋に絶対誤差テストとなる. [in] 誤差テストのための相対誤差許容値. [out] info = 4 で戻ったときには計算を続けるために適当な大きさに調節した値を返す. それ以外の場合には変更されない. [in,out] atol スカラー(itol = 0 の場合) または 配列 atol[latol] (itol = 1 の場合) (latol >= n) (atolまたはatol[i] >= 0) 求める解の精度を指定する絶対誤差許容値. 本パラメータはitolに従ってスカラーまたは配列のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる. 許容値は各ステップにおける局所誤差テストに使われ, y[]の各要素がおおむね次式を満たすようにする (i = 0 〜 n-1).   abs(局所誤差) <= rtol*abs(y[i]) + atol (itol = 0 の場合)     または   abs(局所誤差) <= rtol[i]*abs(y[i]) + atol[i] (itol = 1 の場合) atol = 0 と設定するとその要素については純粋に相対誤差テストとなる. [in] 誤差テストのための絶対誤差許容値. [out] info = 4 で戻ったときには計算を続けるために適当な大きさに調節した値を返す. それ以外の場合には変更されない. [in] itol rtolおよびatolがスカラーか配列かを指定する. 現問題について初回呼び出し後にこれを変更することはできない. = 0: rtol および atolはスカラー. = 1: rtol および atolは配列. (他の値であれば itol = 0 とみなす.) [in] mode 動作モード. = 0: toutで戻る (インターバル・モード). = 1: ステップごとに戻る (中間結果出力モード). (他の値であれば mode = 0 とみなす.) [in] itstop 本ルーチンはtoutを超えて積分を行い内挿によりtoutにおける値を求めることがあるが, ある点において不連続になったり, その点を超えて関数あるいは微分が定義されていなかったりして, その点を超えて積分を行うのが許されないことがある. そのような場合, ストップポイント tstop を設定することができ, その点を超えての積分を行わないようにすることができる. = 0: ストップポイントを設定しない = 1: work[0] = tstop をストップポイントとして設定する (他の値であれば itstop = 0 とみなす.) [in] idjac ヤコビ行列をフル(通常)行列として計算するか帯行列として計算するかを指定する. 一般的には, 2*ml + mu < n であれば帯行列として計算した方が計算時間および記憶容量の点で有利である. ここで, mlは下帯幅, muは上帯幅である. = 0: n x n フル(通常)行列. = 1: 帯行列. 帯幅を iwork[0] = ml, iwork[1] = mu と設定しておくこと. (他の値であれば idjac = 0 とみなす.) [out] work[] 配列 work[lwork] 作業領域. [in] lwork 配列 work[]のサイズ. (lwork >= n^2 + 10*n + 24 (idjac = 0 (ヤコビ行列がフル行列の場合)), (2*ml + mu + 11)*n + 24 (idjac = 1 (ヤコビ行列が帯行列の場合))) [out] iwork[] 配列 iwork[liwork] 作業領域. [in] liwork 配列 iwork[]のサイズ. (liwork >= n + 38) [out] info [in] 制御コード. = 0: 問題の最初の呼び出し時(新たに問題を開始する場合)には info = 0 と設定せよ. 初期化が行われ, 新たな問題についての計算が開始される. = 1, 2, 11〜15: info = 1, 2, 11〜15 で戻った場合, 計算を継続するために, infoの値を変えずに再呼び出しすることができる(下記説明を参照せよ). [out] リターンコード. ユーザーはこの値をチェックし, info = 1, 2, 11〜15 であれば必要により次のアクションとして本ルーチンを再度呼び出すことができる. = -1: 入力パラメータ n の誤り (n < 1) = -2: 入力パラメータ t の誤り (継続呼び出しで t = tout または t != 前回のtout) = -4: 入力パラメータ tout の誤り (積分の方向を変更した, または, toutがtstopと矛盾する) = -5: 入力パラメータ rtol の誤り (rtol < 0 または rtol[i] < 0) = -6: 入力パラメータ atol の誤り (atol < 0 または atol[i] < 0) = -12: 入力パラメータ lwork の誤り (lworkが小さすぎる) = -13: 入力パラメータ iwork の誤り (ml < 0 または ml >= n, または mu < 0 または mu >= n) = -14: 入力パラメータ liwork の誤り (liworkが小さすぎる) = -15: 入力パラメータ info の誤り (info != 0, 1, 2, 11〜15) = -19: 入力パラメータ ldyypd の誤り (ldyypdが小さすぎる) = 1: 正常終了 (t = tout). toutを再設定して再呼び出しすることができる. toutはtと異なる値であること. = 2: 中間結果出力モードのため戻った (まだ toutに達していない). tout方向の次のステップを続行するために再呼び出しすることができる. = 11: 最大ステップ数(10000)を超えた. 再呼び出しして続行することができる. さらに10000ステップ続けることができる. = 12: 誤差許容値が厳しすぎる. 再呼び出しして, 自動的に調節された許容値を使って続行することができる. 許容値を自分で変更してから再呼び出しすることもできる. = 13: 求められた解が0であるため, 純粋な相対誤差テスト(atol = 0)ができない. atolを正の値にして再呼び出しして続行することができる. = 14: 最後の試行ステップにおいて収束テストの失敗を繰り返した. infoを変えずに再呼び出しして続行することができるが, ヤコビ行列が正しくないことが原因かもしれない. = 15: 最後の試行ステップにおいて誤差テストの失敗を繰り返した. infoを変えずに再呼び出しして続行することができるが, 解に特異性があることが原因かもしれない. = 16: 無限ループが検出された. [out] tt irev = 1〜11: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のtの値を返す. [out] yy[] 配列 yy[lyy] (lyy >= n) irev = 1〜11: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のyの値を返す. [in] yyp[] 配列 yyp[lyyp] (lyyp >= n) irev = 1〜11: t(= tt)およびy(= yy[])における微分係数, すなわち, yyp[i] = dyi/dt = fi(tt, yy[0], ..., yy[n-1]) (i = 0 〜 n-1) を計算して, 再呼び出し時に返す. [in] ldyypd = 0: ヤコビ行列は数値計算(差分近似)によりルーチン内で計算される. irev = 30または31で戻ることはない. yypd[][]は参照されない. != 0: 二次元配列yypd[][]の整合寸法. (ldyypd >= n (idjac = 0 の場合), ldyypd >= ml + mu + 1 (idjac = 1 の場合))   irev = 30または31で戻ったとき, ユーザーはヤコビ行列を解析的に計算しyypd[][]に設定すること. [in,out] yypd[][] 配列 yypd[lyypd][ldyypd] (lyypd >= n) irev = 30, 31: 通常行列形式(idjac = 0 の場合)または帯行列形式(idjac = 1 の場合)で, 二次元のヤコビ行列を格納する. 帯行列形式の詳細は下記を参照のこと. yypd[][]の全ての要素はirev = 30または31で戻る前に0に初期化されるので, 0でない要素だけ設定すればよい. [in,out] irev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行いirevを変更せずに再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはinfoをチェックすること. = 1〜11: ttおよびyy[]における微分係数の計算値をyyp[]に設定すること. yyp[]以外の変数を変更してはならない. = 30, 31: ttおよびyy[]におけるヤコビ行列の計算値をyypd[]に設定する. yypd[]以外の変数を変更してはならない. 詳細 次の例は, n = 6, ml = 2, mu = 1 の場合の帯行列形式を表す. 一般行列形式: a11 a12 0 0 0 0 a21 a22 a23 0 0 0 a31 a32 a33 a34 0 0 0 a42 a43 a44 a45 0 0 0 a53 a54 a55 a56 0 0 0 a64 a65 a66 帯行列形式: * a12 a23 a34 a45 a56 a11 a22 a33 a44 a55 a66 a21 a32 a43 a54 a65 * a31 a42 a53 a64 * * *で示された配列要素は使用されない. 出典 SLATEC (DEPAC)