|
|
◆ Odex_r()
| Sub Odex_r |
( |
N As |
Long, |
|
|
T As |
Double, |
|
|
Y() 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, |
|
|
YYp() 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 Mstab As |
Long = 0, |
|
|
Optional Jstab 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 |
|
) |
| |
Initial value problem of ordinary differential equations (extrapolation method (GBS algorithm)) (reverse communication version)
- Purpose
- This routine integrates a system of first order ordinary differential equations of the form
dy/dt = f(t, y), y = y0 at t = t0
Odex is an extrapolation algorithm (GBS) code, based on the explicit midpoint rule with step size control, order selection and dense output.
See for details in the reference below.
Odex_r is the reverse communication version of Odex.
- Parameters
-
| [in] | N | Number of differential equations. (N >= 1) |
| [in,out] | T | This routine integrates from T to Tout. The initial point of the integration is to be given, and the last point of the final step will be returned.
[in] Initial value of the independent variable T.
[out] Last value of the independent variable T of the final step (normally equals to Tout). The solution was successfully advanced to this point. It is possible to continue the integration to new point by recalling this routine with the new Tout value with setting Info = 1. |
| [in,out] | Y() | Array Y(LY - 1) (LY >= N)
[in] Initial values of the dependent variables Y() at initial T.
[out] Computed solution approximation at last T (normally equals to Tout). |
| [in] | Tout | Set Tout to the point at which a solution is desired. Integration either forward in T (Tout > T) or backward in T (Tout < T) is permitted.
The routine advances the solution from T to Tout using step sizes which are automatically selected so as to achieve the desired accuracy. |
| [in] | RTol() | Array RTol(LRTol - 1) (LRTol >= 1) (all components of RTol() >= 0)
The relative error tolerance(s) to tell the code how accurately you want the solution to be computed. This parameter may be a scalar (LRTol = 1) or a vector (LRTol = N). If LRTol = 2, ... or N-1, LRTol = 1 is assumed. If LRTol > N, LRTol = N is assumed. Even if LRTol = N, it is assumed to be 1 if LATol = 1.
The tolerances are used by the code in a local error test at each step which requires roughly that
abs(local error of Y(i)) <= RTol(i)*abs(Y(i)) + ATol(i)
for each component of Y() (i = 0 to LRTol-1).
Setting RTol(i) = 0 results in a pure absolute error test on that component. RTol(i) and ATol(i) should not be zero at the same time (i = 0 to LRTol-1). |
| [in] | ATol() | Array ATol(LATol - 1) (LATol >= 1) (all components of ATol() >= 0)
The absolute error tolerance(s) to tell the code how accurately you want the solution to be computed. This parameter may be a scalar (LATol = 1) or a vector (LATol = N). If LATol = 2, ... or N-1, LATol = 1 is assumed. If LATol > N, LATol = N is assumed. Even if LATol = N, it is assumed to be 1 if LRTol = 1.
The tolerances are used by the code in a local error test at each step which requires roughly that
abs(local error of Y(i)) <= RTol(i)*abs(Y(i)) + ATol(i)
for each component of Y() (i = 0 to LATol-1).
Setting ATol(i) = 0 results in a pure relative error test on that component. RTol(i) and ATol(i) should not be zero at the same time (i = 0 to LRTol-1). |
| [in,out] | RCont() | Array RCont(LRCont - 1) (LRCont >= (2*Km + 5)*N)
Control information for dense output.
(Not referenced if Iout = 0) |
| [in,out] | ICont() | Array ICont(LICont - 1) (LICont >= N)
Control information for dense output.
(Not referenced if Iout = 0) |
| [in,out] | Info | [in]
= 0: Initialize and start computation (Solve new problem).
= 1: Continue computation with new Tout value (Resume computation of previous call).
[out]
= -1: The argument N had an illegal value. (N < 1)
= -3: The argument Y() is invalid.
= -5: The argument RTol() had an illegal value. (RTol(i) < 0, RTol(i) = 0 and ATol(i) = 0)
= -6: The argument ATol() had an illegal value. (ATol(i) < 0)
= -7: The argument RCont() is invalid.
= -8: The argument ICont() is invalid.
= -9: the argument Info had an illegal value (Info <> 0 and Info <> 1)
= -11: The argument YY() is invalid.
= -12: The argument YYp() is invalid.
= -14: the argument IRev had an illegal value. (IRev <> 0, 1 nor 5)
= 1: Successful exit.
= 2: Interrupted by Irtrn (normal return).
= 11: Maximum number of steps exceeded. |
| [out] | TT | IRev = 1: The value of T where the derivative values should be evaluated and given in YYp() in the next call. |
| [out] | YY() | Array YY(LYY - 1) (LYY >= N)
IRev = 1: The value of Y where the derivative values should be evaluated and given in YYp() in the next call. |
| [in] | YYp() | Array YYp(LYYp - 1) (LYYp >= N)
IRev = 1: The computed derivatives at given T (= TT) and Y (= YY()), i.e. YYp(i) = dyi/dt = fi(TT, YY(0), ..., YY(N-1)) (i = 0 to N-1), should be given in the next call. |
| [in,out] | Irtrn | [in] IRev = 5: Do not alter except for the following special cases.
If Irtrn is set to the negative value, the integration will be interrupted and exit with Info = 2. If the numerical solution is altered in Solout, set Irtrn = 3.
[out] IRev = 5: Returns 0, 1 or 2 in the first, intermediate or last return with IRev = 5, respectively. |
| [in,out] | IRev | Control variable for reverse communication.
[in] Before first call, IRev should be initialized to zero. On succeeding calls, IRev should not be altered.
[out] If IRev is not zero, complete the following tasks and call this routine again without changing IRev.
= 0: Computation finished. See return code in Info.
= 1: User should set the computed derivative values at TT and YY() in YYp(). Do not alter any variables other than YYp().
= 5: User may output the intermediate result. Do not alter any variables. (See Iout) |
| [in] | Iout | (Optional)
Specifies if the intermediate result output is required. (default = 0)
= 0: Output is not required. (Never return with IRev = 5)
= 1: Returns after every successful step with IRev = 5 to output the intermediate results. This is same as calling Solout in the case of normal Odex (non reverse communication version). The corresponding information is as follows.
Nr = Naccept + 1, Told = previous T, T = current T, Y() = current Y().
RCont() and ICont() are used in the same way with Solout for the dense output.
Y(i) = Contx1_r(i, T2, RCont(), ICont())
(If other value is specified for Iout, Iout = 0 will be assumed) |
| [out] | Neval | (Optional)
Number of function evaluations. |
| [out] | Nstep | (Optional)
Number of computed steps. |
| [out] | Naccept | (Optional)
Number of accepted steps. |
| [out] | Nreject | (Optional)
Number of rejected steps. (Step rejections in the first step are not counted) |
| [in] | Hinit | (Optional)
Initial step size. (default = 0.0001)
H = 1/||f'||, usually 0.1 or 0.001 is good for initial value.
(If |Hinit| < 0.0001, Hinit = 0.0001 is assumed) |
| [in] | Hmax | (Optional)
Maximal step size. (default = Tout - T)
(If Hmax = 0, the default value will be used) |
| [in] | MaxIter | (Optional)
Maximum number of allowed steps. (default = 100000)
(If MaxIter <= 0, the default value will be used) |
| [in] | Km | (Optional)
The maximum number of columns in the extrapolation table. (Km >= 3) (default = 9)
(If Km < 3, the default value will be used) |
| [in] | Nsequ | (Optional)
Switch for the step size sequence. (default = 1 (if Solout is not used, 4 (if Solout is used))
= 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, ...
= 5: 4, 8, 12, 16, 20, 24, 28, 32, ...
1 to 3 cannot be specified if Solout is used.
(For other values or 1 to 3 when Solout is used, the default value will be used) |
| [in] | Mstab | (Optional)
Stability check is avtivated at most Mstab times in one line of the extrapolation table. (default = 1)
(If Mstab <= 0, the default value will be used) |
| [in] | Jstab | (Optional)
Stability check is avtivated only in the lines 1 to Jstab of the extrapolation table. (default = 2)
(If Jstab <= 0, the default value will be used) |
| [in] | Iderr | (Optional)
Switch for error estimate in the dense output formula. (default = 0)
= 0: Activated.
= 1: Not activated.
If Solout is not used, error estimation will not be activated regradless of Iderr.
(For other values, Iderr = 1 is assumed) |
| [in] | Mudif | (Optional)
Parameter to determine the degree of interpolation formula. (1 <= Mudif <= 6) (default = 4)
(If Mudif < 1 or Mudif > 6, the default value will be used) |
| [in] | Fac1 | (Optional) |
| [in] | Fac2 | (Optional)
Parameters for step size selection. (default: Fac1 = 0.02, Fac2 = 4)
The new step size for the j-th diagonal entry is chosen subject to the restriction
Facmin/Fac2 <= j-th Hnew/Hold <= 1/Facmin
where Facmin = Fac1^(1/(2*j - 1)).
(If Fac1 = 0 or Fac2 = 0, the default values will be used respectively) |
| [in] | Fac3 | (Optional) |
| [in] | Fac3 | (Optional)
Parameters for the order selection (default: Fac3 = 0.8, Fac4 = 0.9)
Order is decreased if W(k-1) <= W(k)*Fac3.
Order is increased if W(k) <= W(k-1)*Fac4.
W(k) = A(k)/H(k) is the work per unit step. (A(k) is the number of function evaluations, H(k) is the step size)
(If Fac3 = 0 or Fac4 = 0, the default values will be used respectively) |
| [in] | Safe1 | (Optional) |
| [in] | Safe2 | (Optional)
Safety factors for step control algorithm. (default: Safe1 = 0.65, Safe2 = 0.94)
(If Safe1 = 0 or Safe2 = 0, the default values will be used) |
| [in] | Safe3 | (Optional)
Step size is reduced by this factor if the stability check is negative (default = 0.5)
(If Safe3 <= 0, the default value will be used. if Safe3 > 1, Safe3 = 1 will be used) |
| [in] | Cnt | (Optional)
Specifies when Neval, Nstep, Naccept and Nreject are reset to zero. (default = 0)
= 0: Reset whenever this routine is called.
<> 0: Reset only if this routine is called with Info = 0. |
- Reference
- E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations. Nonstiff Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1993)
- Example Program (1)
- Solve the following initial value problem of ordinary differential equations.
dy1/dt = -2*y1 + y2 - cos(t)
dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t)
(y1 = 1, y2 = 2 at t = 0)
Sub F1(N As Long, T As Double, Y() As Double, Yp() As Double)
Yp(0) = -2 * Y(0) + Y(1) - Cos(T)
Yp(1) = 2 * Y(0) - 3 * Y(1) + 3 * Cos(T) - Sin(T)
End Sub
Sub Ex_Odex_r()
Const N = 2
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 RCont() As Double, ICont() As Long, TT As Double
Dim YY(N - 1) As Double, YYp(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) = 1: Y(1) = 2
Info = 0
Do
Tout = T + 1
IRev = 0
Do
Call Odex_r(N, T, Y(), Tout, RTol(), ATol(), RCont(), ICont(), Info, TT, YY(), YYp(), Irtrn, IRev)
If IRev = 1 Then Call F1(N, TT, YY(), YYp())
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
- Example Results
1 0.367879441180513 0.90818174702144
2 0.135335283238428 -0.280811553314161
3 4.97870683602588E-02 -0.940205428217376
4 1.83156388764711E-02 -0.635327981950352
5 6.73794699120742E-03 0.290400132478066
6 2.47875218511396E-03 0.962649038810137
7 9.11881969599684E-04 0.754814136300769
8 3.35462632178912E-04 -0.145164571189262
9 1.23409800123739E-04 -0.911006852072663
10 4.53999253959451E-05 -0.839026129137953
- Example Program (2)
- Solve the following initial value problem of ordinary differential equations (using dense output).
dy1/dt = -2*y1 + y2 - cos(t)
dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t)
(y1 = 1, y2 = 2 at t = 0)
Sub F1(N As Long, T As Double, Y() As Double, Yp() As Double)
Yp(0) = -2 * Y(0) + Y(1) - Cos(T)
Yp(1) = 2 * Y(0) - 3 * Y(1) + 3 * Cos(T) - Sin(T)
End Sub
Sub Ex_Odex_r_2()
Const N = 2, KM = 9
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 RCont((2 * KM + 5) * N - 1) As Double, ICont(N - 1) As Long, TT As Double
Dim YY(N - 1) As Double, YYp(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) = 1: Y(1) = 2
Tout = 1
Info = 0
IRev = 0
Do
Call Odex_r(N, T, Y(), Tend, RTol(), ATol(), RCont(), ICont(), Info, TT, YY(), YYp(), Irtrn, IRev, Iout)
If IRev = 1 Then
Call F1(N, TT, YY(), YYp())
ElseIf IRev = 5 Then
While T >= Tout
Y0 = Contx1_r(0, Tout, RCont(), ICont())
Y1 = Contx1_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 Odex_r: Info =", Info
End Sub
- Example Results
1 0.367879441157644 0.908181747067181
2 0.135335283236881 -0.280811553311066
3 4.97870683669239E-02 -0.940205428230704
4 1.83156389169758E-02 -0.635327982031361
5 6.73794699664751E-03 0.290400132467187
6 2.47875227629557E-03 0.962649038627776
7 9.11881939436941E-04 0.754814136361094
8 3.35462629214848E-04 -0.145164571183337
9 1.23409825870176E-04 -0.911006852124154
10 4.53998495624352E-05 -0.839026128986288
|
1.9.6