Contex Contex_r Contr5 Contr5_r Contra Contra_r Contro Contro_r Contrp Contrp_r Dassl Dassl_r Debdf Debdf_r Radau Radau5 Radau5_r Radau_r Radaup Radaup_r Rodas Rodas_r Seulex Seulex_r

Debdf() Sub Debdf ( N As  Long, F As  LongPtr, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, Optional Djac As  LongPtr = NullPtr, Optional Ml As  Long = -1, Optional Mu As  Long = -1, Optional Mode As  Long = -1, Optional ITstop As  Long = -1, Optional Tstop As  Double  ) Initial value problem of ordinary differential equations (1~5-th order backward differentiation formula (BDF)) 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. Debdf is a variable order (one through five) backward differentiation formula code. It is the most complicated of the three choices. Debdf is primarily designed to solve stiff differential equations at crude to moderate tolerances. If the problem is very stiff at all, Derkf and Deabm will be quite inefficient compared to Debdf. However, Debdf will be inefficient compared to Derkf and Deabm on non-stiff problems because it uses much more storage, has a much larger overhead, and the low order formulas will not give high accuracies efficiently. Debdf is a driver for a modification of the code LSODE written by A. C. Hindmarsh. 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) Yp(i) = computed derivative at given T and Y() (i = 0 to N-1) End Sub where N is the number of equations, and Yp() is the computed derivatives at given T and Y(), i.e. Yp(i) = dyi/dt = fi(T, Y(0), ..., Y(N-1)) (i = 0 to N-1). The other variables than Yp() should not be altered. [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. [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. [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) = -3: The argument T had an illegal value. (T = Tout or T != previous Tout on continuation call) = -4: The argument Y() is invalid. = -5: The argument Tout had an illegal value. (Direction of integration is changed) = -6: The argument RTol() is invalid or had an illegal value. (RTol(i) < 0) = -7: The argument ATol() is invalid or had an illegal value. (ATol(i) < 0) = -8: The argument Info had an illegal value. (Info <> 0, 1, 2 or 11 to 15) = 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: Repeated convergence test failures on the last attempted step. User may reenter without changing Info to continue, but an inaccurate Jacobian may be the problem. = 15: Repeated error test failures on the last attempted step. User may reenter without changing Info to continue, but a singularity in the solution may be present. = 16: Infinite loop has been detected. [in] Djac (Optional) The user supplied subroutine, which computes the Jacobian (dfi/dyj) analytically, defined as follows. (default = NullPtr) Sub Djac(N As Long, T As Double, Y() As Double, Ypd() As Double) Ypd(i, j) = dfi/dyj (for i = 0 to N-1, j = 0 to N-1) End Sub Ypd() is 2-dimensional full matrix if Ml = N. It is in band matrix form if 0 <= Ml < N. The calculated value of Jacobian for given T and Y() should be set to Ypd(). The other variables than Ypd() should not be altered. If Djac is not provided (if Djac = NullPtr), Jacobian will be computed by numerical differentiation. [in] Ml (Optional) If Ml = N, Jacobian is handled as N x N full matrix. If 0 <= Ml < N, Jacobian is handled in band matrix form (advantageous if 2*Ml + Mu < N). In that case, Ml is the lower bandwidth. (Ml >= 0, Ml <= N) (default = N) [in] Mu (Optional) In the case of banded Jacobian, Mu is the upper bandwidth. If Ml = N, Mu is not referred. (Mu >= 0, Mu <= N) (default = 0) [in] Mode (Optional) Mode of operation. (default = 0) = 0: Return only at Tout (interval mode) = 1: Return at every step (intermediate output mode) (For other values, the default value will be used.) [in] ITstop (Optional) This routine may integrate past tout and interpolate to obtain the result at tout. However, it may not be permissible to integrate past some specific point because a discontinuity occurs there, or the function or its derivative is not defined beyond that point. In that case, the stopping point Tstop can be specified so that the routine will not integrate beyond that point. (default = 0) = 0: Not define the stopping point = 1: Define the stopping point at Tstop (for other values, the default value will be used) [in] Tstop (Optional if ITstop = 0) Stopping point value. Reference SLATEC (DEPAC) Example Program Solve the following initial value problem of ordinary differential equations (stiff problem). dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 1998*y1 - 1999*y2 + 1999*cos(t) - sin(t) (y1 = 1, y2 = 2 at t = 0) Sub F2(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(0) = -2 * Y(0) + Y(1) - Cos(T) Yp(1) = 1998 * Y(0) - 1999 * Y(1) + 1999 * Cos(T) - Sin(T) End Sub Sub Ex_Debdf() 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 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 Call Debdf(N, AddressOf F2, T, Y(), Tout, RTol(), ATol(), Info) If Info <> 1 Then Debug.Print "Error in Debdf: Info =", Info Exit Sub End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub Example Results 1 0.367879442394948 0.908181748259189 2 0.135335284386392 -0.280811552155648 3 0.049787069049436 -0.940205427530908 4 1.83156392354596E-02 -0.635327981611744 5 6.73794718964386E-03 0.290400132587821 6 2.47875225773856E-03 0.96264903888483 7 9.11881999298641E-04 0.75481413632096 8 3.35462642293389E-04 -0.145164571183809 9 1.2340981022923E-04 -0.911006852067105 10 4.53999332058006E-05 -0.839026129132657