debdf debdf_r radaua radaua_r rodasa rodasa_r seulexa seulexa_r

rodasa_r() void rodasa_r ( int  n, int  ifcn, double *  t, double  y[], double  tout, double  tend, double  rtol[], double  atol[], int  itol, int  ijac, int  mljac, int  mujac, int  idfx, int  imas, int  mlmas, int  mumas, int  mode, double  work[], int  lwork, int  iwork[], int  liwork, int *  info, double *  tt, double  yy[], double  yyp[], int  ldyypd, double  yypd[], int *  irev  ) 常微分方程式の初期値問題 (4(3)次 ローゼンブロック法) (リバースコミュニケーション版) 目的 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式) の初期値問題 M * dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, y は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. M は質量マトリックスである. また, t0 および y0 はそれぞれ t および y の既知の初期値である. 本ルーチンは, 4(3)次ローゼンブロック法 プログラム RODAS (文献 (1)) を書き直したものである. 主として密出力機能の使用法を変更した. 本ルーチンは rodasa のリバースコミュニケーション版である. 引数 [in] n 微分方程式の数. (n >= 1) [in] ifcn f(t, y) が t に依存するかどうかを設定する. = 0: f(t, y) は t に依存しない (自励系). = 1: f(t, y) は t に依存する (非自励系). [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] 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[] の各要素が次式を満たすように自動的に選ばれたステップ幅を用いて積分が行われる. (局所誤差) <= rtol*abs(y[i]) + atol (itol = 0 の場合) (局所誤差) <= rtol[i]*abs(y[i]) + atol[i] (itol = 1 の場合(i = 0 〜 n - 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] ijac ヤコビ行列の計算方法を指定する. = 0: ヤコビ行列を有限差分近似により求める. irev = 30 で戻ることはない. = 1: ヤコビ行列を irev = 30 によりユーザーが求める. [in] mljac ヤコビ行列の下帯幅. (0 <= mljac <= n) mljac = n の場合, ヤコビ行列はフル行列形式で格納される. mljac < n の場合, ヤコビ行列は帯行列形式で格納される. [in] mujac ヤコビ行列の上帯幅. (0 <= mujac <= n) mljac = n の場合, mujac は無視される. [in] idfx f(t, y) の t についての偏微分の計算方法を指定する. = 0: 偏微分を有限差分により求める. irev = 35 で戻ることはない. = 1: ifcn = 1 であれば偏微分を irev = 35 によりユーザーが求める. ifcn = 0 の場合, irev = 30 で戻ることはない. [in] imas 質量マトリクス M が単位行列であるかどうかを指定する. = 0: M は単位行列である. irev = 40 で戻ることはない. = 1: M は irev = 40 によりユーザーが設定する. [in] mlmas 質量マトリクス M の下帯幅. (0 <= mljac <= n) mlmas = n の場合, M はフル行列形式で格納される. mlmas < n の場合, M は帯行列形式で格納される. [in] mumas 質量マトリクス M の上帯幅. (0 <= mumas <= n) mlmas = n の場合, mumas は無視される. [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 と同じステップを踏む. [in,out] work[] 配列 work[lwork] 作業領域. ただし, work[0]〜work[19] はプログラムのためのパラメータ領域として使用される. info = 0 の呼び出し時に 0 に設定されている入力パラメータはそのデフォルト値が使用される. [in] work[0]: ステップ幅の初期値. (デフォルト値 = 1.0e-4) work[2]: ステップ幅の最大値 (デフォルト値 = abs(tend - t)) work[4], work[5]: ステップ幅選択パラメータ. (work[4] <= 1, work[5] >= 1) (デフォルト値: work[4] = 0.2, work[5] = 6) work[4] <= hnew/hold <= work[5] となるようにステップ幅が選ばれる. work[8]: ステップ幅推定時の安全係数. (0.001 < work[8] < 1) (デフォルト値 = 0.9) [out] work[1]: 最後に使用されたステップ幅. [in] lwork 配列 work[]のサイズ. (abs(lwork) >= n*(ldjac + 15) + nm1*(ldmas + lde) + 80) ただし, ldjac = nm1 (フル行列の場合) または mljac + mujac + 1 (帯行列の場合), lde = nm1 (フル行列の場合) または 2*mljac + mujac + 1 (帯行列の場合), ldmas = nm1 (フル行列の場合) または mlmas + mumas + 1 (帯行列の場合) または 1 (imas = 0 の場合) nm1 = n - m1 (m1 = iwork[8]) lwork < 0 であればその絶対値を採用し, work[0]〜work[19] をすべて 0 に設定する. [in,out] iwork[] 配列 iwork[liwork] 整数作業領域. ただし, iwork[0]〜iwork[19] はプログラムのためのパラメータ領域として使用される. info = 0 の呼び出し時に 0 に設定されている入力パラメータはそのデフォルト値が使用される. [in] iwork[0]: 係数を指定する. (デフォルト値 = 0) = 0: 文献(1)のP452(邦訳 P426)の方法. = 1: 同じ方法で異なる係数を使用. = 2: Gerd Steinebachの係数を使用した方法. iwork[1]: 許される最大ステップ数. (デフォルト値 = 100000) iwork[7]: ステップ幅戦略を指定する. (デフォルト値 = 1) = 1: モデル予期コントローラ (グスタフソン). = 2: 伝統的なステップ幅コントロール. iwork[8], iwork[9]: 連立微分方程式が次のような特殊な構造であるとする. y(i)' = y(i + m2) (i = 1, ..., m1) ただし, m1 は m2 の倍数であり, かつ残りの方程式が y'(m1), ..., y'(n - 1) に陽に依存しないものとする. このような場合, iwork[8] = m1 (> 0) および iwork[9] = m2 (> 0) (m1 + m2 <= n) と設定することにより計算時間を大幅に短縮することができる. (デフォルト値: iwork[8] = 0, iwork[9] = iwork[8]) iwork[12]: カウンタ iwork[13]〜iwork[19] のリセット方法を指定する. (デフォルト値 = 0) = 0: info = 0 のときにリセットする. = 1: info = 0 であってもリセットしない. [out] 統計情報(カウンタ). iwork[13]: 関数評価回数. iwork[14]: ヤコビ行列評価回数. (有限差分の場合を含む) iwork[15]: 全ステップ数. iwork[16]: 誤差評価により採用されたステップ数. iwork[17]: 誤差評価により不採用だったステップ数. iwork[18]: LU分解の回数. iwork[19]: 解の計算の回数. (ステップ幅選択のためのものはカウントされない) [in] liwork 配列 iwork[]のサイズ. (abs(liwork) >= nm1 + 60, ただし nm1 = n - m1 (m1 = iwork[8])) liwork < 0 であればその絶対値を採用し, 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: (エラー) ステップ幅が小さくなりすぎたため計算が継続できない. = 13: (エラー) 行列が繰り返し特異になった. [out] tt irev = 1〜9, 35 の場合: yyp[]を参照のこと. irev = 30 の場合: yypd[][]を参照のこと. [out] yy[] 配列 yy[lyy] (lyy >= n) irev = 1〜9, 35 の場合: yyp[]を参照のこと. irev = 30 の場合: yypd[][]を参照のこと. [in] yyp[] 配列 yyp[lyyp] (lyyp >= n) irev = 1〜9 の場合: t = tt および y = yy[] における微分値 dy/dt = f(t, y) を計算して, 再呼び出し時に yyp[] に設定しておくこと. irev = 35 の場合: t = tt および y = yy[] における f(t, y) の t についての偏微分を計算して, 再呼び出し時に yyp[] に設定しておくこと. [in] ldyypd irev = 30, 40 の場合: 2次元配列 yypd[][] の整合寸法. (ldyypd >= max(ljac, lmas), ljac および lmas については lwork を参照せよ) [in] yypd[][] 配列 yypd[lyypd][ldyypd] (lyypd >= n) irev = 30 の場合: t = tt および y = yy[] におけるヤコビ行列を計算して, 再呼び出し時に yypd[][] に設定しておくこと. mljac = n の場合はフル行列形式, mljac < n の場合は帯行列形式で格納すること. irev = 40 の場合: 再呼び出し時に質量マトリクスを yypd[][] に設定しておくこと. mlmas = n の場合はフル行列形式, mlmas < n の場合は帯行列形式で格納すること. [in,out] irev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行い irev を変更せずに再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかは info をチェックすること. = 1〜9: t = tt および y = yy[] における微分の計算値を yyp[] に設定する. yyp[] 以外の変数を変更してはならない. = 30: t = tt および y = yy[] におけるヤコビ行列の計算値を yypd[][] に設定する. yypd[][] 以外の変数を変更してはならない. = 35: t = tt および y = yy[] における f(t, y) の t についての偏微分の計算値を yyp[] に設定する. yyp[] 以外の変数を変更してはならない. = 40: 質量マトリックスを 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 * * *で示された配列要素は使用されない. 詳細 特殊構造を持つ問題 (iwork[8] および iwork[9] の設定) の例. 2階の連立方程式 p' = v v' = g(p, v) の場合, m1 = m2 = n/2 とできる. ただし, p および v は長さ n/2 のベクトルである. m1 > 0 の場合, ヤコビ行列と質量マトリクスを以下のように格納する必要がある. ・ヤコビ行列は, 非ゼロ要素 (第 m1 〜 n-1行) のみ格納する. ・mljac = n - m1 の場合, ヤコビ行列は通常行列である. dfi/dyj = d f(i+m1) / d y(j) (i = 1 〜 n-m1, j = 1 〜 n) ・0 <= mljac < n - m1 の場合, ヤコビ行列は帯行列である. m1 = m2 * mm とすると, mm + 1 個の部分行列のみを格納する. df(i - j + mujac + 1)/dy(j + k*m2) = d f(i + m1) / d y(j + k*m2) (i = 1〜mljac + mujac + 1, j = 1 〜 m2, k = 0 〜 mm) mljac はこれら mm + 1 個の部分行列の最大下帯幅である. mujac は mm + 1 個の部分行列の最大上帯幅であり, mljac = n - m1 の場合, 無視される. ・質量マトリクス M は, 単位行列ではない右下ブロックの n - m1 × n - m1 要素のみを格納する. ・mlmas = n - m1 の場合, この部分行列は通常行列である. am[j - 1][i - 1] = M(i + m1, j + m1) (i = 1 〜 n - m1, j = 1 〜 n - m1). ・0 <= mlmas < n - m1 の場合, この部分行列は帯行列である. am[j - 1][i - j + mumas] = M(i + m1, j + m1) (i = 1 〜 n - m1, j = 1〜mlmas + mumas + 1). mlmasは部分行列の下帯幅である. mumasは部分行列の上帯幅であり, mlmas = n - m1 の場合, 無視される. 文献 (1) E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations II. Stiff and differential-algebraic Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1996) (邦訳: 「常微分方程式の数値解法Ⅱ 発展編」スプリンガージャパン (2008))