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

Dverk_r() Sub Dverk_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, Tol As  Double, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, IRev As  Long, Optional Neval As  Long, Optional Naccept As  Long, Optional MaxEval As  Long = 0, Optional Int1 As  Long = 0, Optional Int2 As  Long = 0, Optional Hmin As  Double = 0, Optional Hmax As  Double = 0, Optional Scal As  Double = 0, Optional Hstart As  Double = 0, Optional Weight As  Long = 0, Optional Floor As  Double, Optional Cont As  LongPtr  ) Initial value problem of ordinary differential equations (6(5)-th order Runge-Kutta-Verner method) (reverse communication version) NOTE - THIS PROGRAM IS DEPRECATED AND WILL BE REMOVED IN THE NEXT VERSION. Purpose This routine integrates a system of first order ordinary differential equations of the form dy/dt = f(t, y), 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. Dverk is a 6(5)-th order Runge-Kutta-Verner code. It is efficient for non-stiff systems when derivative evaluations are not expensive. The dense output based on the reference paper below is supported (see Int2 and Example Program (2)). Dverk_r is the reverse communication version of Dverk. 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. [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] Tol Error tolerance. (Tol > 0) This routine attempts to control a norm of the local error in such a way that the global error is proportional to Tol. The norm is a max norm with weights. The default weight for the k-th component is 1/max(1, abs(y(k))), which therefore provides a mixture of absolute and relative error control. [in,out] Info [in] Control variable. Set info = 0 on the initial call. = 0: The routine will be initialized. = 1: Reentry after a normal return with new value of Tout. = 2, 3, 4: Resume from the interruption 1 or 2. [out] Return code. By examining this code, user can call this routine again when info = 1 to 4 as a next action if necessary. = -1: The argument N had an illegal value. (N < 1) = -3: The argument T had an illegal value. (Tout reached, or T changed on reentry) = -6: The argument Tol had an illegal value. (Tol <= 0) = -9: The argument lc had an illegal value. (lc too small) = -11: The argument lwork had an illegal value. (lwork too small) = -12: The argument Info had an illegal value. (Info < 0 or Info > 4) = 1: Successful exit (T = Tout). User can reenter with new value of Tout without changing any other variables. = 2: Interruption 1 (see Int1). = 3: Interruption 2 (see Int2, trial step is accepted). = 4: Interruption 3 (see Int2, trial step is rejected). = 11: The allowed maximum number of function evaluations has been exceeded. = 12: The value of minimum step size is greater than maximum step size, which probably means that the requested Tol is too small. = 13: Unable to satisfy the error requirement with a particular step size that is less than or equal to minimum step size, which may mean that Tol is too 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] 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 derivative values at TT and YY() in YYp(). Do not alter any variables other than YYp(). [out] Neval (Optional) Number of function evaluations. [out] Naccept (Optional) Number of accepted steps. [in] MaxEval (Optional) Maximum number of function evaluations. (default = 5000) (If MaxEval < 0, the absolute value will be used) (If MaxEval = 0, the default value will be used) [in] Int1 (Optional) Interrupt mode 1. (default = 0) = 0: Return only at Tout. (Normal mode) = 1: Interrupt the computation before taking a trial step. (Interrupt mode 1) (For other values, the default value will be used) [in] Int2 (Optional) Interrupt mode 2. (default = 0) = 0: Return only at Tout. (Normal mode) = 1: Interrupt the computation immediately after deciding whether or not to accept the trial step, with Info = 3 if it plans to accept, or Info = 4 if it plans to reject. (Interrupt mode 2) = 2: Interrupt mode 2 with dense output enabled. That is, DverkInt routine can be used to compute the interpolated values within latest step interval if returned with Info = 3 (Info = 1 for last step). One additional function evaluation per accepted step is required. (For other values, the default value will be used) [in] Hmin (Optional) Minimum step size. (default = 10*max(1.0e-50, eps*max((weighted norm of y)/tol, abs(x))) where eps is the machine epsilon) (If Hmin = 0, the default value will be used) [in] Hmax (Optional) Maximal step size is decided from Hmax and Scal as follows. (default = 0) Hmax = 0 and Scal = 0: 2 Hmax <> 0 and Scal <> 0: min(abs(Hmax), 2/abs(Scal)) Hmax <> 0 and Scal = 0: abs(Hmax) Hmax = 0 and Scal <> 0: 2/abs(Scal) [in] Scal (Optional) Reliability parameter (scaling factor of the problem). (default = 1) To be used for scaling and to decide the maximum step size. Larger values of Scal make results more reliable but require more number of function evaluations. (If Scal = 0, the default value will be used for scaling) [in] Hstart (Optional) Initial value of step size. (default = (max. step size)*Tol^(1/6)) (If Hstart < 0, the absolute value will be used) (If HHstart = 0, the default value will be used) [in] Weight (Optional) Maximum norm of weighted error estimate vector is used for error control. Its weight is specified as follows (default = 0) = 0: weight = 1/max(1, abs(Y(k))) (relative error control with floor value = 1) = 1: weight = 1 (absolute error control) = 2: weight = 1/abs(Y(k)) (relative control) = 3: weight = 1/max(abs(Floor), abs(Y(k))) (relative error control with specified floor value) (For other values, the default value will be used) [in] Floor (Optional if Weight <> 3) Floor value in the case of Weight = 3. [out] Cont (Optional) Dense output information to be used by DverkInt. Returned only when Info = 3 (Info = 1 for last step) and Int2 = 2. Reference (1) netlib/ode (2) W H Enright et al. "Interpolants for Runge-Kutta Formulas" ACM Transactions on Mathematical Software Vol.12, No.3, 1986, pp.193-218 Example Program (1) Solve the following initial value problem of ordinary differential equations. dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t) (y1 = 1, y2 = 2 at t = 0) Sub F1(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(0) = -2 * Y(0) + Y(1) - Cos(T) Yp(1) = 2 * Y(0) - 3 * Y(1) + 3 * Cos(T) - Sin(T) End Sub Sub Ex_Dverk_r() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim Tol As Double, Info As Long Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, IRev As Long Tol = 0.0000000001 '1.0e-10 T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Info = 0 Do Tout = T + 1 IRev = 0 Do Call Dverk_r(N, T, Y(), Tout, Tol, Info, TT, YY(), YYp(), IRev) If IRev = 1 Then Call F1(N, TT, YY(), YYp()) Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Dverk_r: Info =", Info Exit Sub End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub Example Results 1 0.367879441171258 0.908181747039973 2 0.135335283236816 -0.28081155331094 3 4.97870683680296E-02 -0.940205428232924 4 1.83156388888794E-02 -0.635327981975177 5 6.73794699886025E-03 0.290400132462795 6 2.47875217650124E-03 0.962649038827387 7 9.11881965405485E-04 0.754814136309172 8 3.35462627940704E-04 -0.145164571180804 9 1.23409804260653E-04 -0.911006852080961 10 4.53999299169892E-05 -0.839026129147014 Example Program (2) Solve the following initial value problem of ordinary differential equations (using dense output). dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t) (y1 = 1, y2 = 2 at t = 0) Sub F1(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(0) = -2 * Y(0) + Y(1) - Cos(T) Yp(1) = 2 * Y(0) - 3 * Y(1) + 3 * Cos(T) - Sin(T) End Sub Sub Ex_Dverk_r_2() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim Tol As Double, Info As Long Dim Y1(N - 1) As Double, Int2 As Long, Cont As LongPtr Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, IRev As Long Tol = 0.0000000001 '1.0e-10 T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Int2 = 2 Tout = T + 1 Info = 0 Do IRev = 0 Do Call Dverk_r(N, T, Y(), Tend, Tol, Info, TT, YY(), YYp(), IRev, Int2:=Int2, Cont:=Cont) If IRev = 1 Then Call F1(N, TT, YY(), YYp()) Loop While IRev <> 0 While (Info = 1 Or Info = 3) And T >= Tout Call DverkInt(N, Tout, Y1(), Cont) Debug.Print Tout, Y1(0), Y1(1) Tout = Tout + 1 Wend Loop While Tout <= Tend End Sub Example Results 1 0.367879441761895 0.9081817458561 2 0.135335283109193 -0.280811553054101 3 4.97870683242647E-02 -0.940205428144945 4 0.018315638476038 -0.635327981148217 5 6.73794741978772E-03 0.290400131618177 6 2.4787522740674E-03 0.962649038631452 7 9.11882292284201E-04 0.754814135654083 8 3.35462206406771E-04 -0.14516457033392 9 1.23409537546899E-04 -0.911006851545974 10 4.53999299134007E-05 -0.839026129147008