Contd5 Contd5_r Contd8 Contd8_r Contx1 Contx1_r Contx2 Contx2_r Deabm Deabm_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

Retard_r() Sub Retard_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Ngrid As  Long, TGrid() 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 Nrdens As  Long = 0, Optional Hinit As  Double = 0, Optional Hmax As  Double = 0, Optional MaxIter As  Long = 0, Optional Nstiff As  Long = 0, Optional Safe As  Double = 0, Optional Fac1 As  Double = 0, Optional Fac2 As  Double = 0, Optional Beta As  Double = 0, Optional Cnt As  Long = 0  ) Initial value problem of delay differential equations (5(4)-th order Dorman-Prince method) (reverse communication version) Purpose This routine integrates a system of first order delay ordinary differential equations of the form dy/dt = f(t, y(t), y(t-τ), ...), y = y0 at t = t0 where t0 and y0 are the given initial values of t and y, respectively. y may be a vector if the above is a system of differential equations. Retard is the explicit Runge-Kutta code based on the 5(4)-th order Dormand-Prince method. It is provided with the step control algorithm and the dense output feature. See for details in the reference below. Retard_r is the reverse communication version of Retard. 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] Ngrid Number of grid points specified by Tgrid(). (Ngrid >= 1) The solution may have the discontinuities in its derivatives (grid points). In such case, by specifying the points, the precision and the computation time can be improved. The final point, at which the solution is computed, is included. That is, Ngrid = 1 even if no point of discontinuity is specified. (If Ngrid <= 0, Ngrid = 0 is assumed) [in] Tgrid() Array Tgrid(LTgrid - 1) (LTgrid >= Ngrid) Values of t at grid points should be set to Tgrid(0), Tgrid(1), ..., Tgrid(Ngrid - 1). The final point (= Tout), at which the solution is computed, should be given in the last component Tgrid(Ngrid - 1). If no point of discontinuity is specified (Ngrid = 1), Tgrid(0) is automatically set to Tout. [out] RCont() Array RCont(LRCont - 1) (LRCont >= Mxst*(5*Nrdens) Control information for back-values of the solution. Mxst is the number of back steps to be stored (back steps available through Ylag). If RTol and ATol are not very small, typical LRCont value can be decided to satisfy Mxst = 500. [out] ICont() Array ICont(LICont - 1) (LICont >= Nrdens) Control information for back-values of the solution. [in,out] Info [in] = 0: Initialize and start computation (Solve new problem). = 1: Continue computation with new Tout value (Resume computation of previous call). Valid only when Ngrid = 1. [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) = -8: The argument Tgrid() is invalid. = -9: The argument RCont() is invalid. = -10: The argument ICont() is invalid. = -11: the argument Info had an illegal value. (Info <> 0 and Info <> 1) = -16: 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. = 12: Step size becomes too small. = 13: Problem is probably stiff (interrupted). = 14: Computation interrupted by Ylag. (RTol and ATol may be too small, or Mxst may be small) [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 Irtrn unless user want to interrupt the integration. 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: YYp() should be set to the derivatives at given TT and YY(), i.e. Yp(i) = dyi/dt = fi(t, y(t), y(t-τ), ... ) (i = 0 to N-1). The other variables should not be altered. Ylag = y(t - τ) can be obtained by using Ylag_r as follows. IRev = 0 Do Call Ylag_r(I, T, RCont(), ICont(), YY, IRev) If IRev = 1 Then YY = value of Y(I) at T Loop While IRev <> 0 where I is the component number (0 to N - 1) and T is t - τ. RCont() and ICont() must be passed without changing those values. The result is returned in YY when IRev = 0. If t0 - τ <= T <= t0 (where t0 is the initial value of t), IRev = 1 is returned. Then the value of Y(I) at T must be given in YY. = 5: User may output the intermediate result. Do not alter any variables. (See Iout) [in] Iout (Optional) Whether the intermediate result output is required. (default = 0) = 0: Output is not required. (Not return with IRev = 5) = 1: Output is required after every successful step (Return with IRev = 5). This is same as calling Solout in the case of normal Retard (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 when calling Ylag_r for the dense output. (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] Nrdens (Optional) Number of components to be stored (elements available through Ylag). (default = N) Components 0 to Nrdens-1 will be stored for Ylag. (If Nrdens <= 0 or Nrdens > N, the default value will be used) [in] Hinit (Optional) Initial step size. (default = to be estimated from the initial function values) (If Hinit = 0, the default value 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] Nstiff (Optional) Test for stiffness is activated after every Nstiff steps. (default = 1000) If Nstiff < 0, the stiffness test is not activated. (If Nstiff = 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] Beta (Optional) The parameter for the stabilized step size control. (Beta <= 0.2) (default = 0.04) (If Beta < 0 or Beta > 0.2, 0 is assumed) [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 delay differential equations. y'1(t) = -y1(t)*y2(t - 1) + y2(t - 10) y'2(t) = y1(t)*y2(t - 1) - y2(t) y'3(t) = y2(t) - y2(t - 10) (y1(t) = 5, y2(t) = 0.1, y3(t) = 1 at t <= 0) (discontinuous at t = 1, 2, ..., 10, 20 and 40) Sub Ex_Retard_r() Const N = 3, NTout = 12, Ngrid = 0 Const Mxst = 500, LRCont = Mxst * (5 * N + 2), LICont = N Dim T As Double, Y(N - 1) As Double Dim Tout(NTout - 1) As Double, Tgrid(0) As Double Dim RTol(0) As Double, ATol(0) As Double, ITol As Long, Info As Long Dim RCont(LRCont - 1) As Double, ICont(LICont - 1) As Long, TT As Double Dim YY(N - 1) As Double, YYp(N - 1) As Double, Irtrn As Long, IRev As Long Dim Y2L1 As Double, Y2L10 As Double, I As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) T = 0: Y(0) = 5: Y(1) = 0.1: Y(2) = 1 For I = 0 To 9 Tout(I) = I + 1 Next Tout(10) = 20 Tout(11) = 40 Info = 0 For I = 0 To NTout - 1 IRev = 0 Do Call Retard_r(N, T, Y(), Tout(I), RTol(), ATol(), Ngrid, Tgrid(), RCont(), ICont(), Info, TT, YY(), YYp(), Irtrn, IRev) If IRev = 1 Then Y2L1 = YlagR(1, TT - 1, RCont(), ICont()) Y2L10 = YlagR(1, TT - 10, RCont(), ICont()) YYp(0) = -YY(0) * Y2L1 + Y2L10 YYp(1) = YY(0) * Y2L1 - YY(1) YYp(2) = YY(1) - Y2L10 End If Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Retard_r: Info =", Info Exit For End If Debug.Print T, Y(0), Y(1), Y(2) Next End Sub Function YlagR(II As Long, T As Double, RCont() As Double, ICont() As Long) As Double Dim YY As Double, IRev As Long IRev = 0 Do Call Ylag_r(II, T, RCont(), ICont(), YY, IRev) If IRev = 1 Then YY = 0.1 Loop While IRev <> 0 YlagR = YY End Function Example Results 1 4.61934967214384 0.338647989701885 1.14200233815428 2 3.71355732102849 0.798576426134181 1.58786625283732 3 2.22539076042276 1.33022795319548 2.54438128638175 4 0.832713255497743 1.39799979205327 3.86928695244899 5 0.253384515785908 0.904747157403704 4.94186832681039 6 0.139846290897889 0.456681488763738 5.50347222033837 7 0.148082371938012 0.223097926838407 5.72881970122358 8 0.193976398804964 0.115041118545248 5.79098248264979 9 0.258076652158558 6.43200487911403E-02 5.7776032990503 10 0.332855516159152 3.91808850328527E-02 5.72796359880799 20 0.170673966811781 0.864389003690028 5.06493702949818 40 9.12491209293587E-02 0.020299500264818 5.98845137880582 Example Program (2) Solve the following delay differential equations (using grid and dense output). y'1(t) = -y1(t)*y2(t - 1) + y2(t - 10) y'2(t) = y1(t)*y2(t - 1) - y2(t) y'3(t) = y2(t) - y2(t - 10) (y1(t) = 5, y2(t) = 0.1, y3(t) = 1 at t <= 0) (discontinuous at t = 1, 2, ..., 10, 20 and 40) Sub Ex_Retard_r_2() Const N = 3, Ngrid = 12 Const Mxst = 500, LRCont = Mxst * (5 * N + 2), LICont = N Dim T As Double, Y(N - 1) As Double Dim Tend As Double, Tout As Double, Tgrid(Ngrid - 1) As Double Dim RTol(0) As Double, ATol(0) As Double, ITol As Long, Info As Long Dim RCont(LRCont - 1) As Double, ICont(LICont - 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, Y0 As Double, Y1 As Double, Y2 As Double Dim Y2L1 As Double, Y2L10 As Double, I As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) Iout = 1 T = 0: Y(0) = 5: Y(1) = 0.1: Y(2) = 1 For I = 0 To 9 Tgrid(I) = I + 1 Next Tgrid(10) = 20 Tgrid(11) = 40 Tend = 40 Tout = 5 Info = 0 IRev = 0 Do Call Retard_r(N, T, Y(), Tend, RTol(), ATol(), Ngrid, Tgrid(), RCont(), ICont(), Info, TT, YY(), YYp(), Irtrn, IRev, Iout) If IRev = 1 Then Y2L1 = YlagR(1, TT - 1, RCont(), ICont()) Y2L10 = YlagR(1, TT - 10, RCont(), ICont()) YYp(0) = -YY(0) * Y2L1 + Y2L10 YYp(1) = YY(0) * Y2L1 - YY(1) YYp(2) = YY(1) - Y2L10 ElseIf IRev = 5 Then While T >= Tout Y0 = YlagR(0, Tout, RCont(), ICont()) Y1 = YlagR(1, Tout, RCont(), ICont()) Y2 = YlagR(2, Tout, RCont(), ICont()) Debug.Print Tout, Y0, Y1, Y2 Tout = Tout + 5 Wend End If Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Retard_r: Info =", Info End Sub Function YlagR(II As Long, T As Double, RCont() As Double, ICont() As Long) As Double Dim YY As Double, IRev As Long IRev = 0 Do Call Ylag_r(II, T, RCont(), ICont(), YY, IRev) If IRev = 1 Then YY = 0.1 Loop While IRev <> 0 YlagR = YY End Function Example Results 5 0.253384515785908 0.904747157403704 4.94186832681039 10 0.332855516159152 3.91808850328527E-02 5.72796359880799 15 4.40303382758445 0.163106426937532 1.53385974547801 20 0.170673966811781 0.864389003690028 5.06493702949818 25 0.214955222878298 1.64093322939858E-02 5.86863544482771 30 4.87247653187963 7.33384924459654E-02 1.1541849756744 35 0.422930643367461 1.3593087991098 4.31776055752274 40 9.12491209293587E-02 0.020299500264818 5.98845137880582