|
|
◆ Radau_r()
| Sub Radau_r |
( |
N As |
Long, |
|
|
T As |
Double, |
|
|
Y() As |
Double, |
|
|
Tout As |
Double, |
|
|
RTol() As |
Double, |
|
|
ATol() As |
Double, |
|
|
Cont() As |
Double, |
|
|
Info As |
Long, |
|
|
TT As |
Double, |
|
|
YY() As |
Double, |
|
|
YYp() As |
Double, |
|
|
YYpd() As |
Double, |
|
|
Irtrn As |
Long, |
|
|
IRev As |
Long, |
|
|
Optional Nsmin As |
Long, |
|
|
Optional Nsmax As |
Long, |
|
|
Optional Nsus As |
Long, |
|
|
Optional Iout As |
Long, |
|
|
Optional Ijac As |
Long, |
|
|
Optional Mljac As |
Long = -1, |
|
|
Optional Mujac As |
Long, |
|
|
Optional Imas As |
Long, |
|
|
Optional Mlmas As |
Long = -1, |
|
|
Optional Mumas As |
Long, |
|
|
Optional Neval As |
Long, |
|
|
Optional Njac As |
Long, |
|
|
Optional Nstep As |
Long, |
|
|
Optional Naccept As |
Long, |
|
|
Optional Nreject As |
Long, |
|
|
Optional M1 As |
Long, |
|
|
Optional M2 As |
Long, |
|
|
Optional Hinit As |
Double, |
|
|
Optional Hmax As |
Double, |
|
|
Optional MaxIter As |
Long, |
|
|
Optional MaxNit As |
Long, |
|
|
Optional StartN As |
Long, |
|
|
Optional Nind1 As |
Long, |
|
|
Optional Nind2 As |
Long, |
|
|
Optional Nind3 As |
Long, |
|
|
Optional Pred As |
Long, |
|
|
Optional Hess As |
Long, |
|
|
Optional Safe As |
Double, |
|
|
Optional Thet As |
Double, |
|
|
Optional Quot1 As |
Double, |
|
|
Optional Quot2 As |
Double, |
|
|
Optional Facl As |
Double, |
|
|
Optional Facr As |
Double, |
|
|
Optional Vitu As |
Double, |
|
|
Optional Vitd As |
Double, |
|
|
Optional Hhou As |
Double, |
|
|
Optional Hhod As |
Double, |
|
|
Optional Cnt As |
Long |
|
) |
| |
常微分方程式の初期値問題 (次数(5, 9, 13次)自動選択の陰的ルンゲ・クッタ法 (ラダウIIA法)) (リバースコミュニケーション版)
- 目的
- 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式)の初期値問題
M * dy/dt = f(t, y), ただし t = t0 において y = y0
radauはラダウIIA法に基づく次数(5, 9, 13次)自動選択の陰的ルンゲ・クッタ法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている.
詳細については下記参考文献を参照せよ.
Radau_rはRadauのリバースコミュニケーション版である.
- 引数
-
| [in] | N | 微分方程式の数. (N >= 1) |
| [in,out] | T | 本ルーチンはTからToutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される.
[in] 独立変数 t の初期値.
[out] ステップ終了時の最終値(通常Toutに等しい). 解がこの点まで正常に求められたことを示す. 新たなToutで Info = 1 として再呼び出しすることにより, 新たな点まで積分を継続することができる. |
| [in,out] | Y() | 配列 Y(LY - 1) (LY >= N)
[in] Tの初期値における従属変数Y()の初期値.
[out] 最終のT(通常Toutに等しい)において求められた解(の近似値). |
| [in] | Tout | 解を求めたい点を設定する. 積分を行うのはtについて前進方向(Tout > T)でも後退方向(Tout < T)でもよい.
要求精度を満たすように自動的に選ばれたステップ幅を用いてTからToutに向かって解が求められる. |
| [in] | RTol() | 配列 RTol(LRTol - 1) (LRTol >= 1) (RTol()の全要素 >= 0)
求める解の精度を指定する相対誤差許容値. 本パラメータはスカラー(LRTol = 1)またはベクトル(LRTol = N)のどちらでもよい. LRTol = 2, ... または N-1の場合, LRTol = 1 とみなす. LRTol > Nの場合, LRTol = N とみなす. LRTol = N であっても LATol = 1 であれば LRTol = 1 とみなす.
許容値は各ステップにおける局所的な誤差テストに使われ, Y(i)の各要素がおおむね次式を満たすようにする (i = 0 〜 LRTol-1).
abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)
RTol(i) = 0 と設定するとその要素については純粋に絶対誤差テストとなる. RTol(i)とATol(i) (i = 0 〜 LRTol-1) が同時に0になってはいけない. |
| [in] | ATol() | 配列 ATol(LATol - 1) (LATol >= 1) (ATol()の全要素 >= 0)
求める解の精度を指定する絶対誤差許容値. 本パラメータはスカラー(LATol = 1)またはベクトル(LATol = N)のどちらでもよい. LATol = 2, ... または N-1の場合, LATol = 1 とみなす. LATol > Nの場合, LATol = N とみなす. LATol = N であっても LRTol = 1 であれば LATol = 1 とみなす.
許容値は各ステップにおけるローカルな誤差テストに使われ, Y()の各要素がおおむね次式を満たすようにする (i = 0 〜 LATol-1).
abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)
ATol(i) = 0 と設定するとその要素については純粋に相対誤差テストとなる. RTol(i)とATol(i) (i = 0 〜 LRTol-1) が同時に0になってはいけない. |
| [out] | Cont() | 配列 Cont(LCont - 1) (LCont >= (Nsmax + 1)*N)
密出力のための制御情報.
(Iout = 0 のときにも必要) |
| [in,out] | Info | [in]
= 0: 初期化して計算を開始する(新たな問題を解く).
= 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する).
[out]
= -1: パラメータ N の誤り. (N < 1)
= -3: パラメータ Y() の誤り.
= -5: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0)
= -6: パラメータ ATol() の誤り. (ATol(i) < 0)
= -7: パラメータ Cont() の誤り.
= -8: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1)
= -10: パラメータ YY() の誤り.
= -11: パラメータ YYp() の誤り.
= -12: パラメータ YYpd() の誤り.
= -14: パラメータ IRev の誤り. (IRev <> 0, 1, 2, 4, 5)
= -23: パラメータ Mljac あるいは Mlmas の誤り. (Mlmas > Mljac)
= -24: パラメータ Mujac あるいは Mumas の誤り. (Mumas > Mujac)
= -30: パラメータ M1 の誤り. (M1 < 0)
= -31: パラメータ M1 または M2 の誤り. (M2 < 0 または M1 + M2 > N)
= -37: パラメータ Nind1, Nind2 あるいは Nind3 の誤り. (Nind1 + Nind2 + Nind3 <> N)
= 1: 正常終了.
= 2: Irtrn による中断 (正常終了).
= 11: 最大ステップ数を超えた.
= 12: ステップサイズが小さくなり過ぎた.
= 13: 行列が特異になった. |
| [out] | TT | IRev = 1: 再呼び出し時に微分係数値を求めてYYp()に入れるべき点のTの値を返す. |
| [out] | YY() | 配列 YY(LYY - 1) (LYY >= N)
IRev = 1: 再呼び出し時に微分係数値を求めてYYp()に入れるべき点のYの値を返す. |
| [in] | YYp() | 配列 YYp(LYYp - 1) (LYYp >= N)
IRev = 1: T(= TT)およびY(= YY())における微分係数, すなわち, YYp(i) = dyi/dt = fi(TT, YY(0), ..., YY(N-1)) (i = 0 〜 N-1) を計算して, 再呼び出し時に入力する. |
| [in] | YYpd() | 配列 YYpd(LYYpd1 - 1, LYYpd2 - 1) (LYYpd1 >= max(Ljac, Lmas), LYYpd2 >= N)
IRev = 2: TTおよびYY()におけるヤコビ行列(dfi/dyj)を再呼び出し時に入力する.
注 - Mljac = N の場合は通常の2次元配列(Ljac = N), Mljac < N の場合は帯行列形式(Ljac = Mljac + Mujac + 1)で格納する. Ijac = 0 ならば Ljac = 0 とする.
IRev = 4: 質量マトリクスMを再呼び出し時に入力する.
注 - Mlmas = N の場合は通常の2次元配列(Lmas = N), Mlmas < N の場合は帯行列形式(Lmas = Mlmas + Mumas + 1)で格納する. Imas = 0 ならば Lmas = 0 とする. |
| [in,out] | Irtrn | [in] IRev = 5: 処理を中断したい場合以外変更してはならない. Irtrnを負の値に設定すると処理を中断して Info = 2 で終了する.
[out] IRev = 5: 初回, 中間, 最後にIRev = 5 で戻るとき, それぞれ 0, 1, 2 を返す. |
| [in,out] | IRev | リバースコミュニケーションの制御変数
[in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない.
[out] 0 以外の場合, 下記処理を行いIRevを変更せずに再び本ルーチンを呼び出すこと.
= 0: 処理終了. 正常終了かどうかはInfoをチェックすること.
= 1: TTおよびYY()における微分係数の計算値をYYp()に設定する. YYp()以外の変数を変更してはならない.
= 2: TTおよびYY()におけるヤコビ行列を計算しYYpd()に設定する. YYpd()以外の変数を変更してはならない.
= 4: 質量マトリックスMをYYpd()に設定する. YYpd()以外の変数を変更してはならない.
= 5: 途中結果を出力する. 通常はどの変数も変更してはならない(Ioutを参照せよ). |
| [in] | Nsmin,Nsmax | (省略可)
最小および最大段数. (1, 3, 5 または 7) (Nsmin <= Nsmax) (省略時: Nsmin = 3, Nsmax = 7)
(上記以外の値であれば省略時の既定値とみなす) |
| [in,out] | Nsus | (省略可)
[in] 最初のステップの段数. (1, 3, 5 または 7) (Nsmin <= Nsus <= Nsmax) (省略時 = Nsmin)
(上記以外の値であれば省略時の既定値とみなす)
[out] 最後のステップの段数. |
| [in] | Iout | (省略可)
計算の途中結果を出力するかどうか指定する. (省略時 = 0)
= 0: 出力しない. (IRev = 5で戻らない)
= 1: 途中結果出力のために成功したステップごとに IRev = 5 で戻る. リバースコミュニケーション版でないRadaupのSoloutの呼び出しに相当する. 対応する情報は次のとおりである.
Nr = Naccept + 1, Told = 前回のT, T = 現在のT, Y() = 現在のY().
Cont()はSoloutと同様に密出力のために使用する.
Y(i) = Contra_r(i, T2, Cont())
(Ioutが上記以外の値であれば省略時の既定値とみなす) |
| [in] | Ijac | (省略可)
ヤコビ行列の計算方法の選択. (省略時 = 0)
= 0: ヤコビ行列を有限差分により求める. (IRev = 2で戻らない)
= 1: ヤコビ行列を外部でユーザが計算する. (計算が必要なときにIRev = 2で戻る)
(他の値であれば省略時の既定値とみなす) |
| [in] | Mljac | (省略可)
ヤコビ行列の下帯幅. (0 <= Mljac <= N) (省略時 = N)
Mljac = N の場合, ヤコビ行列はN×N行列として格納される. Mljac < N の場合, ヤコビ行列は帯行列として格納される.
(Mljac < 0 または Mljac > N であれば省略時の既定値とみなす) |
| [in] | Mujac | (省略可)
ヤコビ行列の上帯幅. (0 <= Mujac <= N) (省略時 = 0)
Mljac = N の場合, Mujacは無視される.
(Mujac < 0 または Mujac > N であれば省略時の既定値とみなす) |
| [in] | Imas | (省略可)
質量マトリクスMの計算方法の選択. (省略時 = 0)
= 0: Mは単位行列とする. (IRev = 4で戻らない)
= 1: 質量マトリクスMを外部でユーザーが計算する. (計算が必要なときにIRev = 4で戻る)
(他の値であれば省略時の既定値とみなす) |
| [in] | Mlmas | (省略可)
質量マトリクスMの下帯幅. (0 <= Mljac <= N) (省略時 = N)
Mlmas = N の場合, MはN×N行列として格納される. Mlmas < N の場合, Mは帯行列として格納される.
(Mlmas < 0 または Mlmas > N であれば省略時の既定値とみなす) |
| [in] | Mumas | (省略可)
質量マトリクスMの上帯幅. (0 <= Mumas <= N) (省略時 = 0)
Mlmas = N の場合, Mumasは無視される.
(Mumas < 0 または Mumas > N であれば省略時の既定値とみなす) |
| [out] | Neval | (省略可)
関数評価回数. (ヤコビ行列の計算のためのものは含まない) |
| [out] | Njac | (省略可)
ヤコビ行列評価回数. (有限差分の場合を含む) |
| [out] | Nstep | (省略可)
ステップ数. |
| [out] | Naccept | (省略可)
採用されたステップ数. |
| [out] | Nreject | (省略可)
不採用だったステップ数. (最初のステップはカウントされない) |
| [in] | M1,M2 | (省略可)
最初のM1本の方程式が次の形式の場合を考える.
y'(i) = y(i + M2) (i = 0 〜 M1-1)
M1はM2の整数倍とし, かつ残りの方程式がy'(M1), ..., y'(N-1)に陽に依存しない場合, パラメータM1とM2を0以外に設定すると効率よく計算される. (M1 > 0, M2 > 0, M1 + M2 <= N) (省略時は M1 = M2 = 0)
パラメータを0以外に設定としたときには, ヤコビ行列の非ゼロ要素(行 M1〜N-1)だけを(N - M1)×Nの配列に入れなければならない. また, 質量マトリックスMの(N - M1)次元右下ブロックだけを(N - M1)×(N - M1)の配列に入れなければならない. |
| [in] | Hinit | (省略可)
ステップ幅の初期値. (省略時 = 1.0e-6)
H = 1/||f'||である. 初期トランジェントのあるスティフな方程式の場合, 通常は1.0e-3〜1.0e-5 が適切とされる.
(Hinit = 0 であれば1.0e-6とする) |
| [in] | Hmax | (省略可)
ステップ幅の最大値. (省略時 = Tout - T)
(Hmax = 0 であれば省略時の既定値とみなす) |
| [in] | MaxIter | (省略可)
許されるステップ数の最大値. (省略時 = 100000)
(MaxIter <= 0 であれば省略時の既定値とみなす) |
| [in] | MaxNit | (省略可)
陰形式の場合の各ステップにおける解の計算のためのニュートン法の最大反復回数. (省略時 = 7)
(MaxNit <= 0 であれば省略時の既定値とみなす) |
| [in] | StartN | (省略可)
ニュートン法の出発値. (省略時 = 0)
= 0: 外挿した解を使う.
<> 0: 0 を使う.
ニュートン法の収束に問題がある場合には後者が推奨される. |
| [in] | Nind1,Nind2,Nind3 | (省略可)
指数1, 指数2, 指数3変数それぞれの次数. (Nind1 > 0, Nind1 + Nind2 + Nind3 = N) (省略時: Nind1 = N, Nind2 = 0, Nind3 = 0)
指数 > 1 の微分代数方程式(DAE)に関するパラメータである. 関数を定義するサブルーチンは, 指数1, 2および3の変数がこの順番に現れるように作成されていなけらばならない. 常微分方程式(ODE)の場合, Nind1は方程式の次数に一致する. |
| [in] | Pred | (省略可)
ステップ幅戦略の選択. (省略時 = 1)
= 1: モデル予期コントローラ(グスタフソン).
= 2: 伝統的なステップ幅コントロール.
(上記以外の値であれば省略時の既定値とみなす) |
| [in] | Hess | (省略可)
ヤコビ行列をヘッセンベルグ形に変換するかどうか指定する. (省略時 = 0)
= 0: 変換しない.
(本機能はサポートされない. 0 以外の値であれば 0 を指定したものとみなす.) |
| [in] | Safe | (省略可)
ステップ幅推定時の安全係数. (0.001 < Safe < 1) (省略時 = 0.9)
(Safe <= 0.001 あるいは Safe >= 1 であれば省略時の既定値とみなす) |
| [in] | Thet | (省略可)
ヤコビ行列を再計算するかどうか指定する. (Thet < 1) (省略時 = 0.001)
ヤコビ行列の計算に時間がかかる場合には値を大きくすればよい(例えば 0.1). 小規模な場合には小さくすればよい(例えば 0.001). 負の値にすると, ステップが採用されるたびに強制的にヤコビ行列が再計算される.
(Thet = 0 あるいは Thet >= 1 であれば省略時の既定値とみなす) |
| [in] | Quot1,Quot2 | (省略可)
Quot1 < Hnew/Hold < Quot2 であればステップ幅は変更されない. (Quot1 <= 1, Quot2 >= 1) (省略時: Quot1 = 1, Quot2 = 1.2)
小規模な場合 Quot1 = 1, Quot2 = 1.2, 大規模な場合 Quot1 = 0.99, Quot2 = 2 とすればよい.
(Quot1 = 0 または Quot1 > 1, Quot2 = 0 または Quot2 < 1 であればそれぞれ省略時の既定値とみなす) |
| [in] | Facl,Facr | (省略可)
ステップ幅選択パラメータ. (Facl <= 1, Facr >= 1) (省略時: Facl = 0.2, Facr = 8)
Facl < Hnew/Hold < Facr となるようにステップ幅が選ばれる.
(Facl = 0 または Facl > 1, Facr = 0 または Facr < 1 であればそれぞれ省略時の既定値とみなす) |
| [in] | Vitu | (省略可)
縮小因子がこの値より小さいときに次数を増やす. (省略時 = 0.002)
(Vitu = 0 であれば省略時の既定値とみなす) |
| [in] | Vitd | (省略可)
縮小因子がこの値より大きいときに次数を減らす. (省略時 = 0.8)
(Vitd = 0 であれば省略時の既定値とみなす) |
| [in] | Hhou,HHod | (省略可)
ステップ幅の比が Hhod <= hnew/h <= Hhou を満たしたときにのみ次数を減らす. (省略時: Hhou = 1.2, HHod = 0.8)
(Hhou = 0 あるいは Hhod = 0 であれば省略時の既定値とみなす) |
| [in] | Cnt | (省略可)
Neval, Njac, Nstep, Naccept, Nrejectをリセットするタイミングを指定する. (省略時 = 0)
= 0: 本ルーチンが呼び出されるたびにリセットする.
<> 0: Info = 0 で呼び出されたときにだけリセットする. |
- 出典
- 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)
- 使用例 (1)
- 次の常微分方程式の初期値問題(スティフな問題)を解く.
dy1/dt = -2*y1 + y2 - cos(t)
dy2/dt = 1998*y1 - 1999*y2 + 1999*cos(t) - sin(t)
(t = 0 において y1 = 1, y2 = 2)
Sub F2(N As Long, T As Double, Y() As Double, Yp() As Double)
Yp(0) = -2 * Y(0) + Y(1) - Cos(T)
Yp(1) = 1998 * Y(0) - 1999 * Y(1) + 1999 * Cos(T) - Sin(T)
End Sub
Sub Ex_Radau_r()
Const N = 2, Nsmax = 7
Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double
Dim RTol(0) As Double, ATol(0) As Double, Info As Long
Dim Cont((Nsmax + 1) * N) As Double
Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, YYpd(0) As Double
Dim Irtrn As Long, IRev As Long, Iout As Long
RTol(0) = 0.0000000001 '1.0e-10
ATol(0) = RTol(0)
T = 0: Tend = 10: Y(0) = 1: Y(1) = 2
Info = 0
Do
Tout = T + 1
IRev = 0
Do
Call Radau_r(N, T, Y(), Tout, RTol(), ATol(), Cont(), Info, TT, YY(), YYp(), YYpd(), Irtrn, IRev)
If IRev = 1 Then Call F2(N, TT, YY(), YYp())
Loop While IRev <> 0
If Info <> 1 Then
Debug.Print "Error in Radau_r: Info =", Info
Exit Do
End If
Debug.Print T, Y(0), Y(1)
Loop While Tout < Tend
End Sub
- 実行結果
1 0.367879441175306 0.90818174703231
2 0.135335283238768 -0.280811553284687
3 4.97870683693447E-02 -0.940205428204427
4 1.83156388903694E-02 -0.635327981888968
5 6.73794700027563E-03 0.290400132479769
6 2.47875217968704E-03 0.962649034049152
7 9.11881980009936E-04 0.754814110251569
8 3.35462628458022E-04 -0.145164571251503
9 1.23409804293368E-04 -0.911006852111282
10 4.539992985725E-05 -0.839026129195415
- 使用例 (2)
- 次の常微分方程式の初期値問題(スティフな問題)を解く(密出力を使用).
dy1/dt = -2*y1 + y2 - cos(t)
dy2/dt = 1998*y1 - 1999*y2 + 1999*cos(t) - sin(t)
(t = 0 において y1 = 1, y2 = 2)
Sub F2(N As Long, T As Double, Y() As Double, Yp() As Double)
Yp(0) = -2 * Y(0) + Y(1) - Cos(T)
Yp(1) = 1998 * Y(0) - 1999 * Y(1) + 1999 * Cos(T) - Sin(T)
End Sub
Sub Ex_Radau_r_2()
Const N = 2, Nsmax = 7
Dim T As Double, Y(N - 1) As Double, Tend As Double
Dim RTol(0) As Double, ATol(0) As Double, Info As Long
Dim Cont((Nsmax + 1) * N) As Double
Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, YYpd(0) As Double
Dim Irtrn As Long, IRev As Long, Iout As Long, Tout As Double
Dim Y0 As Double, Y1 As Double
RTol(0) = 0.0000000001 '1.0e-10
ATol(0) = RTol(0)
Iout = 1
T = 0: Tend = 10: Y(0) = 1: Y(1) = 2
Tout = 1
Info = 0
IRev = 0
Do
Call Radau_r(N, T, Y(), Tend, RTol(), ATol(), Cont(), Info, TT, YY(), YYp(), YYpd(), Irtrn, IRev, Iout:=Iout)
If IRev = 1 Then
Call F2(N, TT, YY(), YYp())
ElseIf IRev = 5 Then
While T >= Tout
Y0 = Contra_r(0, Tout, Cont())
Y1 = Contra_r(1, Tout, Cont())
Debug.Print Tout, Y0, Y1
Tout = Tout + 1
Wend
End If
Loop While IRev <> 0
If Info <> 1 Then Debug.Print "Error in Radau_r: Info =", Info
End Sub
- 実行結果
1 0.367879435674892 0.908181745212123
2 0.135335280234036 -0.280811557022534
3 0.049787071501545 -0.940205423647481
4 1.83156397566473E-02 -0.63532803510914
5 6.7379470015861E-03 0.290400132815444
6 2.47875208736002E-03 0.962649036642242
7 9.11882020545543E-04 0.754814109823405
8 3.35462904939124E-04 -0.145166180438686
9 1.23411589758396E-04 -0.911009687401114
10 4.53999483279316E-05 -0.839026168472856
|
1.9.6