Contd5 Contd5_r Contd8 Contd8_r Contx1 Contx1_r Contx2 Contx2_r Derkf Derkf_r DerkfInt Dop853 Dop853_r Dopri5 Dopri5_r Doprin Doprin_r Dverk Dverk_r DverkInt Odex Odex2 Odex2_r Odex_r Retard Retard_r Ylag Ylag_r

Odex2() Sub Odex2 ( N As  Long, F2 As  LongPtr, T As  Double, Y() As  Double, Yp() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, Optional Solout2 As  LongPtr = NullPtr, 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  ) Initial value problem of ordinary differential equations (extrapolation method) (for second order differential equations) NOTE - THIS PROGRAM IS DEPRECATED AND WILL BE REMOVED IN THE NEXT VERSION. Purpose This routine integrates a system of second order ordinary differential equations of the form d2y/dt2 = f(t, y), y = y0, y' = y'0 at t = t0 where t0, y0 and y'0 are the given initial values of t, y and y', respectively. y and y' may be a vector if the above is a system of differential equations. Odex2 is an extrapolation algorithm code, based on the Stoermer's rule with step size control, order selection and dense output. See for details in the reference below. Parameters [in] N Number of differential equations. (N >= 1) [in] F2 The user supplied subroutine, which calculates the derivatives of the differential equations, defined as follows. Sub F2(N As Long, T As Double, Y() As Double, Yp2() As Double) Yp2(i) = computed second order derivative at given T and Y() (i = 0 to N-1) End Sub where N is the number of equations, and Yp2() is the computed second order derivatives at given T and Y(), i.e. Yp2(i) = d2yi/dt2 = fi(T, Y(0), ..., Y(N-1)) (i = 0 to N-1). The other variables than Yp2() should not be altered. [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,out] Yp() Array Yp(LYp - 1) (LYp >= N) [in] Initial derivative values y'1, ..., y'n at initial T. [out] Computed approximations of derivatives at last T (normally euqals 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 = 2N). If LRTol = 2, ..., 2N-1, LRTol = 1 is assumed. If LRTol > 2N, LRTol = 2N is assumed. Even if LRTol = 2N, 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 to roughly satisfy the following conditions for each element of Y() and Yp().   In the case of scalar：     abs(local error of Y(i)) <= RTol*abs(Y(i)) + ATol     abs(local error of Yp(i)) <= RTol*abs(Y(i)) + ATol   In the case of vector (i = 0 〜 N-1)：     abs(local error of Y(i)) <= RTol(i)*abs(Y(i)) + ATol(i)     abs(local error of Yp(i)) <= RTol(i+N)*abs(Y(i)) + ATol(i+N) 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 = 2N). If LATol = 2, ..., 2N-1, LATol = 1 is assumed. If LATol > 2N, LATol = 2N is assumed. Even if LATol = 2N, 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 to roughly satisfy the following conditions for each element of Y() and Yp().   In the case of scalar：     abs(local error of Y(i)) <= RTol*abs(Y(i)) + ATol     abs(local error of Yp(i)) <= RTol*abs(Y(i)) + ATol   In the case of vector (i = 0 〜 N-1)：     abs(local error of Y(i)) <= RTol(i)*abs(Y(i)) + ATol(i)     abs(local error of Yp(i)) <= RTol(i+N)*abs(Y(i)) + ATol(i+N) 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) = -4: The argument Y() is invalid. = -5: The argument Yp() 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) = 1: Successful exit. = 2: Interrupted by Solout2 (normal return). = 11: Maximum number of steps exceeded. [in] Solout2 (Optional) The user supplied subroutine to print out the intermediate solutions, which is called after every successful step, defined as follows. (default = NullPtr) Sub Solout2(Nr As Long, Told As Double, T As Double, Y() As Double, Yp() As Double, N As Long, RCont As Double, ICont As Long, Irtrn As Long) Output the Y() and Yp() 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 Solout2, respectively. Irtrn also serves to interrupt the integration. If Irtrn is set to the negative value in Solout2, the integration will be interrupted and exit with Info = 2. If the numerical solution is altered in Solout2, set Irtrn = 3. Dense output is supported by the control information RCont and ICont. 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) = Contx2(i, T2, RCont, ICont) End Sub If Solout2 is not provided (if Solout2 = NullPtr), the intermediate solutions will not be printed out. [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 Solout2 is not used, 4 (if Solout2 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, ... 1 to 3 cannot be specified if Solout2 is used. (For other values or 1 to 3 when Solout2 is used, 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 Solout2 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 <= 8) (default = 6) (If Mudif < 1 or Mudif > 8, 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 second order ordinary differential equations. y1'' = -y1/R y2'' = -y2/R where R = (y1^2 + y2^2)^(3/2) (y1 = 0.5, y2 = 0, y1' = 0, y2' = √3 at t = 0) 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() 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 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 Call Odex2(N, AddressOf F2, T, Y(), Yp(), Tout, RTol(), ATol(), Info) If Info <> 1 Then Debug.Print "Error in Odex2: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub Example Results 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 Example Program (2) Solve the following initial value problem of second order ordinary differential equations (using dense output). y1'' = -y1/R y2'' = -y2/R where R = (y1^2 + y2^2)^(3/2) (y1 = 0.5, y2 = 0, y1' = 0, y2' = √3 at t = 0) 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_2() Const N = 2 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 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 Call Odex2(N, AddressOf F2, T, Y(), Yp(), Tend, RTol(), ATol(), Info, AddressOf SoloutX2) If Info <> 1 Then Debug.Print "Error in Odex2: Info =", Info End Sub Sub SoloutX2(Nr As Long, Told As Double, T As Double, Y() As Double, Yp() 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 = Contx2(0, Tout, RCont, ICont) Y1 = Contx2(1, Tout, RCont, ICont) Debug.Print Tout, Y0, Y1 Tout = Tout + 1 Wend End Sub Example Results 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