_dassl _dassl_r _debdf _debdf_r _radau _radau5 _radau5_r _radau_r _radaup _radaup_r _rodas _rodas_r _seulex _seulex_r

_seulex() void _seulex ( int  n, void(*)(int, double, double *, double *)  f, int  ifcn, double *  t, double  y[], double  tout, double *  rtol, double *  atol, int  itol, void(*)(int, double, double *, int, double *)  jac, int  ijac, int  mljac, int  mujac, void(*)(int, int, double *)  mas, int  imas, int  mlmas, int  mumas, void(*)(int, double, double, double *, int, double *, int *, int *)  solout, int  iout, double  work[], int  lwork, int  iwork[], int  liwork, double  rcont[], int  lrcont, int  icont[], int  licont, int *  info  ) 常微分方程式の初期値問題 (補外法 (線形陰的オイラー法)) 目的 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式)の初期値問題 M * dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. Mは質量マトリックスである. 方程式は線形陰的(M = I(単位行列) でない)でも, 陽的(M = I)でもよい. seulexは線形陰的オイラー法に基づく補外法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている. 詳細については下記参考文献を参照せよ. 引数 [in] n 微分方程式の数. (n >= 1) [in] f 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. void f(int n, double t, const double y[], double yp[]) { yp[i] = tおよびy[]における微分係数の計算値 (i = 0 〜 n-1) } ただし, nは方程式の数, yp[]は与えられたtおよびy[]における微分係数, すなわち, yp[i] = dyi/dt = fi(t, y[0], ..., y[n-1]) (i = 0 〜 n-1). Mは通常単位行列であるが, 必要であればmasで定義することができる. [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] jac ヤコビ行列を計算するユーザー定義サブルーチンで, 次のように定義すること. void jac(int n, double t, double y[], int ldypd, double ypd[][ldypd]) { tおよびy[]におけるヤコビ行列をypd[][]に設定. } ypd[]には, 2次元配列の整合寸法(行数)をldypdとして, mljac = n の場合は通常の2次元配列, mljac < n の場合は帯行列形式でヤコビ行列の2次元データを格納すること. 帯行列形式の詳細は下記を参照のこと. [in] ijac ヤコビ行列の計算方法の選択. = 0: ヤコビ行列を有限差分により内部で求め, jacは呼び出されない. = 1: ヤコビ行列をjacにより計算する. (他の値であれば ijac = 0 とみなす.) [in] mljac ヤコビ行列の下帯幅. (0 <= mljac <= n) mljac = n の場合, ヤコビ行列はn×n行列として格納される. mljac < n の場合, ヤコビ行列は帯行列として格納される. [in] mujac ヤコビ行列の上帯幅. (0 <= mujac <= n) mljac = n の場合, mujacは無視される. [in] mas 質量マトリクスMを与えるユーザー定義サブルーチンで, 次のように定義すること. void mas(int n, int ldam, double am[][ldam]) { 質量マトリクスをam[][]に設定. } am[]には, 2次元配列の整合寸法(行数)をldamとして, mlmas = n の場合は通常の2次元配列, mlmas < n の場合は帯行列形式で質量マトリクスの2次元データを格納すること. 帯行列形式の詳細は下記を参照のこと. [in] imas 質量マトリクスMの選択. = 0: Mは単位行列とする. masは呼び出されない. = 1: Mはmasにより与えられる. (他の値であれば 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] solout 中間結果を出力するユーザー定義サブルーチンで, ステップが成功するごとに呼び出される. 次のように定義すること. void solout(int nr, double told, double t, double y[], int n, double rcont[], int icont[], int *irtrn) { nr番目のステップtでのy[]の値が渡されるのでこれを出力する. toldは前回のtの値, nは方程式の次数である. irtrnの入力値は, 初回, 中間, または, 最終呼び出しのとき, それぞれ 0, 1 または 2である. irtrnを使用して積分を中断することもできる. soloutでirtrnを負の値に設定すると積分を中断して info = 2 で終了する. iout = 1 の場合には, 制御情報 rcont[] および icont[] により密出力を行うことができる. 区間[told, t]内の任意の点t2での解y[i] (0 <= i <= n-1)を次の関数呼び出しにより求めることができる. y[i] = contex(i, t2, rcont, icont); } [in] iout soloutの呼び出し指定. = 0: soloutの呼び出しを行わない = 1: soloutの呼び出しを行う = 2: 密出力のためにsoloutの呼び出しを行う. この場合, iwork[7]に密出力を行う要素数(nrdens)を設定しなければならない. (他の値であれば 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[2]: ヤコビ行列を再計算するかどうかを示す. 計算に時間がかかる場合には, 例えば, 0.01に増やすとよい. 小規模な問題の場合には小さくする. (デフォルト値 = min(0.0001, rtol[0])) work[3] および work[4]: ステップ幅選択パラメータ. j番目の対角エントリーについての新しいステップ幅は次式の範囲で選ばれる. facmin/work[4] <= j番目のhnew/hold <= 1/facmin ただし, facmin = work[3]^(1/(j - 1)) (デフォルト値: work[3] = 0.1, work[4] = 4) work[5], work[6]: 次数選択パラメータ. w[k-2] <= w[k-1]*work[5] の場合, ステップ幅を減らす. w[k-1] <= w[k-2]*work[6] の場合, ステップ幅を増やす. (デフォルト値: work[5] = 0.7, work[6] = 0.9) work[7], work[8]: ステップ幅制御の安全係数. hnew = h*work[8]*(work[7]*tol/err)^(1/(j - 1)) (デフォルト値: work[7] = 0.6, work[8] = 0.93) work[10], work[11], work[12], work[13]: fcn, jac, dec, solそれぞれの呼び出しのワークの推定値. (デフォルト値: work[10] = 1, work[11] = 5, work[12] = 1, work[13] = 1) [out] work[0]: 最後に採用されたステップのステップ幅 [in] lwork 配列 work[]のサイズ. (lwork >= n*(ljac + km + ldypd + 10) + nm1*(lmas + le) + 4*km + km2*nrdens + 20) ただし, km = iwork[2], km2 = km*(km + 1)/2, nrdens = iwork[7] 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]) ldypd = max(ljac, lmas) 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 であれば, ヤコビ行列行列はヘッセンベルグ形に変換される. これは, 特に大規模でヤコビ行列が密行列の場合に有効である. ヤコビ行列が帯行列の場合(mljac < n)と陰的な方程式の場合(imas = 1)には無効である. iwork[1]: 許される最大ステップ数. (デフォルト値 = 100000) (iwork[1] < 0 の場合, デフォルト値を使用する) iwork[2]: 補外表の最大列数. (iwork[2] >= 3) (デフォルト値 = 12) (iwork[2] < 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) (iwork[4] != 0 かつ iwork[4] != 1 の場合, デフォルト値を使用する.) iwork[7]: nrdens = 密出力を行う要素数. iwork[7] < 0 であれば, nrdens = abs(iwork[7]) とする. (デフォルト値 = n) nrdensは iout = 2 の場合のみ有効である. その他の場合, nrdens = 0 とみなす. 0 < nrdens < n の場合, icont[] を設定しなければならない(下記参照). 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 >= 2*n + km + 20 (km: iwork[2]を参照)) liwork < 0 であれば abs(liwork) を使用し, iwork[0]〜iwork[19]を 0 に初期化する. [out] rcont[] 配列 rcont[lcont] 密出力の制御情報のための領域. iout = 2 の場合にのみ使用される. [in] lrcont 配列 cont[]のサイズ. (lrcont >= (km + 2)*nrdens (km: iwork[2]を参照, nrdens: iwork[7]を参照)) [in,out] icont[] 配列 icont[licont] 密出力を行う要素番号のリスト. nrdens = n または iwork[7] < 0 の場合, 自動的に 1, ..., nrdens と設定される. 0 < nrdens < n の場合, 密出力を行う要素番号を設定しておかなければならない. nrdensについてはiwork[7]を参照せよ. [in] licont 配列 icont[]のサイズ. (licont >= nrdens (nrdens: iwork[7]を参照)) [in,out] info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). [out] = -1: 入力パラメータ n の誤り (n < 1) = -7: 入力パラメータ rtol の誤り (rtol < 0 または rtol[i] < 0) = -7: 入力パラメータ rtol または atol の誤り (rtol = 0 かつ atol = 0, または rtol[i] = 0 かつ atol[i] = 0) = -8: 入力パラメータ atol の誤り (atol < 0 または atol[i] < 0) = -16: 入力パラメータ mljac あるいは mlmas の誤り (mlmas > mljac) = -17: 入力パラメータ mujac あるいは mumas の誤り (mumas > mujac) = -21: 入力パラメータ lwork の誤り (lworkが小さすぎる) = -22: 入力パラメータ iwork[8] あるいは iwork[9] の誤り (iwork[8] < 0 または iwork[9] < 0 または iwork[8] + iwork[9] > n) = -23: 入力パラメータ liwork の誤り (liworkが小さすぎる) = -25: 入力パラメータ lrcont の誤り (lrcontが小さすぎる) = -27: 入力パラメータ licont の誤り (licontが小さすぎる) = -28: 入力パラメータ info の誤り (info != 0 かつ info != 1) = 1: 正常終了 = 2: soloutによる中断 (正常終了) = 11: 最大ステップ数を超えた 詳細 次の例は, 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[6]およびiwork[7]の設定)の例. 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)