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

Derkf_r() Sub Derkf_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, IRev As  Long, Optional Mode As  Long = -1, Optional Cont As  LongPtr  ) Initial value problem of a system of first order ordinary differential equations (5(4)-th order Runge-Kutta-Fehlberg 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. This is the code in the package of differential equation solvers DEPAC, consisting of the codes Derkf, Deabm, and Debdf. Derkf is a 5(4)-th order Runge-Kutta-Fehlberg code. It is the simplest of the three choices, both algorithmically and in the use of the code. It is primarily designed to solve non-stiff and mildly stiff differential equations when derivative evaluations are not expensive. It should generally not be used to get high accuracy results nor answers at a great many specific points. Because Derkf has very low overhead costs, it will usually result in the least expensive integration when solving problems requiring a modest amount of accuracy and having equations that are not costly to evaluate. Derkf is a driver for a modification of the code RKF45 written by H. A. Watts and L. F. Shampine. The dense output based on the reference paper below is supported (see Mode and Example Program (2)). Derkf_r is the reverse communication version of Derkf. 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. It is possible to continue the integration to get result at new Tout. This is the interval mode of operation. T should be equal to the previous Tout on continuation call. It is also possible for the routine to return with the solution at each intermediate step on the way to Tout. This is the intermediate-output mode of operation. This mode is a good way to proceed if you want to see the behavior of the solution. The mode of operation is specified by the parameter Mode. [in] Initial value of the independent variable T. [out] Last value of the independent variable T of the final step (normally equals to Tout in interval mode). The solution was successfully advanced to this point. [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 (= Tout in interval mode). [in] Tout Set Tout to the point at which a solution is desired. You can take Tout = T, in which case the code will evaluate the derivative of the solution at T and return. Integration either forward in T (Tout > T) or backward in T (Tout < T) is permitted. It is, however, not allowed to change the direction of integration without restarting. The code advances the solution from T to Tout using step sizes which are automatically selected so as to achieve the desired accuracy. If you wish, the code will return with the solution and its derivative following each intermediate step (intermediate-output mode) so that you can monitor them, but you still must provide Tout in accord with the basic aim of the code. [in] RTol() 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 an array (LRTol = N). The array can be used to control the error test more precisely. If LRTol = 2, ... or N-1, LRTol = 1 is assumed. If LRTol > N, LRTol = N is assumed. The tolerances are used by the code in a local error test at each step which requires roughly that   abs(local error) <= RTol(0)*abs(Y(i)) + ATol(0) (if LRTol or LATol = 1)     or   abs(local error) <= RTol(i)*abs(Y(i)) + ATol(i) (if LRTol and LATol >= N) for each component of Y() (i = 0 to N-1). Setting RTol(i) = 0 results in a pure absolute error test on that component. If you want relative accuracies smaller than about 1.0e-8, you should not ordinarily use Derkf. The code Deabm in DEPAC obtains stringent accuracies more efficiently. [in] ATol() 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 an array (LATol = N). The array can be used to control the error test more precisely. If LATol = 2, ... or N-1, LATol = 1 is assumed. If LATol > N, LATol = N is assumed. The tolerances are used by the code in a local error test at each step which requires roughly that   abs(local error) <= RTol(0)*abs(Y(i)) + ATol(0) (if LRTol or LATol = 1)     or   abs(local error) <= RTol(i)*abs(Y(i)) + ATol(i) (if LRTol and LATol >= N) for each component of Y() (i = 0 to N-1). Setting ATol(i) = 0 results in a pure relative error test on that component. If you want relative accuracies smaller than about 1.0e-8, you should not ordinarily use Derkf. The code Deabm in DEPAC obtains stringent accuracies more efficiently. [in,out] Info [in] Control code. = 0: Set Info = 0 on the initial call for the problem (to start new problem). The routine will be initialized and the computation for the new problem will be started. = 1, 2 or 11 to 15: When returned from the routine with Info = 1, 2 or 11 to 15, user can reenter without changing Info to continue computation. See descriptions below. [out] Return code. By examining this code, user can call this routine again when Info = 1, 2 or 11 to 15 as a next action if necessary. = -1: The argument N had an illegal value. (N < 1) = -2: The argument T had an illegal value. (T = Tout or T != previous Tout on continuation call) = -3: The argument Y() is invalid. = -4: The argument Tout had an illegal value. (Direction of integration is changed) = -5: The argument RTol() is invalid or had an illegal value. (RTol(i) < 0) = -6: The argument ATol() is invalid or had an illegal value. (ATol(i) < 0) = -7: The argument Info had an illegal value. (Info <> 0, 1, 2 or 11 to 15) = -9: The argument YY() is invalid. = -10: The argument YYp() is invalid. = 1: Successful exit (T = Tout). User can reenter with a new Tout, which must be different from T. = 2: Interruption in intermediate-output mode (still not reached Tout). User can reenter to resume for another step in the direction of Tout. = 11: Maximum number of steps (10000) exceeded. User can reenter to resume. Additional 10000 steps will be allowed. = 12: Error tolerances are too stringent. User can reenter to continue with relaxed tolerances. = 13: Pure relative error test (ATol = 0) is impossible because computed solution is zero. User can reenter with positive ATol value to continue computation. = 14: Problem is probably stiff. User can reenter to continue, but the other programs such as Debdf in DEPAC which handle stiff problems efficiently are recommended. = 15: Natural step size is being restricted by too frequent output. User can reenter to continue, but other programs with dense output feature are recommended. = 16: Infinite loop has been detected. [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(). [in] Mode (Optional) Mode of operation. (default = 0) = 0: Return only at Tout (interval mode) = 1: Return at every step (intermediate output mode) = 2: Return at every step (intermediate output mode). Dense output is enabled, i.e. DerkfInt routine can be used to compute the interpolated values within latest step interval when returned with Info = 2 (Info = 1 for last step). (For other values, the default value will be used.) [out] Cont (Optional) Dense output information to be used by DerkfInt. Returned only when Info = 2 (Info = 1 for last step) and Mode = 2. Reference (1) SLATEC (DEPAC) (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_Derkf_r() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, Info As Long Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, IRev As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Info = 0 Do Tout = T + 1 IRev = 0 Do Call Derkf_r(N, T, Y(), Tout, RTol(), ATol(), 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 Derkf_r: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub Example Results 1 0.367879441136182 0.908181747107717 2 0.135335283245885 -0.280811553333642 3 4.97870683958675E-02 -0.940205428291068 4 1.83156389160862E-02 -0.635327982031259 5 6.73794699054023E-03 0.290400132477777 6 2.47875214648626E-03 0.9626490388872 7 9.11881936224063E-04 0.754814136368097 8 3.35462623309683E-04 -0.145164571170354 9 1.23409831706352E-04 -0.91100685213574 10 4.5399959961099E-05 -0.839026129207695 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_Derkf_r_2() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim Y1(N - 1) As Double, Mode As Long, Cont As LongPtr Dim RTol(0) As Double, ATol(0) As Double, Info As Long Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, IRev As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Mode = 2 Tout = T + 1 Info = 0 Do IRev = 0 Do Call Derkf_r(N, T, Y(), Tend, RTol(), ATol(), Info, TT, YY(), YYp(), IRev, Mode, Cont) If IRev = 1 Then Call F1(N, TT, YY(), YYp()) Loop While IRev <> 0 If Info <> 1 And Info <> 2 Then Debug.Print "Error in Derkf_r: Info =", Info Exit Do End If While T >= Tout Call DerkfInt(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.367879441127794 0.908181747124133 2 0.135335283265252 -0.280811553374618 3 4.97870684046614E-02 -0.940205428309061 4 1.83156389289378E-02 -0.635327982057511 5 6.73794698644643E-03 0.290400132486541 6 2.4787521426282E-03 0.962649038895081 7 9.11881915710959E-04 0.754814136410194 8 3.35462644956918E-04 -0.145164571216893 9 1.23409854714012E-04 -0.911006852183254 10 4.53999601058471E-05 -0.839026129207973