debdf debdf_r radaua radaua_r rodasa rodasa_r seulexa seulexa_r

seulexa_r() void seulexa_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  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  ) 常微分方程式の初期値問題 (補外法 (線形陰的オイラー法)) (リバースコミュニケーション版) 目的 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式) の初期値問題 M * dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, y は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. M は質量マトリックスである. また, t0 および y0 はそれぞれ t および y の既知の初期値である. 本ルーチンは, 線形陰的オイラー法に基づく補外法 プログラム SEULEX (文献 (1)) を書き直したものである. 主として密出力機能の使用法を変更した. 本ルーチンは seulexa のリバースコミュニケーション版である. 引数 [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] 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-6) work[2]: ステップ幅の最大値 (デフォルト値 = abs(tend - t)) work[3]: ヤコビ行列を再計算するかどうかを指定する. (デフォルト値 = min(0.0001, rtol[0])) ヤコビ行列の計算に時間がかかる場合には値を大きくすればよい (例えば 0.01). 小規模な場合には小さくする. work[4], work[5]: ステップ幅選択パラメータ. (デフォルト値: work[4] = 0.1, work[5] = 4) facmin/work[5] <= j番目のhnew/hold <= 1/facmin となるようにステップ幅が選ばれる (ただし, facmin = work[4]^(1/(j - 1))). work[6]: 次数選択パラメータ (次数を減らす判定の係数). (デフォルト値 = 0.7) work[7]: 次数選択パラメータ (次数を増やす判定の係数). (デフォルト値 = 0.9) work[8], work[9]: ステップ幅制御の安全係数. (デフォルト値: work[8] = 0.6, work[9] = 0.93) hnew = h*work[9]*(work[8]*tol/err)^(1/(j - 1)) とする. work[11], work[12], work[13], work[14]: f, fjac, LU分解, 解の計算 それぞれの呼び出しの仕事量の見積り. (デフォルト値: work[11] = 1, work[12] = 5, work[13] = 1, work[14] = 1) [out] work[1]: 最後に使用されたステップ幅 [in] lwork 配列 work[]のサイズ. (abs(lwork) >= n*(ldjac + 2*km + 10) + nm1*(ldmas + lde) + 4*km + n*(km*(km + 1))/2 + 50) ただし, 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]). km = iwork[2]. 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]: 補外表の最大列数 km. (km >= 3) (デフォルト値 = 12) (km < 3 の場合, デフォルト値が使用される.) iwork[3]: ステップ数列の選択. (デフォルト値 = 2) = 1: 1, 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, ... = 2: 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, ... = 3: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, ... = 4: 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ... (他の値であれば, デフォルト値が使用される) iwork[4]: 密出力のパラメータλ. (0 または 1) (デフォルト値 = 0) (0 または 1 以外であればデフォルト値が使用される.) 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) >= 2*n + km + 70, ただし, km = iwork[2]) 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: (エラー) 最大ステップ数を超えた. [out] tt irev = 1〜20 の場合: yyp[]を参照のこと. irev = 30 の場合: yypd[][]を参照のこと. [out] yy[] 配列 yy[lyy] (lyy >= n) irev = 1〜20 の場合: yyp[]を参照のこと. irev = 30 の場合: yypd[][]を参照のこと. [in] yyp[] 配列 yyp[lyyp] (lyyp >= n) irev = 1〜20 の場合: t = tt および y = yy[] における微分値 dy/dt = f(t, y) を計算して, 再呼び出し時に 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〜20: t = tt および y = yy[] における微分の計算値を yyp[] に設定する. yyp[] 以外の変数を変更してはならない. = 30: t = tt および y = yy[] におけるヤコビ行列の計算値を yypd[][] に設定する. yypd[][] 以外の変数を変更してはならない. = 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))