|
|
◆ dassl()
| void dassl |
( |
int |
n, |
|
|
void(*)(int, double, double *, double *, double *, int *) |
res, |
|
|
double * |
t, |
|
|
double |
y[], |
|
|
double |
yp[], |
|
|
double |
tout, |
|
|
int |
iopt[], |
|
|
double * |
rtol, |
|
|
double * |
atol, |
|
|
void(*)(int, double, double *, double *, int, double *, double) |
jac, |
|
|
double |
work[], |
|
|
int |
lwork, |
|
|
int |
iwork[], |
|
|
int |
liwork, |
|
|
int * |
info |
|
) |
| |
微分代数方程式(DAE) (1〜5次 後退微分公式 (BDF))
注 - 本ルーチンは次バージョンで廃止予定です.
- 目的
- 本ルーチンは陰的な微分代数方程式(DAE)
f(t, y, y') = 0, ただし t = t0 において y = y0 および y' = y0'
本ルーチンでは, 微分y'を1〜5次の後退微分公式(BDF)によって近似し, 各時間ステップにおいて得られた非線形連立方程式をニュートン法より解く.
- 引数
-
| [in] | n | 微分方程式の数. (n >= 1) |
| [in] | res | 微分代数方程式の残差 delta = f(t, y, y') を求めるユーザー定義サブルーチンで, 次のように定義すること. void res(int n, double t, double y[], double yp[], double delta[], int *ires)
{
delta[i] = t, y[] および yp[] における残差 (i = 0 〜 n-1)
}
iresは入力時常に0の整数フラグである. サブルーチンresでは, y[]の値が正しくなかった場合あるいは停止条件に会った場合にのみiresを変更しなけらばならない. 入力値が正しくなかった場合 ires = -1 に設定すると, dasslは ires = -1 を避けて問題を解こうとする. ires = -2 とすると, dasslは呼び出し元プログラムに info = 11 で戻る.
変数 n, t, y[] および yp[] をresの中で変更してはならない. |
| [in,out] | t | 本ルーチンはtからtoutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される.
新たなtoutにおける解を求めるために積分を継続することができる(インターバル・モードとよぶ). その場合, 2回目以降の呼び出しでは, tは前回のtoutに等しくなければならない.
また, toutまでの各中間ステップにおいて戻ることもできる(中間結果出力モードとよぶ). 本モードは解の挙動をみたい場合に使うとよい.
モードはパラメータiopt[2]で指定できる.
[in] 独立変数tの初期値.
[out] 独立変数tの最終ステップの最後の点の値 (インターバル・モードでは通常toutに等しい). 解がこの点まで正常に求められたことを示す. |
| [in,out] | y[] | 配列 y[ly] (ly >= n)
[in] tの初期値における従属変数yの初期値.
[out] 最終のt(インターバル・モードにおいてはtoutに等しい)において求められた解(の近似値). |
| [in,out] | yp[] | 配列 yp[lyp] (lyp >= n)
[in] tの初期値における微分値y'の初期値.
[out] 最終のt(インターバル・モードにおいてはtoutに等しい)において求められた微分値(の近似値). |
| [in] | tout | 解を求めたい点を設定する. tout = t とはできない. 積分を行うのはtについて前進方向(tout > t)でも後退方向(tout < t)でもよい. ただし, 初期化時でなければ積分方向を変更することはできない.
要求精度を満たすように自動的に選ばれたステップ幅を用いてtからtoutに向かって解が求められる. 必要であれば, 各中間ステップにおける解とその微分係数を確認するために戻ることができる(中間結果出力モード). ただし, その場合であっても従来どおりtoutを指定しなければならない. |
| [in,out] | iopt[] | 配列y iopt[liopt] (liopt >= 15)
本配列は問題の解き方についてより詳細に制御するために使われる. 一般的にはすべて 0に設定しておけばよい.
[in]
iopt[0]: プログラムが自動的に設定するパラメータ. ユーザーが設定する必要はない.
iopt[1]: rtolおよびatolがスカラーか配列かを指定する. 現問題について初回呼び出し後にこれを変更することはできない.
= 0: rtol および atolはスカラー.
= 1: rtol および atolは配列.
iopt[2]: 動作モード.
= 0: toutで戻る (インターバル・モード).
= 1: ステップごとに戻る (中間結果出力モード).
iopt[3]: 本ルーチンはtoutを超えて積分を行い内挿によりtoutにおける値を求めることがあるが, ある点において不連続になったり, その点を超えて関数あるいは微分が定義されていなかったりして, その点を超えて積分を行うのが許されないことがある. そのような場合, ストップポイント tstop を設定することができ, その点を超えての積分を行わないようにすることができる.
= 0: ストップポイントを設定しない
= 1: work[0] = tstop をストップポイントとして設定する
iopt[4]: 微分方程式の偏微分を解析的に求めるためのサブルーチン(下記参照)が提供されない場合, ルーチン内で数値計算(差分近似)を行う.
= 0: 偏微分を数値計算(差分近似)により自動的に求める.
= 1: 偏微分を解析的に求めるためのサブルーチンjacをユーザーが用意する.
iopt[5]: 偏微分行列をフル(通常)行列として計算するか帯行列として計算するかを指定する. 一般的には, 2*ml + mu < n であれば帯行列として計算した方が計算時間および記憶容量の点で有利である. ここで, mlは下帯幅, muは上帯幅である.
= 0: n x n フル(通常)行列.
= 1: 帯行列. 帯幅を iwork[0] = ml, iwork[1] = mu と設定すること. (0 <= ml <= n かつ 0 <= mu <= n)
iopt[6]: ユーザーは最大ステップ幅(絶対値)を指定することができる.
= 0: 最大ステップ幅の選択をプログラムが行う.
= 1: 最大ステップ幅の選択を work[1] = hmax と設定することにより行う.
iopt[7]: ユーザーはステップ幅の初期値(絶対値)を指定することができる.
= 0: ステップ幅の初期値の選択をプログラムが行う.
= 1: ステップ幅の初期値の選択を work[2] = h0 と設定することにより行う.
iopt[8]: ユーザーは最大次数を指定することができる.
= 0: 最大次数のデフォルト値(maxord = 5)を使用する.
= 1: 最大次数の選択を iwork[2] = maxord と設定することにより行う (1 <= maxord <= 5).
iopt[9]: 方程式の解が常に非負の場合, このパラメータを設定するとよい.
= 0: 非負の制約を課することなく問題を解く.
= 1: 特別に非負の制約を課して問題を解く.
iopt[10]: t, y および y'の初期値は通常は矛盾のない値である. すなわち, 初期値について f(t, y, y') = 0 が成り立っている. もし微分の初期値が正確にわからない場合, プログラムでそれを解くように試行することができる.
= 0: t, y および y'の初期値は矛盾がない.
= 1: プログラムでy'を計算してみるようにする. yp[]に初期近似を設定しておくこと. もし全く不明であれば, yp[]に0を設定せよ.
[out]
iopt[0]: 本パラメータの値は計算を継続する際に必要である. 変更してはならない. |
| [in,out] | rtol | スカラー(itol = 0 の場合) または 配列 rtol[lrtol] (itol = 1 の場合) (lrtol >= n) (rtolまたはrtol[i] >= 0)
求める解の精度を指定する相対誤差許容値. 本パラメータはiopt[1]に従ってスカラーまたは配列のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる.
許容値は各ステップにおける局所誤差テストに使われ, y[]の各要素がおおむね次式を満たすようにする (i = 0 〜 n-1).
abs(局所誤差) <= rtol*abs(y[i]) + atol (itol = 0 の場合)
または
abs(局所誤差) <= rtol[i]*abs(y[i]) + atol[i] (itol = 1 の場合)
rtol = 0 と設定するとその要素については純粋に絶対誤差テストとなる.
[in] 誤差テストのための相対誤差許容値.
[out] info = 4 で戻ったときには計算を続けるために適当な大きさに調節した値を返す. それ以外の場合には変更されない. |
| [in,out] | atol | スカラー(itol = 0 の場合) または 配列 atol[latol] (itol = 1 の場合) (latol >= n) (atolまたはatol[i] >= 0)
求める解の精度を指定する絶対誤差許容値. 本パラメータはiopt[1]に従ってスカラーまたは配列のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる.
許容値は各ステップにおける局所誤差テストに使われ, y[]の各要素がおおむね次式を満たすようにする (i = 0 〜 n-1).
abs(局所誤差) <= rtol*abs(y[i]) + atol (itol = 0 の場合)
または
abs(局所誤差) <= rtol[i]*abs(y[i]) + atol[i] (itol = 1 の場合)
atol = 0 と設定するとその要素については純粋に相対誤差テストとなる.
[in] 誤差テストのための絶対誤差許容値.
[out] info = 4 で戻ったときには計算を続けるために適当な大きさに調節した値を返す. それ以外の場合には変更されない. |
| [in] | jac | 偏微分行列を解析的に求めるユーザー定義サブルーチンで, 次のように定義すること. void jac(int n, double t, double y[], double yp[], int ldypd, double ypd[][ldypd], double cj)
{
df/dy + cj*df/dy' をypd[][]に設定.
}
ypd[][]の全ての要素はdjacの呼び出し前に0に初期化されるので, 0でない要素だけ設定すればよい.
iopt[4] = 0 が指定された場合, 偏微分行列は数値計算(差分近似)によりルーチン内で計算される.
変数n, t, y[], yp[], ldypdおよびcjをjacの中で変更してはならない. |
| [in,out] | work[] | 配列 work[lwork]
作業領域.
いくつかの要素はプログラムのためのパラメータ値として使用される.
[in]
work[0]: iopt[3] = 1 の場合, ストップポイントを work[0] = tstop と設定する.
work[1]: iopt[6] = 1 の場合, 最大ステップ幅を work[1] = hmax と設定する.
work[2]: iopt[7] = 1 の場合, 初期ステップ幅を work[2] = h0 と設定する.
[out]
work[2]: 次のステップで試されるステップ幅.
work[3]: 独立変数の現在の値. info = 1 であっても内挿が行われた場合にはtと異なる.
work[6]: 直近の成功したステップで使われたステップ幅. |
| [in] | lwork | 配列 work[]のサイズ.
iopt[4] = 0 (jacを使用しない) の場合,
lwork >= n^2 + (maxord + 5)*n + 40 (偏微分がフル行列の場合),
lwork >= (2*ml + mu + 1)*n + (maxord + 5)*n + 40 + 2*(n/(ml + mu + 1) + 1) (偏微分が帯行列の場合);
iopt[4] = 1 (jacを使用する) の場合,
lwork >= 2*n^2 + (maxord + 5)*n + 40 (偏微分がフル行列の場合),
lwork >= 2*(2*ml + mu + 1)*n + (maxord + 5)*n + 40 (偏微分が帯行列の場合);
mlおよびmuについてはiopt[5]を, maxordについてはiopt[8]を参照せよ. |
| [in,out] | iwork[] | 配列 iwork[liwork]
作業領域.
いくつかの要素はプログラムのためのパラメータ値として使用される.
[in]
iwork[0]およびiwork[1]: iopt[5] = 1 の場合, 偏微分の行列の下および上帯幅を iwork[0] = ml および iwork[1] = mu と設定する.
iwork[2]: iopt[8] = 1 の場合, 最大次数を iwork[2] = maxord と設定する.
[out]
iwork[6]: 次のステップで試される次数.
iwork[7]: 直近のステップで使われた次数.
iwork[10]: ステップ数.
iwork[11]: resの呼び出し回数.
iwork[12]: 偏微分行列の評価回数.
iwork[13]: 誤差テストの失敗回数.
iwork[14]: 収束テストの失敗回数 (特異な反復行列による失敗を含む). |
| [in] | liwork | 配列 iwork[]のサイズ. (liwork >= n + 20) |
| [out] | info | [in] 制御コード.
= 0: 問題の最初の呼び出し時(新たに問題を開始する場合)には info = 0 と設定せよ. 初期化が行われ, 新たな問題についての計算が開始される.
= 1, 2, 11〜13: info = 1, 2, 11〜13 で戻った場合, 計算を継続するために, infoの値を変えずに再呼び出しすることができる(下記説明を参照せよ).
[out] リターンコード. ユーザーはこの値をチェックし, info = 1, 2, 11〜13 であれば必要により次のアクションとして本ルーチンを再度呼び出すことができる.
= -1: 入力パラメータ n の誤り (n < 1)
= -6: 入力パラメータ tout の誤り (tout = t)
= -6: 入力パラメータ tout の誤り (toutの値がtに対して積分方向の反対側にある)
= -6: 入力パラメータ tout の誤り (toutがtに近すぎる)
= -7: 入力パラメータ iopt[] の誤り (0または1でないものがある)
= -8: 入力パラメータ rtol または atol の誤り (rtolとatolの全要素が0)
= -8: 入力パラメータ rtol または atol の誤り (重みベクトル (= rtol[i]*abs(y[i]) + atol[i]) の要素に0以下のものがある)
= -8: 入力パラメータ rtol の誤り (rtol < 0 または rtol[i] < 0)
= -9: 入力パラメータ atol の誤り (atol < 0 または atol[i] < 0)
= -11: 入力パラメータ work[0] の誤り (iopt[3] = 1 かつ tstopがtに対して積分方向の反対側にある)
= -11: 入力パラメータ work[0] の誤り (iopt[3] = 1 かつ tstopがtoutに対して積分方向の反対側にある)
= -11: 入力パラメータ work[1] の誤り (hmax < 0)
= -11: 入力パラメータ work[2] の誤り (iopt[7] = 1 かつ h0 = 0)
= -12: 入力パラメータ lwork の誤り (lworkが小さすぎる)
= -13: 入力パラメータ iwork[0] の誤り (ml < 0 または ml > n)
= -13: 入力パラメータ iwork[1] の誤り (mu < 0 または mu > n)
= -13: 入力パラメータ iwork[2] の誤り (maxordが範囲外)
= -14: 入力パラメータ liwork の誤り (liworkが小さすぎる)
= -15: 入力パラメータ info の誤り (info != 0, 1, 2, 11〜13)
= 1: 正常終了 (t = tout). toutを再設定して再呼び出しすることができる. toutはtと異なる値であること.
= 2: 中間結果出力モードのため戻った (まだ toutに達していない). tout方向の次のステップを続行するために再呼び出しすることができる.
= 11: 最大ステップ数(500)を超えた. 再呼び出しして続行することができる. さらに500ステップ続けることができる.
= 12: 誤差許容値が厳しすぎる. 再呼び出しして, 自動的に調節された許容値を使って続行することができる. 許容値を自分で変更してから再呼び出しすることもできる.
= 13: 求められた解が0であるため, 純粋な相対誤差テスト(atol = 0)ができない. atolを正の値にして再呼び出しして続行することができる.
= 14: 誤差テストの失敗を繰り返した.
= 15: 修正子の収束テストの失敗を繰り返した.
= 16: 偏微分行列が特異である.
= 17: 最後の試行ステップにおいて収束テストの失敗を繰り返した.
= 18: ires = -1 のため修正子が収束しない.
= 19: ires = -2 のため終了した.
= 20: 微分の初期値の計算に失敗した.
= 21: 無限ループが検出された. |
- 詳細
- 次の例は, n = 6, ml = 2, mu = 1 の場合の, 本ルーチンで使われる2*ml+mu+1行×n列の帯行列形式を表す.
一般行列形式:
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 * *
+で示された最初の2行は作業領域として使用される. *で示された配列要素は使用されない.
- 出典
- SLATEC (DASSL)
|
1.9.7