dassl dassl_r radau radau5 radau5_r radau_r radaup radaup_r rodas rodas_r seulex seulex_r

rodas_r() void rodas_r ( int  n, int  ifcn, double *  t, double  y[], double  tout, double  rtol[], double  atol[], int  itol, int  ijac, int  mljac, int  mujac, int  idfx, int  imas, int  mlmas, int  mumas, int  iout, double  work[], int  lwork, int  iwork[], int  liwork, double  cont[], int  lcont, int *  info, double *  tt, double  yy[], double  yyp[], int  ldyypd, double  yypd[], int *  irtrn, int *  irev  ) 常微分方程式の初期値問題 (4(3)次 ローゼンブロック法) (リバースコミュニケーション版) 注 - 本ルーチンは次バージョンで廃止予定です. 目的 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式)の初期値問題 M * dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. Mは質量マトリックスである. 方程式は線形陰的(M = I(単位行列) でない)でも, 陽的(M = I)でもよい. rodasは4(3)次のローゼンブロック法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている. 詳細については下記参考文献を参照せよ. rodas_rはrodasのリバースコミュニケーション版である. 引数 [in] n 微分方程式の数. (n >= 1) [in] ifcn f(t, y)がtに依存するかどうかを示す. = 0: f(t, y)はtに依存しない (自励系). = 1: f(t, y)はtに依存する (非自励系). (他の値であれば ifcn = 0 とみなす.) [in,out] t 本ルーチンはtからtoutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. [in] 独立変数 t の初期値. [out] ステップ終了時の最終値(通常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] rtol スカラー(itol = 0 の場合) または 配列 rtol[lrtol] (itol = 1 の場合) (lrtol >= n) (rtolまたはrtol[i] >= 0) 相対誤差許容値. y[i]の局所誤差がおおむね次式以下になるようにする.   rtol*abs(y[i]) + atol (i = 0〜n-1) (itol = 0 の場合)     または   rtol[i]*abs(y[i]) + atol[i] (i = 0〜n-1) (itol = 1 の場合). rtolとatol, または, rtol[i]とatol[i] (i = 0〜n-1) が同時に0になってはいけない. [in] atol スカラー(itol = 0 の場合) または 配列 atol[latol] (itol = 1 の場合) (latol >= n) (atolまたはatol[i] >= 0) 絶対誤差許容値. y[i]の局所誤差がおおむね次式以下になるようにする.   rtol*abs(y[i]) + atol (i = 0〜n-1) (itol = 0 の場合)     または   rtol[i]*abs(y[i]) + atol[i] (i = 0〜n-1) (itol = 1 の場合). rtolとatol, または, rtol[i]とatol[i] (i = 0〜n-1) が同時に0になってはいけない. [in] itol rtolおよびatolがスカラーか配列かを指定する. = 0: rtol および atolはスカラー. = 1: rtol および atolは配列. (他の値であれば itol = 0 とみなす.) [in] ijac ヤコビ行列の計算方法の選択. = 0: ヤコビ行列を有限差分により内部で求め, irev = 30 で戻ることはない. = 1: ヤコビ行列計算するために irev = 30 で戻る. (他の値であれば ijac = 0 とみなす.) [in] mljac ヤコビ行列の下帯幅. (0 <= mljac <= n) mljac = n の場合, ヤコビ行列はn×n行列として格納される. mljac < n の場合, ヤコビ行列は帯行列として格納される. [in] mujac ヤコビ行列の上帯幅. (0 <= mujac <= n) mljac = n の場合, mujacは無視される. [in] idfx f(t, y)のtについての偏微分の計算方法の選択. = 0: 偏微分を有限差分により内部で求め, irev = 35 で戻ることはない. = 1: 偏微分を計算するために irev = 35 で戻る. (他の値であれば idfx = 0 とみなす.) [in] imas 質量マトリクスMの選択. = 0: Mは単位行列とする. irev = 40 で戻ることはない. = 1: Mを与えるために irev = 40 で戻る. (他の値であれば imas = 0 とみなす.) [in] mlmas 質量マトリクスMの下帯幅. (0 <= mljac <= n) mlmas = n の場合, Mはn×n行列として格納される. mlmas < n の場合, Mは帯行列として格納される. [in] mumas 質量マトリクスMの上帯幅. (0 <= mumas <= n) mlmas = n の場合, mumasは無視される. [in] iout 中間結果を出力するために irev = 50または51で戻るかどうか指定. = 0: 戻らない = 1: 中間結果出力のために戻る (他の値であれば iout = 0 とみなす.) [in,out] work[] 配列 work[lwork] 作業領域. work[0]〜work[19]はプログラムのためのパラメータ領域として使用される. 入力パラメータを0に設定すると, パラメータ値としてはそれぞれのデフォルト値を使用する. [in] work[0]: ステップ幅の初期値. = 1/(f'のノルム) であり, 初期トランジェントのあるスティフな方程式の場合, 通常1.0e-3〜1.0e-5が適切. この値の選択はシビアではない. (デフォルト値 = 1.0e-4) work[1]: ステップ幅の最大値. (デフォルト値 = tout - t) work[3] および work[4]: ステップ幅選択パラメータ. (work[3] <= 1, work[4] >= 1) work[3] <= hnew/hold <= work[4] となるようにステップ幅が選ばれる. (デフォルト値: work[3] = 0.2, work[4] = 6) work[7]: ステップ幅推定時の安全係数 (0.001 < work[7] < 1) (デフォルト値 = 0.9) [out] work[0]: 最後に採用されたステップのステップ幅 [in] lwork 配列 work[]のサイズ. (lwork >= n*(ljac + 11) + nm1*(lmas + le) + 20) ただし, ljac = nm1 (mljac = nm1 の場合), mljac + mujac + 1 (mljac < nm1 の場合) lmas = 0 (imas = 0 の場合), nm1 (imas = 1 かつ mlmas = nm1 の場合), mlmas + mumas + 1 (imas = 1 かつ mlmas < nm1 の場合) le = nm1 (mljac = nm1 の場合), 2*mljac + mujac + 1 (mljac < nm1 の場合) nm1 = n - m1 (m1 = iwork[8]) lwork < 0 であれば abs(lwork) を使用し, work[0]〜work[19]を 0 に初期化する. [in,out] iwork[] 配列 iwork[liwork] 作業領域. iwork[0]〜iwork[19]はプログラムのためのパラメータ領域として使用される. 入力パラメータを0に設定すると, パラメータ値としてはそれぞれのデフォルト値を使用する. [in] iwork[0]: 係数の選択. (デフォルト値 = 1) = 1: 出典のP452(邦訳 P426)の方法 = 2: 同じ方法で異なる係数を使用 = 3: Gerd Steinebachの係数を使用した方法 (他の値の場合, デフォルト値を使用する) iwork[1]: 許される最大ステップ数. (デフォルト値 = 100000) (iwork[1] < 0 の場合, デフォルト値を使用する) 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 のときにリセットする. [out] iwork[13]: nfcn = 関数評価回数 (ヤコビ行列の計算のためのものは含まない) iwork[14]: njac = ヤコビ行列評価回数 (解析的および数値計算) iwork[15]: nstep = ステップの計算回数 iwork[16]: naccept = 採用されたステップ数 iwork[17]: nreject = 不採用だったステップ数 (誤差評価の結果による) (最初のステップの不採用はカウントされない) iwork[18]: ndec = LU分解の回数 (n次元行列) iwork[19]: nsol = 前進・後退代入の回数 [in] liwork 配列 iwork[]のサイズ. (liwork >= nm1 + 20, ただし nm1 = n - m1 (m1 = iwork[8])) liwork < 0 であれば abs(liwork) を使用し, iwork[0]〜iwork[19]を 0 に初期化する. [out] cont[] 配列 rcont[lcont] 密出力の制御情報のための領域. iout = 1 の場合にのみ使用される. [in] lcont 配列 cont[]のサイズ. (lcont >= 4*n) [in,out] info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). [out] = -1: 入力パラメータ n の誤り (n < 1) = -6: 入力パラメータ rtol の誤り (rtol < 0 または rtol[i] < 0) = -6: 入力パラメータ rtol または atol の誤り (rtol = 0 かつ atol = 0, または rtol[i] = 0 かつ atol[i] = 0) = -7: 入力パラメータ atol の誤り (atol < 0 または atol[i] < 0) = -14: 入力パラメータ mljac あるいは mlmas の誤り (mlmas > mljac) = -15: 入力パラメータ mujac あるいは mumas の誤り (mumas > mujac) = -17: 入力パラメータ work[3] あるいは work[4] の誤り (work[3] > 1 または work[4] < 1) = -17: 入力パラメータ work[7] の誤り (work[7] <= 0.001 または work[7] >= 1) = -18: 入力パラメータ lwork の誤り (lworkが小さすぎる) = -19: 入力パラメータ iwork[0] の誤り (iwork[0] < 0 または iwork[0] > 3) = -19: 入力パラメータ iwork[8] あるいは iwork[9] の誤り (iwork[8] < 0 または iwork[9] < 0 または iwork[8] + iwork[9] > n) = -20: 入力パラメータ liwork の誤り (liworkが小さすぎる) = -22: 入力パラメータ lcont の誤り (lcontが小さすぎる) = -23: 入力パラメータ info の誤り (info != 0 かつ info != 1) = -27: 入力パラメータ ldyypd の誤り (ldyypd < max(ljac, lmas)) = -30: 入力パラメータ irev の誤り (irevの値が無効) = 1: 正常終了 = 2: irtrn による中断 (正常終了) = 11: 最大ステップ数を超えた = 12: ステップサイズが小さくなり過ぎた = 13: 行列が特異になった [out] tt irev = 1〜9: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のtの値を返す. [out] yy[] 配列 yy[lyy] (lyy >= n) irev = 1〜9: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のyの値を返す. [in] yyp[] 配列 yyp[lyyp] (lyyp >= n) irev = 1〜9: t(= tt)およびy(= yy[])における微分係数, すなわち, yyp[i] = dyi/dt = fi(tt, yy[0], ..., yy[n-1]) (i = 0 〜 n-1) を計算して, 再呼び出し時に設定すること. irev = 35: ttおよびyy[]におけるf(t, y)のtについての偏微分を再呼び出し時に設定すること. [in] ldyypd irev = 30 または 40: 2次元配列yypd[][]の整合寸法. (ldyypd >= max(ljac, lmas), ljacおよびlmasについてはlworkを参照せよ) [in] yypd[][] 配列 yypd[lyypd][ldyypd] (lyypd >= n) irev = 30: mljac = n の場合は通常の2次元配列, mljac < n の場合は帯行列形式でヤコビ行列の2次元データを格納すること. 帯行列形式の詳細は下記を参照のこと. irev = 40: mlmas = n の場合は通常の2次元配列, mlmas < n の場合は帯行列形式で質量マトリクスの2次元データを格納すること. 帯行列形式の詳細は下記を参照のこと. [in,out] irtrn [in] irev = 50, 51: 積分を中断したい場合以外変更してはならない. irtrnを負の値に設定すると, 積分を中断して info = 2 で終了する. [out] irev = 50, 51: irev = 50 または 51で最初に戻ったときに0, 最後に戻ったときに2, それ以外のときに1を返す. [in,out] irev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはinfoをチェックすること. = 1〜9: ttおよびyy[]における微分係数の計算値をyyp[]に設定すること. yyp[]以外の変数を変更してはならない. = 30: ttおよびyy[]におけるヤコビ行列を解析的に計算しyypd[][]に設定すること. yypd[][]以外の変数を変更してはならない. = 35: ttおよびyy[]におけるf(t, y)のtについての偏微分を計算しyyp[]に設定すること (idfx = 1 かつ ifcn = 1 の場合). yyp[]以外の変数を変更してはならない. = 40: 質量マトリックスをyypd[][]に設定すること(imass = 1の場合). yypd[][]以外の変数を変更してはならない. = 50, 51: iout = 1の場合, 中間結果を出力するためにステップが成功するごとにこのコードで戻る. tにおけるy[]の値が渡される. ttは前回のtの値である.   制御情報 cont[] を用いて密出力を行うことができる. 区間[tt, t]内の任意の点t2での解y[i] (0 <= i <= n-1)を次の関数呼び出しにより求めることができる.     y[i] = contro(i, t2, cont); 詳細 次の例は, 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[3]およびiwork[4]の設定)の例. 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 の場合, 無視される. 出典 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)