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

Dassl() Sub Dassl ( N As  Long, Res As  LongPtr, T As  Double, Y() As  Double, Yp() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, Optional Jac As  LongPtr = NullPtr, Optional Ml As  Long = -1, Optional Mu As  Long, Optional Neval As  Long, Optional Njac As  Long, Optional Nstep As  Long, Optional Mode As  Long, Optional ITstop As  Long, Optional Tstop As  Double, Optional Hmax As  Double, Optional H0 As  Double, Optional Maxord As  Long, Optional NonNeg As  Long, Optional NoInit As  Long  ) Solution of differential algebraic equation (DAE) (1~5-th order backward differentiation formula (BDF)) NOTE - THIS PROGRAM IS DEPRECATED AND WILL BE REMOVED IN THE NEXT VERSION. Purpose This routine solves a implicit system of differential algebraic equations (DAEs) of the form f(t, y, y') = 0, y = y0 and y' = y'0 at t = t0 where t0, y0 and y'0 are the given initial values of t, y and y', respectively. y and y' may be vectors if the above is a system of differential algebraic equations. In this routine, the derivatives y' are approximated by backward differentiation formulae (BDF) of orders 1 through 5, and the resulting nonlinear system at each time-step is solved by Newton's method. Parameters [in] N Number of differential equations. (N >= 1) [in] Res The user supplied subroutine, which calculates the residual of the differential algebraic system delta = f(t, y, y'), defined as follows. Sub Res(N As Long, T As Double, Y() As Double, Yp() As Double, Delta() As Double, IRes As Long) Delta(i) = residual at T, Y() and Yp() (i = 0 〜 N-1) End Sub where N is the number of equations, Y() is the function values at T and Yp() is the derivatives at T. IRes is an integer flag which is always equal to 0 on input. Subroutine Res should alter IRes only if it encounters an illegal value of Y() or a stop condition. Set IRes = -1 if an input value is illegal, and Dassl will try to solve the problem without getting IRes = -1. If IRes = -2, Dassl will return control to the calling program with Info = 11. Do not alter variables N, T, Y() and Yp() in Res. [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 at initial T. [out] Computed solution approximation at last T (= Tout in interval mode). [in,out] Yp() Array Yp(LYyp - 1) (LYp >= N) [in] Initial values of the derivatives of independent variables at initial T. [out] Computed derivatives of the 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 13: When returned from the routine with Info = 1, 2 or 11 to 13, 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 13 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 Yp() is invalid. = -6: The argument Tout had an illegal value. (Direction of integration is changed) = -7: The argument RTol() is invalid or had an illegal value. (RTol(i) < 0) = -8: The argument ATol() is invalid or had an illegal value. (ATol(i) < 0) = -9: The argument Info had an illegal value. (Info <> 0, 1, 2 nor 11 to 13) = -18: The argument Tstop had an illegal value. (Tstop is behind T or Tout) = 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 automatically relaxed tolerances. User can also reenter with manually changed 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 error test failure = 15: Repeated corrector convergence test failure = 16: The matrix of partial derivatives is singular = 17: Multiple convergence test failure in this step = 18: The corrector could not converge because IRes = -1 = 19: IRes = -2 was encountered = 20: Failed to compute the initial yprime = 21: Infinite loop has been detected. [in] Jac (Optional) The user supplied subroutine, which computes the matrix of partial derivatives, defined as follows. (default = NullPtr) Sub Jac(N As Long, T As Double, Y() As Double, Yp() As Double, Ypd() As Double, Cj As Double) Store dfi/dyj + Cj*dfi/dyj' in Ypd() (i = 0 to N-1, j = 0 to N-1). End Sub The matrix of partial derivatives at T, Y() and Tp() should be stored in Ypd() in N x N general matrix form if Ml = N, or 2Ml+Mu+1 x N band matrix form if 0 <= Ml < N. See below for further details of band matrix form. Since all elements of Ypd() are initialized to 0 before calling Jac, only the non-zero elements can be stored. Do not alter variables N, T, Y(), Yp() and Cj in Jac. If Jac is not provided (if Jac = NullPtr), the matrix of partial derivatives will be computed by finite difference approximation. [in] Ml (Optional) If Ml = N, the matrix of partial derivatives is handled as N x N full matrix. If 0 <= Ml < N, the matrix 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 the band matrix, Mu is the upper bandwidth. If Ml = N, Mu is not referred. (Mu >= 0, Mu <= N) (default = 0) [out] Neval (Optional) Number of function evaluations. (Those for partial derivatives evaluations are not included) [out] Njac (Optional) Number of evaluations of the matrix of partial derivatives. (Those by finite differences are included) [out] Nstep (Optional) Number of computed steps. [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. [in] Hmax (Optional) Maximum step size. (Hmax > 0) (default = decided by the program) (If Hmax <= 0, the default value will be used) [in] H0 (Optional) Initial step size. (H0 > 0) (default = decided by the program) (If H0 <= 0, the default value will be used) [in] Maxord (Optional) Maximum order. (1 <= Maxord <= 5) (default = 5) (Maxord <= 0 or Maxord >= 6, the default value will be used) [in] NonNeg (Optional) Specify if the solutions to the equations will always be nonnegative. (default = 0) = 0: Solve the problem without invoking any special nonnegativity constraints. = 1: Solve the problem with invoking special nonnegativity constraints. (For other values, the default value will be used.) [in] NoInit (Optional) Specify if the initial t, y, and y' are consistent, i.e. f(t, y, y') = 0 at the initial time. (default = 0) = 0: The initial t, y, and y' are consistent. = 1: The precise initial derivatives are unknown. Request to try to compute y' by the program. Set Yp() to the initial approximation values. If no idea, set Yp() to 0. (For other values, the default value will be used.) Further Details The special band matrix form with 2*Ml+Mu+1 rows x N columns used in this routine is illustrated by the following example, when N = 6, Ml = 2 and Mu = 1. General matrix form: a11 a12 0 0 0 0 a21 a22 a23 0 0 0 a31 a32 a33 a34 0 0 0 a42 a43 a44 a45 0 0 0 a53 a54 a55 a56 0 0 0 a64 a65 a66 Band matrix form: * * * + + + * * + + + + * a12 a23 a34 a45 a56 a11 a22 a33 a44 a55 a66 a21 a32 a43 a54 a65 * a31 a42 a53 a64 * * The first ml rows marked + are used as the work area. The elements marked * are not used by the routine. Reference SLATEC (DASSL) Example Program Solve the following differential algebraic equations (Robertson problem). -0.04*y1 + 10000*y2*y3 - y1' = 0 0.04*y1 - 10000*y2*y3 - 30000000*y2^2 - y2' = 0 y1 + y2 + y3 -1 = 0 (y1 = 1, y2 = 0, y3 = 0, y1' = y2' = y3' = 0 at t = 0) Sub R2(N As Long, T As Double, Y() As Double, Yp() As Double, Delta() As Double, Ires As Long) Delta(0) = -0.04 * Y(0) + 10000 * Y(1) * Y(2) - Yp(0) Delta(1) = 0.04 * Y(0) - 10000 * Y(1) * Y(2) - 30000000 * Y(1) ^ 2 - Yp(1) Delta(2) = Y(0) + Y(1) + Y(2) - 1 End Sub Sub Ex_Dassl() Const N = 3 Dim T As Double, Y(N - 1) As Double, Yp(N - 1) As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, ITol As Long Dim Tend As Double, Info As Long '--- Parameters RTol(0) = 0.0000001: ATol(0) = 0.0000001: ITol = 0 '--- Initial values T = 0 Y(0) = 1: Y(1) = 0: Y(2) = 0 Yp(0) = 0: Yp(1) = 0: Yp(2) = 0 '--- Start calculation Tend = 10 Info = 0 Do Tout = T + 1 Call Dassl(N, AddressOf R2, T, Y(), Yp(), Tout, RTol(), ATol(), Info) If T >= Tout Then Debug.Print T, Y(0), Y(1), Y(2) End If If T >= Tend Then Exit Do Loop While Info = 1 If Info <> 1 Then Debug.Print "Error in Dassl: Info =", Info End Sub Example Results 1 0.966459990474207 3.07463058949132E-05 3.35092632198979E-02 2 0.941609675690738 2.7017862636127E-05 5.83633064466271E-02 3 0.921884807943058 2.43833786001128E-05 7.80908086783409E-02 4 0.905518993577723 2.24047935986599E-05 9.44586016286797E-02 5 0.891518109563595 2.08527025402239E-05 0.108461037733863 6 0.87927028562438 1.95949526270386E-05 0.120710119422982 7 0.868373377809956 1.85498046053848E-05 0.13160807238542 8 0.85854920043329 1.76638124059415E-05 0.141433135754299 9 0.849597711635835 1.69004734510369E-05 0.150385387890728 10 0.841370267657359 1.62339372132848E-05 0.158613498405457