|
|
◆ Seulex()
| Sub Seulex |
( |
N As |
Long, |
|
|
F As |
LongPtr, |
|
|
Ifcn As |
Long, |
|
|
T As |
Double, |
|
|
Y() As |
Double, |
|
|
Tout As |
Double, |
|
|
RTol() As |
Double, |
|
|
ATol() As |
Double, |
|
|
Info As |
Long, |
|
|
Optional Solout As |
LongPtr = NullPtr, |
|
|
Optional Jac As |
LongPtr = NullPtr, |
|
|
Optional Mljac As |
Long = -1, |
|
|
Optional Mujac As |
Long, |
|
|
Optional Mas As |
LongPtr = NullPtr, |
|
|
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 Km As |
Long, |
|
|
Optional Nsequ As |
Long, |
|
|
Optional Lambda As |
Long, |
|
|
Optional Hess As |
Long, |
|
|
Optional Thet As |
Double, |
|
|
Optional Fac1 As |
Double, |
|
|
Optional Fac2 As |
Double, |
|
|
Optional Fac3 As |
Double, |
|
|
Optional Fac4 As |
Double, |
|
|
Optional Safe1 As |
Double, |
|
|
Optional Safe2 As |
Double, |
|
|
Optional Cnt As |
Long |
|
) |
| |
常微分方程式の初期値問題 (補外法 (線形陰的オイラー法))
- 目的
- 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式)の初期値問題
M * dy/dt = f(t, y), ただし t = t0 において y = y0
Seulexは線形陰的オイラー法に基づく補外法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている.
詳細については下記参考文献を参照せよ.
- 引数
-
| [in] | N | 微分方程式の数. (N >= 1) |
| [in] | F | 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub F(N As Long, T As Double, Y() As Double, Yp() As Double)
Yp(i) = TおよびY()における微分係数の計算値 (i = 0 〜 N-1)
End Sub
Mは通常単位行列であるが, 必要であればMasで定義することができる. |
| [in] | Ifcn | f()がtに依存するかどうかを示す.
= 0: tに依存しない dy/dt = f(y) (自励系).
= 1: tに依存する dy/dt = f(y, t) (非自励系). |
| [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になってはいけない. |
| [in,out] | Info | [in]
= 0: 初期化して計算を開始する(新たな問題を解く).
= 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する).
[out]
= -1: パラメータ N の誤り. (N < 1)
= -5: パラメータ Y() の誤り.
= -7: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0)
= -8: パラメータ ATol() の誤り. (ATol(i) < 0)
= -9: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1)
= -15: パラメータ Mljac あるいは Mlmas の誤り. (Mlmas > Mljac)
= -16: パラメータ Mujac あるいは Mumas の誤り. (Mumas > Mujac)
= -22: パラメータ M1 の誤り. (M1 < 0)
= -23: パラメータ M1 または M2 の誤り. (M2 < 0 または M1 + M2 > N)
= 1: 正常終了.
= 2: Soloutによる中断 (正常終了).
= 11: 最大ステップ数を超えた. |
| [in] | Solout | (省略可)
中間結果を出力するユーザー定義サブルーチンで, ステップが成功するごとに呼び出される. 次のように定義すること. (省略時 = NullPtr) Sub Solout(Nr As Long, Told As Double, T As Double, Y() As Double, N As Long, RCont As Double, ICont As Long, Irtrn As Long)
Nr番目のステップTでのY()の値が渡されるのでこれを出力する.
Toldは前回のTの値, Nは方程式の次数である.
Irtrnの値は, 初回, 中間, または, 最終呼び出しのとき, それぞれ 0, 1 または 2 である.
Irtrnを使用して積分を中断することができる. SoloutでIrtrnを負の値に設定して戻ると積分を中断して Info = 2 で終了する.
制御情報 RCont および ICont を使用して密出力を行うことができる.
区間[Told, T]内の任意の点T2での解Y(i) (0 <= i <= N-1)を次の関数呼び出しにより得ることができる.
Y(i) = Contex(i, T2, RCont, ICont);
End Sub
|
| [in] | Jac | (省略可)
ヤコビ行列 dfi(t, y)/dyj を計算するユーザー定義サブルーチンで, 次のように定義すること. (省略時 = NullPtr) Sub Jac(N As Long, T As Double, Y() As Double, Ypd() As Double)
Ypd(i, j) = TおよびY()におけるdfi/dyjの計算値 (i = 0〜N-1, j = 0〜N-1).
End Sub
これを省略した場合(Jac = NullPtrの場合), ヤコビ行列は有限差分により求められる. |
| [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] | Mas | (省略可)
質量マトリクスMを与えるユーザー定義サブルーチンで, 次のように定義すること. (省略時 = NullPtr) Sub Mas(N As Long, Am() As Double)
Am(i, j) = 質量マトリクスMijの値 (i = 0〜N-1, j = 0〜N-1).
End Sub
これを省略した場合(Mas = NullPtrの場合), 質量マトリクスMは単位行列とみなす. |
| [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以外に設定としたときには, Jacにおいてヤコビ行列の非ゼロ要素(行 M1〜N-1)だけを(N - M1)×Nの配列に入れなければならない. また, Masにおいて質量マトリックスMの(N - M1)次元右下ブロックだけを(N - M1)×(N - M1)の配列に入れなければならない. |
| [in] | Hinit | (省略可)
ステップ幅の初期値. (省略時 = 1.0e-6)
H = 1/||f'||である. 初期トランジェントのあるスティフな方程式の場合, 通常は1.0e-2〜1.0e-3 が適切とされる.
(Hinit = 0 であれば1.0e-6とする) |
| [in] | Hmax | (省略可)
ステップ幅の最大値. (省略時 = Tout - T)
(Hmax = 0 であれば省略時の既定値とみなす) |
| [in] | MaxIter | (省略可)
許されるステップ数の最大値. (省略時 = 100000)
(MaxIter <= 0 であれば省略時の既定値とみなす) |
| [in] | Km | (省略可)
補外表の最大列数. (Km >= 3) (省略時 = 12)
(Km < 3 であれば省略時の既定値とみなす) |
| [in] | Nsequ | (省略可)
ステップ数列の選択. (省略時 = 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, ...
(上記以外の値であれば省略時の既定値とみなす) |
| [in] | Lambda | (省略可)
密出力のパラメータλ. (Lambda = 0 または 1) (省略時 = 0)
(Lambda <> 0 かつ Lambda <> 1 であれば省略時の既定値とみなす) |
| [in] | Hess | (省略可)
ヤコビ行列をヘッセンベルグ形に変換するかどうか指定する. (省略時 = 0)
= 0: 変換しない.
<> 0: 変換する.
ヤコビ行列が大規模で密な場合に有効である. 帯行列形式の場合および陰形式の場合(Mas <> NullPtr)には変換を行わない. |
| [in] | Thet | (省略可)
ヤコビ行列を再計算するかどうか指定する. (省略時 = min(0.0001, RTol))
ヤコビ行列の計算に時間がかかる場合には値を大きくすればよい(例えば 0.01). 小規模な場合には小さくすればよい.
(Thet = 0 であれば省略時の既定値とみなす) |
| [in] | Fac1,Fac2 | (省略可)
ステップ幅選択パラメータ. (省略時: Fac1 = 0.1, Fac2 = 4)
j番目の対角エントリーについての新しいステップ幅を次式の範囲で選ぶ.
Facmin/Fac2 <= Hnew(j)/Hold <= 1/Facmin
ただし, Facmin = Fac1^(1/(j - 1)).
(Fac1 = 0 あるいは Fac2 = 0 であればそれぞれ省略時の既定値とみなす) |
| [in] | Fac3,Fac4 | (省略可)
次数選択パラメータ. (省略時: Fac3 = 0.7, Fac4 = 0.9)
W(k-1) <= W(k)*Fac3: ステップ幅を減らす.
W(k) <= W(k-1)*Fac4: ステップ幅を増やす.
(Fac3 = 0 あるいは Fac4 = 0 であればそれぞれ省略時の既定値とみなす) |
| [in] | Safe1,Safe2 | (省略可)
ステップ幅制御の安全係数. (省略時: Safe1 = 0.6, Safe2 = 0.93)
Hnew(j) = H*Safe2*(Safe1*Tol/Err)^(1/(j - 1))
(Safe1 = 0 あるいは Safe2 = 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_Seulex()
Const N = 2
Dim Ifcn As Long, 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, I As Long
RTol(0) = 0.0000000001 '1.0e-10
ATol(0) = RTol(0)
Ifcn = 1
T = 0: Tend = 10: Y(0) = 1: Y(1) = 2
Info = 0
Do
Tout = T + 1
Call Seulex(N, AddressOf F2, Ifcn, T, Y(), Tout, RTol(), ATol(), Info)
If Info <> 1 Then
Debug.Print "Error in Seulex: Info =", Info
Exit Do
End If
Debug.Print T, Y(0), Y(1)
Loop While Tout < Tend
End Sub
- 実行結果
1 0.367879441171429 0.908181747019678
2 0.135335283238996 -0.280811553295101
3 4.97870683828775E-02 -0.940205428219208
4 1.83156388962834E-02 -0.635327981969723
5 6.73794700708669E-03 0.290400132470383
6 2.47875218151303E-03 0.962649038831868
7 9.11881988566409E-04 0.754814136332881
8 3.35462640858368E-04 -0.145164571166839
9 1.23409810850164E-04 -0.911006852073842
10 4.53999366716977E-05 -0.839026129140713
- 使用例 (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_Seulex_2()
Const N = 2
Dim Ifcn As Long, T As Double, Y(N - 1) As Double, Tend As Double
Dim RTol(0) As Double, ATol(0) As Double, Info As Long, I As Long
RTol(0) = 0.0000000001 '1.0e-10
ATol(0) = RTol(0)
Ifcn = 1
T = 0: Tend = 10: Y(0) = 1: Y(1) = 2
Info = 0
Call Seulex(N, AddressOf F2, Ifcn, T, Y(), Tend, RTol(), ATol(), Info, AddressOf SoloutEx)
If Info <> 1 Then Debug.Print "Error in Seulex: Info =", Info
End Sub
Sub SoloutEx(Nr As Long, Told As Double, T As Double, Y() As Double, N As Long, RCont As Double, ICont As Long, Irtrn As Long)
Dim Y0 As Double, Y1 As Double
Static Tout As Double
If Nr = 1 Then Tout = 1
While T >= Tout
Y0 = Contex(0, Tout, RCont, ICont)
Y1 = Contex(1, Tout, RCont, ICont)
Debug.Print Tout, Y0, Y1
Tout = Tout + 1
Wend
End Sub
- 実行結果
1 0.367879441171669 0.90818174654094
2 0.135335283289468 -0.280811548008303
3 4.97870684964011E-02 -0.940205424371012
4 1.83156530207606E-02 -0.635328169713382
5 6.73794698286098E-03 0.290400132401757
6 2.47877855763725E-03 0.962655193543705
7 9.11881827350402E-04 0.754814133775546
8 3.35470947611422E-04 -0.145189025251434
9 1.23409759086634E-04 -0.911006868342305
10 4.53999113068494E-05 -0.839026129166472
|
1.9.6