Deabm Deabm_r Derkfa Derkfa_r Dop853a Dop853a_r Dopn1210 Dopn1210_r Dopn43 Dopn43_r Dopn64 Dopn64_r Dopn86 Dopn86_r Dopri5a Dopri5a_r Dverka Dverka_r Odex2a Odex2a_r Odexa Odexa_r Retarda Retarda_r Ylaga

Retarda_r() Sub Retarda_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, Tend As  Double, RTol() As  Double, ATol() As  Double, Mode As  Long, Grid() As  Double, Cont() As  Double, ICont() As  Long, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, IRev As  Long, Optional Neval As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional MaxIter As  Long = 0, Optional Nstiff As  Long = 0, Optional Ngrid As  Long = 0, Optional Mxst As  Long = 0, Optional Cnt As  Long = 0, Optional Hinit As  Double = 0, Optional Hmax As  Double = 0, Optional Fac1 As  Double = 0, Optional Fac2 As  Double = 0, Optional Safe As  Double = 0, Optional Beta As  Double = 0  ) Initial value problem of delay differential equations (5(4)-th order Dorman-Prince method) (Reverse communication version) Purpose This program 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 y is the n vector and the system consists of n equations. t0 and y0 are the given initial values of t and y, respectively. f depends on the values at previous points y(t - τ), y(t - τ2), ... as well as y(t). This program is the rewritten version of the 5(4)-th order Dormand-Prince code for delay equations, RETARD (Reference (1)). Parameters [in] N Number of differential equations. (N >= 1) [in,out] T Independent variable t. This program integrates from the initial value of t to Tend. Depending on the setting of Mode, this program may return at Tout or every successful step to provide an intermediate result. [in] An initial value of t (t0). [out] T = Tend if the integration has been completed. Or T = Tout or T of the latest step depending on Mode in the case of intermediate return. [in,out] Y() Array Y(LY - 1) (LY >= N) Dependent variable y. [in] Initial value of y (y0) at initial t (t0). [out] The value of y (computed numerical solution) at t = T. [in] Tout Set Tout to the point (t) at which an intermediate result is desired in mode = 2 or 3. Tout is not referenced in Mode = 0 or 1. In Mode 2 or 3, y at Tout is computed and returned with Info = 2 or 3 respectively. To continue the integration to get result at new Tout, it is possible to call this program again with new Tout and without changing other variables including Info. T < Tout <= Tend is required. If Tend is reached before Tout, the integration is terminated immediately and T = Tend and Info = 0 is returned. [in] Tend The point at which an integration is completed. If Tend is reached, the integration is terminated and T = Tend and Info = 0 is returned. The integration is possible only for forward (Tend > T) direction. [in] Rtol() Array Rtol(LRtol - 1) (LRtol >= 1) (all elements of Rtol() >= 0) The relative error tolerance(s) to specify how accurately compute the solution. This parameter may be a scalar or an array according to the values of LRtol and LAtol (Scalar if LRtol < N or LAtol < N, array if LRtol >= N and LAtol >= N). Rtol is used together with Atol in a local error test at each step. The integration is performed using step sizes which are automatically selected so as to achieve the following criteria. Scalar case (i = 0 to N - 1): (local error of Y(i)) <= Rtol(0)*Abs(Y(i)) + Atol(0) Both Rtol(0) and Atol(0) must not be 0 at the same time. Array case (i = 0 to N - 1): (local error of Y(i)) <= Rtol(i)*Abs(Y(i)) + Atol(i) Both Rtol(i) and Atol(i) must not be 0 at the same time. [in] ATol() Array Atol(LAtol - 1) (LAtol >= 1) (all elements of Atol() >= 0) The absolute error tolerance(s) to specify how accurately compute the solution. This parameter may be a scalar or an array according to the values of LRtol and LAtol (Scalar if LRtol < N or LAtol < N, array if LRtol >= N and LAtol >= N). Atol is used together with Rtol in a local error test at each step (Refer to Rtol above). [in] Mode Mode of operation. This program integrates from t0 to Tend. Four modes of operation are provided depending on the timing of returning imtermediate results. When Tend is reached, the integration is terminated and Info = 0 is returned in any of four modes. = 0: Not return until Tend. Tout is not referenced. = 1: Returns at every successful step (Info = 1 is returned). Tout is not referenced. The end point of the last step is returned in T. In the final step, the step size is adjusted to fit into Tend. = 2: Returns at Tout during the ingegration to provide an intermediate result (T = Tout and Info = 2 are returned). To continue the integration to get result at new Tout, it is possible to call this program again with new Tout. In the final step, the step size is adjusted to fit into Tout. = 3: Returns at Tout during the ingegration to provide an intermediate result (T = Tout and Info = 3 are returned). To continue the integration to get result at new Tout, it is possible to call this program again with new Tout. Different from Mode = 2, the value of Y() at Tout is computed by interpolation, and the step size is not adjusted even in the last step before Tout. The steps through integration are same with those of Mode = 1. [in] Grid() Array Grid(LGrid - 1) (LGrid >= Ngrid) If the solution has the discontinuities in its derivatives within the integration interval, by specifying those points (grid points), the precision and the computation time can be improved. If Ngrid > 0, the grid points must be stored in Grid(0), ..., Grid(Ngrid - 1) in ascending order (except starting and end points). If Ngrid = 0, Grid() is not referenced. [in,out] Cont() Array Cont(LCont - 1) Control information area. The contents of this array must not be altered. This array is also used to share information with Ylaga. When this routine is called with Info = -1, the necessary LCont value is computed and returned in Info without any other computation. [in,out] ICont() Array Icont(LICont - 1) (LICont >= 50) Integer control information area. The contents of this array must not be altered. This array is also used to share information with Ylaga. [in,out] Info [in] Control code. = 0: Set Info = 0 on the initial call to start new problem. All variables will be initialized and the computation will begin. = 1, 2, 3: When returned with Info = 1, 2 or 3, it is possible to call this program again with new Tout and without changing Info to continue the integration. [out] Return code. = 0: Successful exit. Integration to Tend completed. < 0: The (-Info)-th argument is invalid. = 1: Returned to provide intermediate result in Mode = 1. It is possible to call this program again to forward to the next step. = 2: Returned at T = Tout to provide intermediate result in Mode = 2. It is possible to call this program again with new Tout. = 3: Returned at T = Tout to provide intermediate result in Mode = 3. It is possible to call this program again with new Tout. = 11: (Error) Maximum number of steps exceeded. = 12: (Error) Step size becomes too small. = 13: (Error) Problem is probably stiff. [out] TT IRev = 1: Refer to YYp(). [out] YY() Array YY(LYY - 1) (LYY >= N) IRev = 1: Refer to YYp(). [in] YYp() Array YYp(LYYp - 1) (LYYp >= N) IRev = 1: Compute derivatives at given T = TT and Y = YY(), dy/dt = f(t, y(t), y(t - τ), y(t - τ2), ...), and store in YYp() in the next call. y(t - τ), y(t - τ2), ... can be computed by using Ylaga. [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 program again without changing IRev. = 0: Computation finished. Check return code in Info. = 1: Computed derivatives at T = TT and Y = YY() should be stored in YYp(). Do not alter any variables other than YYp(). [out] Neval (Optional) Number of function evaluations. [out] Nstep (Optional) Number of all computed steps. [out] Naccept (Optional) Number of accepted steps. [out] Nreject (Optional) Number of rejected steps. [in] Maxiter (Optional) Maximum number of allowed steps. (default = 10000) [in] Nstiff (Optional) Test for stiffness is activated after every Nstiff steps. (default = 1000)\ If Nstiff < 0, the stiffness test is not activated. [in] Ngrid (Optional) Number of grid points which are stored in Grid(). (default = 0) [in] Mxst (Optional) Specifies for how many previous steps the information is retained. Ylaga uses this information to compute a value y(t - τ) in previous step. (default = 1000) [in] Cnt (Optional) Specifies when counters (Neval, Nstep, Naccept and Nreject) are reset to zero. (default = 0) = 0: Reset if Info = 0. = 1: Do not reset even if Info = 0. [in] Hinit (Optional) Initial step size. (default = estimated by the program) [in] Hmax (Optional) Maximaum step size. (default = Abs(Tend - T)) [in] Fac1,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. [in] Safe (Optional) The safety factor in step size prediction. (0.0001 < Safe < 1) (default = 0.9) [in] Beta (Optional) The parameter for the stabilized step size control (Beta <= 0.2). (default = 0.04 (0 if Beta < 0 is specified)) Reference (1) 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 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 Phi(I As Long, N As Long, T As Double, Y() As Double, Info As Long) If I = 1 Then Y(0) = 5 If I = 2 Then Y(0) = 0.1 If I = 3 Then Y(0) = 1 End Sub Sub Ex_Retarda_r() Const N = 3, Ngrid = 11 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, Grid(Ngrid) As Double Dim Mode As Long, Neval As Long, Info As Long, I As Long Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, IRev As Long Dim Cont() As Double, ICont() As Long, Y2L1(0) As Double, Y2L10(0) As Double RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) For I = 0 To 9 Grid(I) = I + 1 Next Grid(10) = 20 Mode = 2 Info = -1 Call Retarda_r(N, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Grid(), Cont(), ICont(), Info, TT, YY(), YYp(), IRev, Ngrid:=Ngrid) ReDim Cont(Info - 1), ICont(50 - 1) T = 0: Y(0) = 5: Y(1) = 0.1: Y(2) = 1 Tend = 40 Info = 0 Do Tout = T + 5 IRev = 0 Do Call Retarda_r(N, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Grid(), Cont(), ICont(), Info, TT, YY(), YYp(), IRev, Neval, Ngrid:=Ngrid) If IRev = 1 Then Call Ylaga(2, N, TT - 1, Y2L1(), AddressOf Phi, Cont(), ICont()) Call Ylaga(2, N, TT - 10, Y2L10(), AddressOf Phi, Cont(), ICont()) YYp(0) = -YY(0) * Y2L1(0) + Y2L10(0) YYp(1) = YY(0) * Y2L1(0) - YY(1) YYp(2) = YY(1) - Y2L10(0) End If Loop While IRev <> 0 Debug.Print T, Y(0), Y(1), Y(2) Loop While Info >= 1 And Info <= 3 Debug.Print Neval, Info End Sub Note - If Mode is chenged to 3 instead of 2 in the same program, the dense outout (interpolation) is used. Example Results 5 0.253384515787142 0.904747157404 4.94186832680886 10 0.332855516159117 3.91808850328078E-02 5.72796359880807 15 4.40303383018292 0.163106427890271 1.5338597419268 20 0.170673965282798 0.864389000478564 5.06493703423863 25 0.214955223950597 1.64093322434754E-02 5.86863544380592 30 4.87247654002283 7.33384930344983E-02 1.15418496694268 35 0.42293063684405 1.35930879446553 4.31776056869042 40 9.12491214440587E-02 2.02995001131839E-02 5.98845137844276 4438 0