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

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  ) Initial value problem of ordinary differential equations (7(6)-th order Runge-Kutta-Nystrom 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. This routine is the code of the Runge-Kutta-Nystrom method of order 7(6) by Dormand and Prince with step size control. 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, Ypp() As Double) Ypp(i) = computed second derivative at given T and Y() (i = 0 to N-1) End Sub where N is the number of equations, and Ypp() is the computed second derivatives at given T and Y(), i.e. Ypp(i) = d2yi/dt2 = fi(T, Y(0), ..., Y(N-1)) (i = 0 to 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. [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 = 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) = -3: The argument Y() is invalid. = -4: The argument Yp() is invalid. = -6: The argument RTol() had an illegal value. (RTol(i) < 0, RTol(i) = 0 and ATol(i) = 0) = -7: The argument ATol() had an illegal value. (ATol(i) < 0) = -8: the argument Info had an illegal value (Info <> 0 and Info <> 1) = -10: The argument YY() is invalid. = -11: The argument YYpp() is invalid. = -13: the argument IRev had an illegal value. (IRev <> 0, 1 nor 5) = 1: Successful exit. = 11: Convergence failure (maximum number of steps exceeded, etc.). [out] TT IRev = 1: The value of T where the second derivative values should be evaluated and given in YYpp() in the next call. IRev = 5: The current value of T. [out] YY() Array YY(LYY - 1) (LYY >= N) IRev = 1: The value of Y where the second derivative values should be evaluated and given in YYpp() in the next call. IRev = 5: The value of Y() at current T. [in,out] YYpp() Array YYpp(LYYpp - 1) (LYYpp >= N) [in] IRev = 1: The computed second derivatives at given T (= TT) and Y (= YY()), i.e. YYpp(i) = d2yi/dt2 = fi(TT, YY(0), ..., YY(N-1)) (i = 0 to N-1), should be given in the next call. [out] IRev = 5: The value of Yp() at current T. [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. [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 second derivative values at TT and YY() in YYpp(). Do not alter any variables other than YYpp(). = 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. (Not 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 Doprin (non reverse communication version). The values t, y and yp are provided in TT, YY() and YYpp(). (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] Hmax (Optional) Maximal step size. (Hmax >= 0) (default = Tout - T) (If Hmax = 0, the default value will be used) [in] MaxIter (Optional) Maximum number of allowed steps. (MaxIter >= 0) (default = 2000) (If MaxIter = 0, the default value will be used) [in] Safe (Optional) The safety factor in step size prediction. (0.0001 < Safe < 1) (default = 0.9) (If Safe <= 0.0001 or Safe >= 1, the default value will be used) [in] Fac1 (Optional) [in] Fac2 (Optional) Parameters for step size selection. (default: Fac1 = 0.2, Fac2 = 10) The new step size is chosen subject to the restriction   Fac1 <= Hnew/Hold <= Fac2. (If Fac1 = 0 or Fac2 = 0, the default values 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 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) Example Program 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_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 Example Results 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