|
|
◆ Odex2_r()
| Sub Odex2_r |
( |
N As |
Long, |
|
|
T As |
Double, |
|
|
Y() As |
Double, |
|
|
Yp() As |
Double, |
|
|
Tout As |
Double, |
|
|
RTol() As |
Double, |
|
|
ATol() As |
Double, |
|
|
RCont() As |
Double, |
|
|
ICont() As |
Long, |
|
|
Info As |
Long, |
|
|
TT As |
Double, |
|
|
YY() As |
Double, |
|
|
YYp2() As |
Double, |
|
|
Irtrn As |
Long, |
|
|
IRev As |
Long, |
|
|
Optional Iout As |
Long = 0, |
|
|
Optional Neval As |
Long, |
|
|
Optional Nstep As |
Long, |
|
|
Optional Naccept As |
Long, |
|
|
Optional Nreject As |
Long, |
|
|
Optional Hinit As |
Double = 0, |
|
|
Optional Hmax As |
Double = 0, |
|
|
Optional MaxIter As |
Long = 0, |
|
|
Optional Km As |
Long = 0, |
|
|
Optional Nsequ As |
Long = 0, |
|
|
Optional Iderr As |
Long = 0, |
|
|
Optional Mudif As |
Long = 0, |
|
|
Optional Fac1 As |
Double = 0, |
|
|
Optional Fac2 As |
Double = 0, |
|
|
Optional Fac3 As |
Double = 0, |
|
|
Optional Fac4 As |
Double = 0, |
|
|
Optional Safe1 As |
Double = 0, |
|
|
Optional Safe2 As |
Double = 0, |
|
|
Optional Safe3 As |
Double = 0, |
|
|
Optional Cnt As |
Long = 0 |
|
) |
| |
常微分方程式の初期値問題 (補外法) (2階微分方程式用) (リバースコミュニケーション版)
- 目的
- 本ルーチンは2階の常微分方程式の初期値問題
d2y/dt2 = f(t, y), ただし t = t0 において y = y0, y' = y'0
Odex2はシュテルマー公式にもとづく補外法のプログラムである. ステップ制御, 次数選択および密出力機能を持っている.
詳細については下記参考文献を参照せよ.
Odex2_rはOdex2のリバースコミュニケーション版である.
- 引数
-
| [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,out] | Yp() | 配列 Yp(LYp - 1) (LYp >= N)
[in] Tの初期値における微分値y'1〜y'nの初期値.
[out] 最終のT(通常Toutに等しい)において求められた微分値(の近似値). |
| [in] | Tout | 解を求めたい点を設定する. 積分を行うのはtについて前進方向(Tout > T)でも後退方向(Tout < T)でもよい.
要求精度を満たすように自動的に選ばれたステップ幅を用いてTからToutに向かって解が求められる. |
| [in] | RTol() | 配列 RTol(LRTol - 1) (LRTol >= 1) (RTol()の全要素 >= 0)
求める解の精度を指定する相対誤差許容値. 本パラメータはスカラー(LRTol = 1)またはベクトル(LRTol = 2N)のどちらでもよい. LRTol = 2, ..., 2N-1の場合, LRTol = 1 とみなす. LRTol > 2N の場合, LRTol = 2N とみなす. LRTol = 2N であっても LATol = 1 であれば LRTol = 1 とみなす.
許容値は各ステップにおける局所的な誤差テストに使われ, Y()およびYp()の各要素がおおむね次式を満たすようにする.
スカラーの場合：
abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)
abs(Yp(i)の局所誤差) <= RTol(i+N)*abs(Y(i)) + ATol(i+N)
ベクトルの場合 (i = 0 〜 N-1)：
abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)
abs(Yp(i)の局所誤差) <= RTol(i+N)*abs(Y(i)) + ATol(i+N)
RTol(i) = 0 と設定するとその要素については純粋に絶対誤差テストとなる. RTol(i)とATol(i) (i = 0 〜 LRTol-1) が同時に0になってはいけない. |
| [in] | ATol() | 配列 ATol(LATol - 1) (LATol >= 1) (ATol()の全要素 >= 0)
求める解の精度を指定する絶対誤差許容値. 本パラメータはスカラー(LATol = 1)またはベクトル(LATol = 2N)のどちらでもよい. LATol = 2, ..., 2N-1の場合, LATol = 1 とみなす. LATol > 2N の場合, LATol = 2N とみなす. LATol = 2N であっても LRTol = 1 であれば LATol = 1 とみなす.
許容値は各ステップにおける局所的な誤差テストに使われ, Y()およびYp()の各要素がおおむね次式を満たすようにする.
スカラーの場合：
abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)
abs(Yp(i)の局所誤差) <= RTol(i+N)*abs(Y(i)) + ATol(i+N)
ベクトルの場合 (i = 0 〜 N-1)：
abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)
abs(Yp(i)の局所誤差) <= RTol(i+N)*abs(Y(i)) + ATol(i+N)
ATol(i) = 0 と設定するとその要素については純粋に相対誤差テストとなる. RTol(i)とATol(i) (i = 0 〜 LRTol-1) が同時に0になってはいけない. |
| [in,out] | RCont() | 配列 RCont(LRCont - 1) (LRCont >= (2*Km + 6)*N)
密出力のための制御情報.
(Iout = 0 のときには参照されない) |
| [in,out] | ICont() | 配列 ICont(LICont - 1) (LICont >= N)
密出力のための制御情報.
(Iout = 0 のときには参照されない) |
| [in,out] | Info | [in]
= 0: 初期化して計算を開始する(新たな問題を解く).
= 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する).
[out]
= -1: パラメータ N の誤り. (N < 1)
= -3: パラメータ Y() の誤り.
= -4: パラメータ Yp() の誤り.
= -6: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0)
= -7: パラメータ ATol() の誤り. (ATol(i) < 0)
= -8: パラメータ RCont() の誤り.
= -9: パラメータ ICont() の誤り.
= -10: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1)
= -12: パラメータ YY() の誤り.
= -13: パラメータ YYp2() の誤り.
= -15: パラメータ IRev の誤り. (IRev != 0, 1, 5)
= 1: 正常終了.
= 2: Irtrn による中断 (正常終了).
= 11: 最大ステップ数を超えた. |
| [out] | TT | IRev = 1: 再呼び出し時に2次微分係数値を求めてYYp2()に入れるべき点のTの値を返す. |
| [out] | YY() | 配列 YY(LYY - 1) (LYY >= N) IRev = 1: 再呼び出し時に2次微分係数値を求めてYYp2()に入れるべき点のYの値を返す. |
| [in] | YYp2() | 配列 YYp2(LYYp2 - 1) (LYYp2 >= N)
IRev = 1: T(= TT)およびY(= YY())における2次微分係数, すなわち, YYp2(i) = d2yi/dt2 = fi(TT, YY(0), ..., YY(N-1)) (i = 0 〜 N-1) を計算して, 再呼び出し時に入力する. |
| [in,out] | Irtrn | [in] IRev = 5: 以下の特殊ケース以外は変更してはならない.
Irtrnを負の値に設定すると積分を中断して Info = 2 で終了する. また, 解を変更した場合には Irtrn = 3 と設定する必要がある.
[out] IRev = 5: 初回, 中間, 最後にIRev = 5 で戻るとき, それぞれ 0, 1, 2 を返す. |
| [in,out] | IRev | リバースコミュニケーションの制御変数
[in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない.
[out] 0 以外の場合, 下記処理を行いIRevを変更せずに再び本ルーチンを呼び出すこと.
= 0: 処理終了. 正常終了かどうかはInfoをチェックすること.
= 1: TTおよびYY()における2次微分係数の計算値をYYp2()に設定する. YYp2()以外の変数を変更してはならない.
= 5: 途中結果を出力する. 通常はどの変数も変更してはならない(Ioutを参照せよ). |
| [in] | Iout | (省略可)
計算の途中結果を出力するかどうか指定する. (省略時 = 0)
= 0: 出力しない. (IRev = 5では戻らない)
= 1: 途中結果出力のために成功したステップごとに IRev = 5 で戻る. リバースコミュニケーション版でないOdex2のSoloutの呼び出しに相当する. 対応する情報は次のとおりである.
Nr = Naccept + 1, Told = 前回のT, T = 現在のT, Y() = 現在のY().
RCont()およびICont()はSoloutと同様に密出力のために使用する.
Y(i) = Contx1_r(i, T2, RCont(), ICont())
(Ioutが上記以外の値であれば省略時の既定値とみなす) |
| [out] | Neval | (省略可)
関数評価回数. |
| [out] | Nstep | (省略可)
ステップ数. |
| [out] | Naccept | (省略可)
採用されたステップ数. |
| [out] | Nreject | (省略可)
不採用だったステップ数. (最初のステップはカウントされない) |
| [in] | Hinit | (省略可)
ステップ幅の初期値. (省略時 = 0.0001)
H = 1/||f'|| であり, 通常 0.1〜0.001 が適切とされる.
(|Hinit| < 0.0001 であれば Hinit = 0.0001 とみなす) |
| [in] | Hmax | (省略可)
ステップ幅の最大値. (省略時 = Tout - T)
(Hmax = 0 であれば省略時の既定値とみなす) |
| [in] | MaxIter | (省略可)
許されるステップ数の最大値. (省略時 = 100000)
(MaxIter <= 0 であれば省略時の既定値とみなす) |
| [in] | Km | (省略可)
補外表の最大列数. (Km >= 3) (省略時 = 9)
(Km < 3 であれば省略時の既定値とみなす) |
| [in] | Nsequ | (省略可)
ステップ数列の選択. (省略時 = 1 (Solout2不使用時), 4 (Solout2使用時))
= 1: 2, 4, 6, 8, 10, 12, 14, 16, ...
= 2: 2, 4, 8, 12, 16, 20, 24, 28, ...
= 3: 2, 4, 6, 8, 12, 16, 24, 32, ...
= 4: 2, 6, 10, 14, 18, 22, 26, 30, ...
Solout2使用時には1〜3は指定できない.
(上記以外の値, あるいは, Solout2使用時に1〜3であれば省略時の既定値とみなす) |
| [in] | Iderr | (省略可)
密出力の式における誤差推定. (省略時 = 0)
= 0: 行う
= 1: 行わない
Solout2を使用しない場合, Iderrの指定にかかわらず誤差推定は行わない.
(上記以外の値であれば Iderr = 1 とみなす) |
| [in] | Mudif | (省略可)
補外式の次数を決めるパラメータ. (1 <= Mudif <= 8) (省略時 = 6)
(Mudif < 1 または Mudif > 8 であれば省略時の既定値とみなす) |
| [in] | Fac1 | (省略可) |
| [in] | Fac2 | (省略可)
ステップ幅選択パラメータ. (省略時: Fac1 = 0.02, Fac2 = 4)
j番目の要素が Facmin/Fac2 <= j番目のHnew/Hold <= 1/Facmin となるように新しいステップ幅が選ばれる. ただし, Facmin = Fac1^(1/(2*j - 1)) である.
(Fac1 = 0 あるいは Fac2 = 0 であればそれぞれ省略時の既定値とみなす) |
| [in] | Fac3 | (省略可) |
| [in] | Fac4 | (省略可)
次数選択パラメータ. (省略時: Fac3 = 0.8, Fac4 = 0.9)
W(k-1) <= W(k)*Fac3 : 次数を減らす
W(k) <= W(k-1)*Fac4 : 次数を増やす
W(k) = A(k)/H(k)は単位ステップあたりの仕事量である. (A(k)は関数値計算回数, H(k)はステップ幅を表す)
(Fac3 = 0 あるいは Fac4 = 0 であればそれぞれ省略時の既定値とみなす) |
| [in] | Safe1 | (省略可) |
| [in] | Safe2 | (省略可)
ステップ幅制御の安全係数. (省略時: Safe1 = 0.65, Safe2 = 0.94)
(Safe1 = 0 あるいは Safe2 = 0 であればそれぞれ省略時の既定値とみなす) |
| [in] | Safe3 | (省略可)
ステップ幅を減らす際の係数. (省略時 = 0.5)
(Safe3 <= 0 であれば省略時の既定値とみなす. Safe3 > 1 であれば 1 とみなす) |
| [in] | Cnt | (省略可)
Neval, Nstep, Naccept, Nrejectをリセットするタイミングを指定する. (省略時 = 0)
= 0: 本ルーチンが呼び出されるたびにリセットする.
<> 0: Info = 0 で呼び出されたときにだけリセットする. |
- 出典
- E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations. Nonstiff Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1993)
邦訳: 「常微分方程式の数値解法Ⅰ 基礎編」スプリンガージャパン (2007)
- 使用例 (1)
- 次の2階常微分方程式の初期値問題を解く.
y1'' = -y1/R
y2'' = -y2/R
ただし, R = (y1^2 + y2^2)^(3/2)
(t = 0 において y1 = 0.5, y2 = 0, y1' = 0, y2' = √3)
Sub F2(N As Long, T As Double, Y() As Double, Yp2() As Double)
Dim R As Double
R = (Y(0) ^ 2 + Y(1) ^ 2) ^ (3 / 2)
Yp2(0) = -Y(0) / R
Yp2(1) = -Y(1) / R
End Sub
Sub Ex_Odex2_r()
Const N = 2
Dim T As Double, Y(N - 1) As Double, Yp(N - 1) As Double
Dim Tend As Double, Tout As Double
Dim RTol(0) As Double, ATol(0) As Double, Info As Long
Dim RCont() As Double, ICont() As Long, TT As Double
Dim YY(N - 1) As Double, YYpp(N - 1) As Double, Irtrn As Long, IRev As Long
RTol(0) = 0.0000000001 '1.0e-10
ATol(0) = RTol(0)
T = 0: Tend = 10: Y(0) = 0.5: Y(1) = 0: Yp(0) = 0: Yp(1) = Sqr(3)
Info = 0
Do
Tout = T + 1
IRev = 0
Do
Call Odex2_r(N, T, Y(), Yp(), Tout, RTol(), ATol(), RCont(), ICont(), Info, TT, YY(), YYpp(), Irtrn, IRev)
If IRev = 1 Then Call F2(N, TT, YY(), YYpp())
Loop While IRev <> 0
If Info <> 1 Then
Debug.Print "Error in Odex_r: Info =", Info
Exit Do
End If
Debug.Print T, Y(0), Y(1)
Loop While Tout < Tend
End Sub
- 実行結果
1 -0.427967245621361 0.863775701003976
2 -1.20572535244133 0.613566455365295
3 -1.4955436794428 8.16675371786414E-02
4 -1.33475968929777 -0.47684609244579
5 -0.700827262135424 -0.848381581815096
6 0.357480600846299 -0.44558418320875
7 -0.118067376945689 0.800372165255964
8 -1.03762003516232 0.730221556218178
9 -1.45981364707677 0.243039748770668
10 -1.42617024854319 -0.326583069133528
- 使用例 (2)
- 次の2階常微分方程式の初期値問題を解く(密出力を使用).
y1'' = -y1/R
y2'' = -y2/R
ただし, R = (y1^2 + y2^2)^(3/2)
(t = 0 において y1 = 0.5, y2 = 0, y1' = 0, y2' = √3)
Sub F2(N As Long, T As Double, Y() As Double, Yp2() As Double)
Dim R As Double
R = (Y(0) ^ 2 + Y(1) ^ 2) ^ (3 / 2)
Yp2(0) = -Y(0) / R
Yp2(1) = -Y(1) / R
End Sub
Sub Ex_Odex2_r_2()
Const N = 2, KM = 9
Dim T As Double, Y(N - 1) As Double, Yp(N - 1) As Double, Tend As Double
Dim RTol(0) As Double, ATol(0) As Double, Info As Long
Dim RCont((2 * KM + 6) * N - 1) As Double, ICont(N - 1) As Long, TT As Double
Dim YY(N - 1) As Double, YYpp(N - 1) As Double, Irtrn As Long, IRev As Long
Dim Iout As Long, Tout As Double, 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) = 0.5: Y(1) = 0: Yp(0) = 0: Yp(1) = Sqr(3)
Tout = 1
Info = 0
IRev = 0
Do
Call Odex2_r(N, T, Y(), Yp(), Tend, RTol(), ATol(), RCont(), ICont(), Info, TT, YY(), YYpp(), Irtrn, IRev, Iout)
If IRev = 1 Then
Call F2(N, TT, YY(), YYpp())
ElseIf IRev = 5 Then
While T >= Tout
Y0 = Contx2_r(0, Tout, RCont(), ICont())
Y1 = Contx2_r(1, Tout, RCont(), ICont())
Debug.Print Tout, Y0, Y1
Tout = Tout + 1
Wend
End If
Loop While IRev <> 0
If Info <> 1 Then Debug.Print "Error in Odex2_r: Info =", Info
End Sub
- 実行結果
1 -0.427967245557026 0.863775700969804
2 -1.20572535233059 0.613566455161914
3 -1.49554367919972 8.16675368902836E-02
4 -1.33475968872957 -0.476846092763801
5 -0.700827261091603 -0.848381581795345
6 0.357480601918728 -0.445584181833266
7 -0.118067378360736 0.800372165912016
8 -1.03762003630738 0.730221557185247
9 -1.45981364907474 0.243039751005974
10 -1.42617025292365 -0.326583065997367
|
1.9.6