debdf debdf_r radaua radaua_r rodasa rodasa_r seulexa seulexa_r

radaua() void radaua ( int  n, void(*)(int, double, double *, double *)  f, double *  t, double  y[], double  tout, double  tend, double *  rtol, double *  atol, int  itol, void(*)(int, double, double *, int, double *)  fjac, int  ijac, int  mljac, int  mujac, void(*)(int, int, double *)  fmas, int  imas, int  mlmas, int  mumas, int  mode, double  work[], int  lwork, int  iwork[], int  liwork, int *  info  ) 常微分方程式の初期値問題 (5, 9, 13次 可変次数陰的ルンゲ・クッタ法 (ラダウIIA法)) 目的 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式) の初期値問題 M * dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, y は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. M は質量マトリックスである. また, t0 および y0 はそれぞれ t および y の既知の初期値である. 本ルーチンは, 5, 9, 13次 可変次数陰的ルンゲ・クッタ法 (ラダウIIA法) プログラム RADAU (文献 (1)) を書き直したものである. 主として密出力機能の使用法を変更した. 引数 [in] n 微分方程式の数. (n >= 1) [in] f 微分方程式の関数値を求めるユーザー定義サブルーチンで, 次のように定義すること. void f(int n, double t, const double y[], double yp[]) { dy/dt を求め yp[] に入れる. } ただし, n は方程式の数, dy/dt (= f(t, y[])) は与えられた t および y[] における微分の計算値である. [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] fjac ヤコビ行列を求めるユーザー定義サブルーチンで, 次のように定義すること. void fjac(int n, double t, double y[], int ldypd, double ypd[][ldypd]) { t および y[] におけるヤコビ行列を求め ypd[][] に入れる. } ypd[] には, mljac = n の場合にはフル行列形式(通常の n×n 2次元配列), mljac < n の場合には帯行列形式(下記詳細参照)でヤコビ行列を格納すること. [in] ijac ヤコビ行列の計算方法を指定する. = 0: ヤコビ行列を有限差分近似により求める. fjac は参照されない. = 1: ヤコビ行列を fjac により求める. [in] mljac ヤコビ行列の下帯幅. (0 <= mljac <= n) mljac = n の場合, ヤコビ行列はフル行列形式で格納される. mljac < n の場合, ヤコビ行列は帯行列形式で格納される. [in] mujac ヤコビ行列の上帯幅. (0 <= mujac <= n) mljac = n の場合, mujac は無視される. [in] fmas 質量マトリクス M を与えるユーザー定義サブルーチンで, 次のように定義すること. void fmas(int n, int ldam, double am[][ldam]) { 質量マトリクスを am[][] に入れる. } am[] には, mlmas = n の場合にはフル行列形式, mlmas < n の場合には帯行列形式で質量マトリクスを格納すること. [in] imas 質量マトリクス M が単位行列であるかどうかを指定する. = 0: M は単位行列である. fmas は参照されない. = 1: M は fmas により与えられる. [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-6) work[2]: ステップ幅の最大値 (デフォルト値 = abs(tend - t)) work[3]: ヤコビ行列を再計算するかどうかを指定する. (work[3] < 1) (デフォルト値 = 0.001) ヤコビ行列の計算に時間がかかる場合には値を大きくすればよい (例えば 0.1). 小規模な場合には小さくすればよい (例えば 0.001). 負の値にするとステップが採用されるたびに強制的にヤコビ行列が再計算される. work[4], work[5]: ステップ幅選択パラメータ. (work[4] <= 1, work[5] >= 1) (デフォルト値: work[4] = 0.2, work[5] = 8) work[4] <= hnew/hold <= work[5] となるようにステップ幅が選ばれる. work[8]: ステップ幅推定時の安全係数 (0.001 < work[8] < 1) (デフォルト値 = 0.9) work[11], work[12]: work[11] < hnew/hold < work[12] であればステップ幅は変更されない. (work[11] <= 1, work[12] >= 1) (デフォルト値: work[11] = 1, work[12] = 1.2) 大規模な場合には, work[11]= 0.99, work[12] = 2 が推奨される. これに加え work[3] の値を大きくすることによりLU分解と計算時間を節約できる. 小規模な場合にはデフォルト値でよい. work[13]: 縮小因子がこの値より小さいときに次数を増やす. (デフォルト値 = 0.002) work[14]: 縮小因子がこの値より大きいときに次数を減らす. (デフォルト値 = 0.8) work[15], work[16]: ステップ幅の比が work[15] <= hnew/h <= work[16] を満たしたときにのみ次数を減らす. (デフォルト値: work[15] = 1.2, work[16] = 0.8) [out] work[1]: 最後に使用されたステップ幅. [in] lwork 配列 work[]のサイズ. (abs(lwork) >= n*(3 + 3*nsmax + ldjac) + nm1*(nsmax*lde + ldmas) + 3*nsmax + 60 + n*(max(ldjac, ldmas) + 2)) ただし, 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]) nsmax = iwork[10]/10 lwork < 0 であればその絶対値を採用し, work[0]〜work[19] をすべて 0 に設定する. [in,out] iwork[] 配列 iwork[liwork] 整数作業領域. ただし, iwork[0]〜iwork[19] はプログラムのためのパラメータ領域として使用される. info = 0 の呼び出し時に 0 に設定されている入力パラメータはそのデフォルト値が使用される. [in] iwork[0]: ヤコビ行列をヘッセンベルグ形に変換するかどうか指定する. 大規模でフル行列の場合にこの変換は有効である. 帯行列の場合および IMAS = 1 の場合には無効である. = 0: 変換を行う. = 1: 変換を行わない. iwork[1]: 許される最大ステップ数. (デフォルト値 = 100000) iwork[2]: 陰的な方程式を解くために使用されるニュートン法の最大反復回数. (デフォルト値 = 7) iwork[3]: ニュートン法の出発値を指定する. (デフォルト値 = 0) = 0: 外挿した解を出発値として使う. = 1: 0 を出発値として使う. ニュートン法の収束に問題がある場合, 後者が推奨される. iwork[4], iwork[5], iwork[6]: 指数 1, 2 および 3 変数の次数. (iwork[4] > 0, iwork[4] + iwork[5] + iwork[6] = n) (デフォルト値: iwork[4] = n, iwork[5] = 0, iwork[6] = 0) 指数 > 1 の微分代数方程式 (DAE) に関するパラメータ. 関数を定義するサブルーチンは, 指数 1, 2 および 3 の変数がこの順番に現れるように作成されていなけらばならない. 誤差推定の際, 指数 2 変数には h を, 指数 3 変数には h^2 を乗算する. 常微分方程式 (ODE) の場合には iwork[4] = n とする. 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[10]: 最小段数 nsmin および 最大段数 nsmax (= 3, 5 または 7). 10*nsmax + nsmin を設定する. (デフォルト値 = 73 (nsmin = 3, nsmax = 7)) 段数を ns とすると, その値により 2*ns - 1 次の陰的ルンゲ・クッタ法に対応する. iwork[11]: 最初のステップの ns の値. (3, 5 または 7) (デフォルト値 = nsmin) 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) >= (2 + (nsmax - 1)/2)*nm1 + 80, ただし, nsmax = iwork[10]/10, 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: (エラー) 行列が繰り返し特異になった. 詳細 次の例は, 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))