|
|
◆ Doprin_r()
| Sub Doprin_r |
( |
N As |
Long, |
|
|
T As |
Double, |
|
|
Y() As |
Double, |
|
|
Yp() As |
Double, |
|
|
Tout As |
Double, |
|
|
RTol() As |
Double, |
|
|
ATol() As |
Double, |
|
|
Info As |
Long, |
|
|
TT As |
Double, |
|
|
YY() As |
Double, |
|
|
YYpp() 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 Safe As |
Double = 0, |
|
|
Optional Fac1 As |
Double = 0, |
|
|
Optional Fac2 As |
Double = 0, |
|
|
Optional Cnt As |
Long = 0 |
|
) |
| |
常微分方程式の初期値問題 (7(6)次ルンゲ・クッタ・ニュストレム法) (2階微分方程式用) (リバースコミュニケーション版)
- 目的
- 本ルーチンは2階の常微分方程式の初期値問題
d2y/dt2 = f(t, y), ただし t = t0 において y = y0, y' = y'0
本ルーチンはドルマンとプリンスによる7(6)次ルンゲ・クッタ・ニュストレム法のプログラムである. ステップ制御アルゴリズムを持っている.
詳細については下記参考文献を参照せよ.
Doprin_rはDoprinのリバースコミュニケーション版である.
- 引数
-
| [in] | N | 微分方程式の数. (N >= 1) |
| [in,out] | T | 本ルーチンはTからToutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される.
[in] 独立変数 t の初期値.
[out] ステップ終了時の最終値(通常Toutに等しい). 解がこの点まで正常に求められたことを示す. 新たなToutで再呼び出しすることにより, 新たな点まで積分を継続することができる. |
| [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 = 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)
= -3: パラメータ Y() の誤り.
= -4: パラメータ Yp() の誤り.
= -6: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0)
= -7: パラメータ ATol() の誤り. (ATol(i) < 0)
= -8: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1)
= -10: パラメータ YY() の誤り.
= -11: パラメータ YYpp() の誤り.
= -13: パラメータ IRev の誤り. (IRev != 0, 1, 5)
= 1: 正常終了.
= 11: 収束しなかった (ステップ数の最大値を超えたなど). |
| [out] | TT | IRev = 1: 再呼び出し時に2次微分係数値を求めて YYpp() に入れるべき点の T の値を返す.
IRev = 5: 現在の T の値. |
| [out] | YY() | 配列 YY(LYY - 1) (LYY >= N) IRev = 1: 再呼び出し時に2次微分係数値を求めて YYpp() に入れるべき点の Y() の値を返す.
IRev = 5: 現在の T における Y() の値. |
| [in,out] | YYpp() | 配列 YYpp(LYYpp - 1) (LYYpp >= N)
[in] IRev = 1: T (= TT) および Y (= YY()) における2次微分係数, すなわち, YYpp(i) = d2yi/dt2 = fi(TT, YY(0), ..., YY(N-1)) (i = 0 〜 N-1) を計算して, 再呼び出し時に入力する.
IRev = 5: 現在の T における Yp() の値. |
| [in,out] | IRev | リバースコミュニケーションの制御変数
[in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない.
[out] 0 以外の場合, 下記処理を行い IRev を変更せずに再び本ルーチンを呼び出すこと.
= 0: 処理終了. 正常終了かどうかは Info をチェックすること.
= 1: TT および YY() における2次微分係数の計算値を YYpp() に設定する. YYpp() 以外の変数を変更してはならない.
= 5: 途中結果を出力する. 通常はどの変数も変更してはならない(Iout を参照せよ). |
| [in] | Iout | (省略可)
計算の途中結果を出力するかどうか指定する. (省略時 = 0)
= 0: 出力しない. (IRev = 5では戻らない)
= 1: 途中結果出力のために成功したステップごとに IRev = 5 で戻る. リバースコミュニケーション版でない Doprin の Solout の呼び出しに相当する. t, y, yp の値が TT, YY(), YYpp() に返される.
(Ioutが上記以外の値であれば省略時の既定値とみなす) |
| [out] | Neval | (省略可)
関数評価回数. |
| [out] | Nstep | (省略可)
ステップ数. |
| [out] | Naccept | (省略可)
採用されたステップ数. |
| [out] | Nreject | (省略可)
不採用だったステップ数. (最初のステップはカウントされない) |
| [in] | Hmax | (省略可)
ステップ幅の最大値. (Hmax >= 0) (省略時 = Tout - T)
(Hmax = 0 であれば省略時の既定値とみなす) |
| [in] | MaxIter | (省略可)
許されるステップ数の最大値. (MaxIter >= 0) (省略時 = 2000)
(MaxIter = 0 であれば省略時の既定値とみなす) |
| [in] | Safe | (省略可)
ステップ幅推定時の安全係数. (0.0001 < Safe < 1) (省略時 = 0.9)
(Safe <= 0.0001 または Safe >= 1 であれば省略時の既定値とみなす) |
| [in] | Fac1 | (省略可) |
| [in] | Fac2 | (省略可)
ステップ幅選択パラメータ. (省略時: Fac1 = 0.2, Fac2 = 10)
Fac1 <= Hnew/Hold <= Fac2 となるようにステップ幅が選ばれる.
(Fac1=0 あるいは Fac2=0 であればそれぞれ省略時の既定値とみなす) |
| [in] | Cnt | (省略可)
Neval, Nstep, Naccept, Nrejectをリセットするタイミングを指定する. (省略時 = 0)
= 0: 本ルーチンが呼び出されるたびにリセットする.
<> 0: Info = 0 で呼び出されたときにだけリセットする. |
- 参考文献
- E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations I", Springer Series in Computational Mathematics, Springer-Verlag (1987)
- J. R. Dormand and P. J. Prince, "New Runge-Kutta algorithms for numerical simulation in dynamical astoronomy", Celestial Mechanics 18, 223?232 (1978)
- 使用例
- 次の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_Doprin_r()
Const N = 2
Dim T As Double, Y(N - 1) As Double, Yp(N - 1) As Double, Tend As Double, Tout As Double
Dim RTol(0) As Double, ATol(0) As Double, Info As Long, Irtrn As Long
Dim TT As Double, YY(N - 1) As Double, YYpp(N - 1) As Double, 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 Doprin_r(N, T, Y(), Yp(), Tout, RTol(), ATol(), 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 Doprin_r: Info =", Info
Exit Do
End If
Debug.Print T, Y(0), Y(1)
Loop While Tout < Tend
End Sub
- 実行結果
1 -0.427967245622589 0.863775701069405
2 -1.20572535251092 0.61356645560506
3 -1.49554367984035 8.16675377144076E-02
4 -1.33475969024023 -0.476846091796985
5 -0.700827263990535 -0.848381581520395
6 0.357480599005646 -0.445584185957294
7 -0.118067374412737 0.800372164683605
8 -1.03762003373262 0.730221558340758
9 -1.45981364767163 0.243039752399323
10 -1.42617025207063 -0.326583065015208
|
1.9.6