|
|
◆ Rodas()
| Sub Rodas |
( |
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 Dfx As |
LongPtr = NullPtr, |
|
|
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 Meth As |
Long, |
|
|
Optional Pred As |
Long, |
|
|
Optional Safe As |
Double, |
|
|
Optional Fac1 As |
Double, |
|
|
Optional Fac2 As |
Double, |
|
|
Optional Cnt As |
Long |
|
) |
| |
Initial value problem of ordinary differential equations (4(3)-th order Rosenbrock method)
- Purpose
- This routine computes a numerical solution of a stiff (or differential algebraic) system of first order ordinary differential equations of the form
M * dy/dt = f(t, y), y = y0 at t = t0
Rodas is the code based on the 4(3)-th order Rosenbrock method. It is provided with the step control algorithm and the continuous output feature.
See for details in the reference below.
- Parameters
-
| [in] | N | Number of differential equations. (N >= 1) |
| [in] | F | The user supplied subroutine, which calculates the derivatives of the differential equations, defined as follows. Sub F(N As Long, T As Double, Y() As Double, Yp() As Double)
Yp(i) = computed derivative at given T and Y() (i = 0 to N-1)
End Sub
In many cases, M is an identity matrix. However, it can be defined in subroutine Mas if necessary. |
| [in] | Ifcn | Whether f() depends on t or not.
= 0: Independent of t: dy/dt = f(y) (autonomous).
= 1: May depend on t: dy/dt = f(y, t) (non-autonomous). |
| [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] | 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)
= -5: The argument Y() is invalid.
= -7: The argument RTol() had an illegal value. (RTol(i) < 0, RTol(i) = 0 and ATol(i) = 0)
= -8: The argument ATol() had an illegal value. (ATol(i) < 0)
= -9: the argument Info had an illegal value (Info <> 0 and Info <> 1)
= -16: The argument Mljac or Mlmas had an illegal value. (Mlmas > Mljac)
= -17: The argument Mujac or Mumas had an illegal value. (Mumas > Mujac)
= -23: The argument M1 had an illegal value. (M1 < 0)
= -24: The argument M1 or M2 had an illegal value. (M2 < 0 or M1 + M2 > N)
= 1: Successful exit.
= 2: Interrupted by Solout (normal return).
= 11: Maximum number of steps exceeded.
= 12: Step size becomes too small.
= 13: Matrix is repeatedly singular. |
| [in] | Solout | (Optional)
The user supplied subroutine to print out the intermediate solutions, which is called after every successful step, defined as follows. (default = NullPtr) Sub Solout(Nr As Long, Told As Double, T As Double, Y() As Double, N As Long, Cont As Double, Irtrn As Long)
Output the Y() values at Nr-th step T.
Told is the previous value of T. N is the order of equations.
The value of Irtrn will be 0, 1 or 2 in the first, intermediate or last call of Solout, respectively.
Irtrn also serves to interrupt the integration. If Irtrn is set to the negative value in Solout, the integration will be interrupted and exit with Info = 2.
Dense output is supported by the control information Cont.
The solution Y(i) (0 <= i <= N-1) at the arbitrary point T2 in the interval [Told, T] can be computed by the function call
Y(i) = Contro(i, T2, Cont)
End Sub
|
| [in] | Jac | (Optional)
The user supplied subroutine, which computes the Jacobian dfi(t, y)/dyj, defined as follows. (default = NullPtr) Sub Jac(N As Long, T As Double, Y() As Double, Ypd() As Double)
Ypd(i,j) = the calculated value of dfi/dyj at T and Y() (for i = 0 to N-1, j = 0 to N-1).
End Sub
If Jac is not provided (if Jac = NullPtr), the Jacobian will be computed by finite differences. |
| [in] | Mljac | (Optional)
The lower bandwidth of Jacobian. (0 <= Mljac <= N) (default = N)
If Mljac = N, Jacobian is stored as N x N full matrix. If Mljac < N, Jacobian is stored in band matrix form.
(If Mljac < 0 or Mljac > N, the default value will be used) |
| [in] | Mujac | (Optional)
The upper bandwidth of Jacobian. (0 <= Mujac <= N) (default = 0)
If Mljac = N, Mujac is ignored.
(If Mujac < 0 or Mujac > N, the default value will be used) |
| [in] | Dfx | (Optional)
The user supplied subroutine, which computes the partial derivatives of f(t, y) with respect to t, defined as follows. Sub Dfx(N As Long, T As Double, Y() As Double, Fx() As Double)
Fx(i) = dfi/dt (for i = 0 to N-1)
End Sub
If Dfx is not provided (if Dfx = NullPtr), the partial derivatives will be computed by finite differences. |
| [in] | Mas | (Optional)
The user supplied subroutine, which provides the mass matrix M, defined as follows (default = NullPtr). Sub Mas(N As Long, Am() As Double)
Am(i,j) = the value of mass matrix Mij (i = 0 to N-1, j = 0 to N-1).
End Sub
If Mas is not provided (if Mas = NullPtr), the mass matrix M is supposed to be the identity matrix. |
| [in] | Mlmas | (Optional)
The lower bandwidth of mass matrix M. (0 <= Mlmas <= N) (default = N)
If Mlmas = N, M is stored as N x N full matrix. If Mlmas < N, M is stored in band matrix form.
(If Mlmas < 0 or Mlmas > N, the default value will be used) |
| [in] | Mumas | (Optional)
The upper bandwidth of mass matrix M. (0 <= Mumas <= N) (default = 0)
If Mlmas = N, Mumas is ignored.
(If Mumas < 0 or Mumas > N, the default value will be used) |
| [out] | Neval | (Optional)
Number of function evaluations. (Those for Jacobian evaluations are not included) |
| [out] | Njac | (Optional)
Number of Jacobian evaluations. (Those by finite differences are included) |
| [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] | M1,M2 | (Optional)
If the first M1 equations has the following form
y'(i) = y(i + M2) for i = 1 to M1,
with M1 a multiple of M2, and the remaining equations do not explicitly depend on y'(M1), ..., y'(N-1), efficient computation can be achieved by setting parameters M1 and M2 to nonzero values. (M1 > 0, M2 > 0, M1 + M2 <= N) (default M1 = M2 = 0)
When parameters are set to nonzero, only the elements of non-trivial part of the Jacobian (rows M1+1 to N) heve to be stored in (N - M1) x N array in Jac. Also only the elements of right lower block of order N - M1 of the mass matrix M have to be stored in (N - M1) x (N - M1) array in Mas. |
| [in] | Hinit | (Optional)
Initial step size. (default = 1.0e-6)
H = 1/||f'||, usually 1.0e-3 or 1.0e-5 is good for stiff equations with initial transient.
(If Hinit = 0, 1.0e-6 will be used) |
| [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] | Meth | (Optional)
Choice of the coefficients. (default = 1)
= 1: Method in the reference on page 452.
= 2: Same method with different parameters.
= 3: Method with coefficients of Gerd Steinebach.
(For other values, the default value will be used) |
| [in] | Pred | (Optional)
Switch for step size strategy. (default = 1)
= 1: Model predictive controller (Gustafsson).
= 2: Classical step size control.
(For other values, the default value will be assumed) |
| [in] | Safe | (Optional)
The safety factor in step size prediction. (default = 0.9)
(If Safe <= 0.001 or Safe >= 1, the default value will be used) |
| [in] | Fac1,Fac2 | (Optional)
Parameters for step size selection. (Fac1 <= 1, Fac2 >= 1) (default: Fac1 = 0.2, Fac2 = 6)
The new step size is chosen subject to the restriction Fac1 < Hnew/Hold < Fac2.
(If Fac1 = 0 or Fac1 > 1, Fac2 = 0 or Fac2 < 1, the default values will be used) |
| [in] | Cnt | (Optional)
Specifies when Neval, Njac, 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 II. Stiff and differential-algebraic Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1996)
- Example Program (1)
- Solve the following initial value problem of ordinary differential equations (stiff problem).
dy1/dt = -2*y1 + y2 - cos(t)
dy2/dt = 1998*y1 - 1999*y2 + 1999*cos(t) - sin(t)
(y1 = 1, y2 = 2 at t = 0)
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_Rodas()
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 Rodas(N, AddressOf F2, Ifcn, T, Y(), Tout, RTol(), ATol(), Info)
If Info <> 1 Then
Debug.Print "Error in Rodas: Info =", Info
Exit Do
End If
Debug.Print T, Y(0), Y(1)
Loop While Tout < Tend
End Sub
- Example Results
1 0.367879441171443 0.908181747041726
2 0.135335283236614 -0.280811553312466
3 4.97870683678661E-02 -0.940205428235504
4 1.83156388887351E-02 -0.635327981976607
5 6.73794699908466E-03 0.290400132463827
6 2.47875217666488E-03 0.962649038829728
7 9.11881965553467E-04 0.754814136310732
8 3.35462627903455E-04 -0.145164571182498
9 1.2340980408788E-04 -0.911006852082675
10 4.53999297639715E-05 -0.839026129149385
- Example Program (2)
- Solve the following initial value problem of ordinary differential equations (stiff problem) (using dense output).
dy1/dt = -2*y1 + y2 - cos(t)
dy2/dt = 1998*y1 - 1999*y2 + 1999*cos(t) - sin(t)
(y1 = 1, y2 = 2 at t = 0)
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_Rodas_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 Rodas(N, AddressOf F2, Ifcn, T, Y(), Tend, RTol(), ATol(), Info, AddressOf SoloutRo)
If Info <> 1 Then Debug.Print "Error in Rodas: Info =", Info
End Sub
Sub SoloutRo(Nr As Long, Told As Double, T As Double, Y() As Double, N As Long, Cont As Double, Irtrn As Long)
Dim Y0 As Double, Y1 As Double
Static Tout As Double
If Nr = 1 Then Tout = 1
While T >= Tout
Y0 = Contro(0, Tout, Cont)
Y1 = Contro(1, Tout, Cont)
Debug.Print Tout, Y0, Y1
Tout = Tout + 1
Wend
End Sub
- Example Results
1 0.36787944117149 0.908181746948026
2 0.135335283236613 -0.280811553312504
3 4.97870683678624E-02 -0.940205428229929
4 1.83156388886953E-02 -0.635327981896935
5 6.73794699908206E-03 0.290400132468997
6 2.47875217666092E-03 0.962649038837632
7 9.11881965597003E-04 0.754814136223758
8 3.35462627893874E-04 -0.145164571163287
9 1.23409804039225E-04 -0.911006851985428
10 4.53999297635988E-05 -0.839026129148628
|
1.9.6