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

radau5_r() void radau5_r ( int  n, double *  t, double  y[], double  tout, double *  rtol, double *  atol, int  itol, int  ijac, int  mljac, int  mujac, 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  ) 常微分方程式の初期値問題 (5次の陰的ルンゲ・クッタ法 (ラダウIIA法)) (リバースコミュニケーション版) 注 - 本ルーチンは次バージョンで廃止予定です. 目的 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式)の初期値問題 M * dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. Mは質量マトリックスである. 方程式は線形陰的(M = I(単位行列) でない)でも, 陽的(M = I)でもよい. radau5は5次のラダウIIA法に基づく陰的ルンゲ・クッタ法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている. 詳細については下記参考文献を参照せよ. radau5_rはradau5のリバースコミュニケーション版である. 引数 [in] n 微分方程式の数. (n >= 1) [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] 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-6) work[1]: ステップ幅の最大値. (デフォルト値 = tout - t) work[2]: ヤコビ行列を再計算するかどうか指定. (work[2] < 1) ヤコビ行列の計算に時間がかかる場合には値を大きくすればよい(例えば 0.1). 小規模な場合には小さくすればよい(例えば 0.001). 負の値にすると, ステップが採用されるたびに強制的にヤコビ行列が再計算される. (デフォルト値 = 0.001) work[3] および work[4]: ステップ幅選択パラメータ. (work[3] <= 1, work[4] >= 1) work[3] <= hnew/hold <= work[4] となるようにステップ幅が選ばれる. (デフォルト値: work[3] = 0.2, work[4] = 8) work[7]: ステップ幅推定時の安全係数 (0.001 < work[7] < 1) (デフォルト値 = 0.9) work[8] および work[9]: work[8] < hnew/hold < work[9] であればステップ幅は変更されない. (work[8] <= 1, work[9] >= 1) これに加えwork[2]の値を大きくすることにより, 大規模な場合にLU分解と計算時間を節約できる. 小規模な場合, work[8] = 1, work[9] = 1.2 とすればよい. 大規模な場合, work[8]= 0.99, work[9] = 2 とすればよい. (デフォルト値: work[8] = 1, work[9] = 1.2) work[10]: ニュートン法の停止基準で, 通常は work[10] < 1 とする. 小さくすると遅くなるがより安全になる. (デフォルト値 = min(0.03, √rtol[0])) [out] work[0]: 最後に採用されたステップのステップ幅 [in] lwork 配列 work[]のサイズ. (lwork >= n*(ljac + 8) + nm1*(lmas + 3*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]: iwork[0] = 0 とすること. ヤコビ行列をヘッセンベルグ形に変換する機能はサポートされない. iwork[1]: 許される最大ステップ数. (デフォルト値 = 100000) (iwork[1] < 0 の場合, デフォルト値を使用する) iwork[2]: 陰的な方程式の場合の各ステップにおける解の計算のためのニュートン法の最大反復回数. (デフォルト値 = 7) iwork[3]: iwork[3] = 0 であれば外挿した解がニュートン法の出発値として使われる. iwork[3] = 1 であれば 0 が出発値として使われる. ニュートン法の収束に問題がある場合, 後者が推奨される. (デフォルト値 = 0) (他の値の場合, デフォルト値を使用する) iwork[4], iwork[5]およびiwork[6]: 指数1, 2 および 3変数の次数. (iwork[4] > 0, iwork[4] + iwork[5] + iwork[6] = n) 指数 > 1 の微分代数方程式(DAE)に関するパラメータである. 関数を定義するサブルーチンは, 指数1, 2および3の変数がこの順番に現れるように作成されていなけらばならない. 誤差推定の際, 指数2変数にはhを, 指数3変数にはh^2を乗算する. 常微分方程式(ODE)の場合, iwork[4]は方程式の数に一致する. (デフォルト値: iwork[4] = n, iwork[5] = 0, iwork[6] = 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分解の回数 iwork[19]: nsol = 前進・後退代入の回数 (ステップ幅選択のためのものはカウントされない) [in] liwork 配列 iwork[]のサイズ. (liwork >= 3*nm1 + 20, ただし nm1 = n - m1 (m1 = iwork[8])) liwork < 0 であれば abs(liwork) を使用し, iwork[0]〜iwork[19]を 0 に初期化する. [out] cont[] 配列 rcont[lcont] 密出力の制御情報のための領域. iout = 0 であっても必要である. [in] lcont 配列 cont[]のサイズ. (lcont >= 4*n) [in,out] info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). [out] = -1: 入力パラメータ n の誤り (n < 1) = -5: 入力パラメータ rtol の誤り (rtol < 0 または rtol[i] < 0) = -5: 入力パラメータ rtol または atol の誤り (rtol = 0 かつ atol = 0, または rtol[i] = 0 かつ atol[i] = 0) = -6: 入力パラメータ atol の誤り (atol < 0 または atol[i] < 0) = -12: 入力パラメータ mljac あるいは mlmas の誤り (mlmas > mljac) = -13: 入力パラメータ mujac あるいは mumas の誤り (mumas > mujac) = -15: 入力パラメータ work[2] の誤り (work[2] >= 1) = -15: 入力パラメータ work[3] あるいは work[4] の誤り (work[3] > 1 または work[4] < 1) = -15: 入力パラメータ work[7] の誤り (work[7] <= 0.001 または work[7] >= 1) = -15: 入力パラメータ work[8] あるいは work[9] の誤り (work[8] > 1 または work[9] < 1) = -15: 入力パラメータ work[10] の誤り (work[10]が小さすぎる) = -16: 入力パラメータ lwork の誤り (lworkが小さすぎる) = -17: 入力パラメータ iwork[4], iwork[5] あるいは iwork[6] の誤り (iwork[4] + iwork[5] + iwork[6] != n) = -17: 入力パラメータ iwork[8] あるいは iwork[9] の誤り (iwork[8] < 0 または iwork[9] < 0 または iwork[8] + iwork[9] > n) = -18: 入力パラメータ liwork の誤り (liworkが小さすぎる) = -20: 入力パラメータ lcont の誤り (lcontが小さすぎる) = -21: 入力パラメータ info の誤り (info != 0 かつ info != 1) = -25: 入力パラメータ ldyypd の誤り (ldyypdが小さすぎる) = -28: 入力パラメータ irev の誤り (irevの値が無効) = 1: 正常終了 = 2: irtrn による中断 (正常終了) = 11: 最大ステップ数を超えた = 12: ステップサイズが小さくなり過ぎた = 13: 行列が特異になった [out] tt irev = 1〜8: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のtの値を返す. [out] yy[] 配列 yy[lyy] (lyy >= n) irev = 1〜8: 再呼び出し時に微分係数値を求めてyyp[]に入れるべき点のyの値を返す. [in] yyp[] 配列 yyp[lyyp] (lyyp >= n) irev = 1〜8: t(= tt)およびy(= yy[])における微分係数, すなわち, yyp[i] = dyi/dt = fi(tt, yy[0], ..., yy[n-1]) (i = 0 〜 n-1) を計算して, 再呼び出し時に設定すること. [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〜8: ttおよびyy[]における微分係数の計算値をyyp[]に設定すること. yyp[]以外の変数を変更してはならない. = 30: ttおよびyy[]におけるヤコビ行列を解析的に計算しyypd[][]に設定すること. yypd[][]以外の変数を変更してはならない. = 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] = contr5(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[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 の場合, 無視される. 出典 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)