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

Dverka() Sub Dverka ( N As  Long, F As  LongPtr, T As  Double, Y() As  Double, Tout As  Double, Tend As  Double, Tol As  Double, Mode As  Long, Info As  Long, Optional Neval As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional MaxIter As  Long = 0, Optional ErrCntl As  Long = 0, Optional Cnt As  Long = 0, Optional Hinit As  Double = 0, Optional Hmax As  Double = 0, Optional Hmin As  Double = 0, Optional Scal As  Double = 0, Optional Efloor As  Double = 0  ) Initial value problem of ordinary differential equations (6(5)-th order Runge-Kutta-Verner method) Purpose This program integrates a system of first order ordinary differential equations of the form dy/dt = f(t, y), 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. This program is the rewritten version of the 6(5)-th order Runge-Kutta-Verner code DVERK (Reference (1)). The dense output feature by Enright, et al. (Reference (2)) has been embedded. 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, Yp() As Double) Compute dy/dt and store in Yp(). End Sub where N is the number of equations, and dy/dt (= f(t, y)) is the computed derivative at given T and Y(). [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. However, if the integration is in backward direction, Tend <= Tout < T 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 for both forward (Tend > T) and backward (Tend < T) direction. [in] Tol Error tolerance. (Tol > 0) This program 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 that depend on the error control strategy specified by ErrCntl parameter. The default weight is 1/max(1, abs(y(i))), which provides a mixture of absolute and relative error control. [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,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] 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] ErrCntl (Optional) Error control indicator, which specifies the weights of the error estimate vector for error control. (default = 0) = 0: weight = 1/max(1, abs(y(i))) (mixture of absolute and relative error control) = 1: weight = 1 (absolute error control) = 2: weight = 1/abs(y(i)) (relative error control) = 3: weight = 1/max(abs(work[13]), abs(y(i))) (relative error control, unless abs(y(i)) is less than the floor value abs(Efloor)) [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) Parameter to determine maximum step size. Maximum step size = min(abs(Hmax), 2/abs(Scal)) (abs(Hmax) if Scal = 0). (default: Maximum step size = 2/abs(Scal) (2 if scale = 0)) [in] Hmin (Optional) Minimum step size. (default = 10*max(dwarf, rreb*max((weighted norm of y)/Tol, abs(T))), where dwarf is a very small positive machine number (= 1.0e-50) and rreb is the relative roundoff error bound.) [in] Scal (Optional) Measure of the scale of the problem. (default = 1) This parameter is used to determine the maximum step size and to modify the acceptance criterion. Larger values of Scal tend to make the method more reliable, [in] Efloor (Optional) Floor value used when ErrCntl = 3. Reference (1) netlib/ode (2) W. H. Enright, et al.: "Interpolants for Runge-Kutta Formulas", ACM Transactions on Mathematocal Software, Vol.12, No.3. (1986) Example Program 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_Dverka() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim Tol As Double, Mode As Long, Neval As Long, Info As Long Tol = 0.0000000001 '1.0e-10 Mode = 2 T = 0: Y(0) = 1: Y(1) = 2 Tend = 10 Info = 0 Do Tout = T + 1 Call Dverka(N, AddressOf F1, T, Y(), Tout, Tend, Tol, Mode, Info, Neval) Debug.Print T, Y(0), Y(1) 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 1 0.367879441171258 0.908181747039973 2 0.135335283236816 -0.28081155331094 3 4.97870683680297E-02 -0.940205428232924 4 1.83156388888794E-02 -0.635327981975176 5 6.73794699886029E-03 0.290400132462794 6 2.47875217650127E-03 0.962649038827388 7 9.11881965405588E-04 0.754814136309173 8 3.35462627940427E-04 -0.145164571180805 9 1.23409804260466E-04 -0.911006852080962 10 4.53999299169156E-05 -0.839026129147014 2539 0