Debdf Debdf_r Radaua Radaua_r Rodasa Rodasa_r Seulexa Seulexa_r

Radaua_r() Sub Radaua_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, Tend As  Double, RTol() As  Double, ATol() As  Double, Mode As  Long, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, YYpd() As  Double, IRev As  Long, Optional Neval As  Long, Optional Njac As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional Ndec As  Long, Optional Nsol As  Long, Optional Ijac As  Long, Optional Mljac As  Long = -1, Optional Mujac As  Long, Optional Imas As  Long, Optional Mlmas As  Long = -1, Optional Mumas As  Long, Optional Hes As  Long, Optional MaxIter As  Long, Optional Nit1 As  Long, Optional Startn As  Long, Optional Nind1 As  Long, Optional Nind2 As  Long, Optional Nind3 As  Long, Optional Pred As  Long, Optional M1 As  Long, Optional M2 As  Long, Optional Nsmax As  Long, Optional Nsmin As  Long, Optional Ns As  Long, Optional Cnt As  Long, Optional Hinit As  Double, Optional Hmax As  Double, Optional Thet As  Double, Optional Facl As  Double, Optional Facr As  Double, Optional Safe As  Double, Optional Quot1 As  Double, Optional Quot2 As  Double, Optional Vitu As  Double, Optional Vitd As  Double, Optional Hhou As  Double, Optional Hhod As  Double  ) Initial value problem of ordinary differential equations (variable (5, 9, 13-th) order implicit Runge-Kutta method (Radau IIA)) (Reverse communication version) Purpose This program computes a numerical solution of a stiff (or differential algebraic) system of first order ordinary differential equations of the form M * dy/dt = f(t, y), y = y0 at t = t0 where y is the n vector and the system consists of n equations. M is the mass matrix. t0 and y0 are the given initial values of t and y, respectively. This code is the rewritten version of the variable (5, 9, 13-th) order implicit Runge-Kutta method (Radau IIA)) code, RADAU (Reference (1)). This program is the reverse communication version of Radaua. Parameters [in] N Number of differential equations. (N >= 1) [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] Rtol() Array Rtol(LRtol - 1) (LRtol >= 1) (all elements of Rtol() >= 0) The relative error tolerance(s) to specify how accurately compute the solution. This parameter may be a scalar or an array according to the values of LRtol and LAtol (Scalar if LRtol < N or LAtol < N, array if LRtol >= N and LAtol >= N). Rtol is used together with Atol in a local error test at each step. The integration is performed using step sizes which are automatically selected so as to achieve the following criteria. Scalar case (i = 0 to N - 1): (local error of Y(i)) <= Rtol(0)*Abs(Y(i)) + Atol(0) Both Rtol(0) and Atol(0) must not be 0 at the same time. Array case (i = 0 to N - 1): (local error of Y(i)) <= Rtol(i)*Abs(Y(i)) + Atol(i) Both Rtol(i) and Atol(i) must not be 0 at the same time. [in] ATol() Array Atol(LAtol - 1) (LAtol >= 1) (all elements of Atol() >= 0) The absolute error tolerance(s) to specify how accurately compute the solution. This parameter may be a scalar or an array according to the values of LRtol and LAtol (Scalar if LRtol < N or LAtol < N, array if LRtol >= N and LAtol >= N). Atol is used together with Rtol in a local error test at each step (Refer to Rtol above). [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) Matrix is repeatedly singular. [out] TT IRev = 1: Refer to YYp(). IRev = 2: Refer to YYpd(). [out] YY() Array YY(LYY - 1) (LYY >= N) IRev = 1: Refer to YYp(). IRev = 2: Refer to YYpd(). [in] YYp() Array YYp(LYYp - 1) (LYYp >= N) IRev = 1: Compute derivatives at given T = TT and Y = YY(), dy/dt = f(t, y), and store in YYp() in the next call. [in] YYpd() Array YYpd(LYYpd1 - 1, LYYpd2 - 1) (LYYpd1 >= max(Ljac, Lmas), LYYpd2 >= N) IRev = 2: Compute Jacobian at given T = TT and Y = YY(), and store in YYpd() in the next call in full matrix form (general n x n two dimensional array) (Ljac = N) if Mljac = N, or in band matrix form (Ljac = Mljac + Mujac + 1) if Mljac < N. Ljac = 0 if Ijac = 0 IRev = 4: Mass matrix should be set in YYpd() in the next call in full matrix form if Mljac = N, or in band matrix form if Mljac < N. Lmas = 0 if Imas = 0. [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 program again without changing IRev. = 0: Computation finished. Check return code in Info. = 1: Computed derivatives at T = TT and Y = YY() should be stored in YYp(). Do not alter any variables other than YYp(). = 2: Computed Jacobian at T = TT and Y = YY() should be stored in YYpd(). Do not alter any variables other than YYpd(). = 4: Set the mass matrix in YYpd(). Do not alter any variables other than YYpd(). [out] Neval (Optional) Number of function evaluations. [out] Njac (Optional) Number of Jacobian evaluations (either analytically or numerically). [out] Nstep (Optional) Number of all computed steps. [out] Naccept (Optional) Number of accepted steps. [out] Nreject (Optional) Number of rejected steps. [out] Ndec (Optional) Number of LU decompositions. [out] Nsol (Optional) Number of forward-backward substitutions. [in] Ijac (Optional) Specify how to compute the Jacobian. = 0: Jacobian is computed by finite difference approximation. Never returns with IRev = 2. = 1: Jacobian is computed by IRev = 2 by user. [in] Mljac (Optional) The lower bandwidth of Jacobian. (0 <= Mljac <= N) (default = N) If Mljac = N, Jacobian is stored as full matrix. If Mljac < N, Jacobian is stored in band matrix form. [in] Mujac (Optional) The upper bandwidth of Jacobian. (0 <= Mujac <= N) (default = 0) If Mljac = N, Mujac is ignored. [in] Imas (Optional) Specify whether mass matrix M is identity matrix. = 0: M is the identity matrix. Never returns with IRev = 4. = 1: M is supplied by IRev = 4 by user. [in] Mlmas (Optional) The lower bandwidth of mass matrix M. (0 <= Mljac <= N) (default = N) If Mlmas = N, M is stored in full matrix form. If Mlmas < N, M is stored in band matrix form. [in] Mumas (Optional) The upper bandwidth of mass matrix M. (0 <= Mumas <= N) (default = 0) If Mlmas = N, Mumas is ignored. [in] Hes (Optional) Specifies whether to transform the Jacobian matrix to Hessenberg form. It is advantageous for large systems with full Jacobian. Not effective for banded Jacobian nor for the system with IMAS = 1. = 0: Transform Jacobian. = 1: Do not transform Jacobian. [in] Maxiter (Optional) Maximum number of allowed steps. (default = 10000) [in] Nit1 (default) Maximum number of Newton iterations for the solution of the implicit system in each step. (default = 7) [in] Startn (default) Specifies the starting value of Newton iteration. (default = 0) = 0: Use extrapolated solution as starting value. = 1: Use 0 as starting value. The latter is recommended if Newton method has difficulties with convergence. [in] Nind1,Nind2,Nind3 (default) Dimensions of the index 1, 2 and 3 variables. (Nind1 > 0, Nind1 + Nind2 + Nind3 = N) (default: Nind1 = N, Nind2 = 0, Nind3 = 0) These parameters are important for differential algebraic equations (DAEs) of index > 1. The function subroutine should be written such that the index 1, 2 and 3 variables appear in this order. In estimating the error, the index 2 variables are multiplied by h, the index 3 variables by h^2. For ordinary defferential equations (ODEs), Nind1 = N. [in] Pred (default) Specifies the step size strategy. (default = 1) = 1: Model predictive controller (Gustafsson). = 2: Classical step size control. [in] M1,M2 (default) If the differential system has the special structure that y(i)' = y(i + m2) (i = 1, ..., m1) with m1 a multiple of m2, and the remaining equations do not explicitly depend on y'(m1), ..., y'(n - 1), a substantial gain in computer time can be achived by setting the parameters M1 = m1 (> 0) and M2 = m2 (> 0) (m1 + m2 <= n). (default: M1 = 0, M2 = M1) [in] Nsmax,Nsmin (default) Minimum number of stages Nsmin and maximum number of stages Nsmax (= 3, 5 or 7). (default: Nsmin = 3, Nsmax = 7) The number of stages ns correspond to implicit Runge-Kutta method of order 2*ns - 1. [in] Ns (default) Value of ns for the first step. (3, 5 or 7) (default = Nsmin) [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) Maximaum step size. (default = Abs(Tend - T)) [in] Thet (Optional) Decides whether the Jacobian should be recomputed. (Thet < 1) (default = 0.001) Increase this value (e.g. to 0.1) when Jacobian evaluations are costly. For small systems, this value should be smaller (e.g. 0.001). Negative value forces to compute the Jacobian after every accepted step. [in] Facl,Facr (Optional) Parameters for step size selection. (Facl <= 1, Facr >= 1. (default Facl = 0.2, Facr = 8) The new step size is chosen subject to the restriction Facl <= hnew/hold <= Facr. [in] Safe (Optional) The safety factor in step size prediction (0.001 < Safe < 1). (default = 0.9) [in] Quot1,Quot2 (Optional) If Quot1 < hnew/hold < Quot2, the step size is not changed. (Quot1 <= 1, Quot2 >= 1) (default: Quot1 = 1, Quot2 = 1.2) For large systems, Quot1 = 0.99 and Quot2 = 2 might be good. This saves, together with a large Thet, LU decompositions and computing time. For small systems default values might be good. [in] Vitu (Optional) Order is increased if the contractivity factor is larger than this value. (default = 0.002) [in] Vitd (Optional) Order is decreased if the contractivity factor is larger than this value. (default = 0.8) [in] Hhou,Hhod (Optional) Order is decreased only if the step size ratio safisfies Hhou <= hnew/h <= Hhod. (default: Hhou = 1.2, Hhod = 0.8) Reference (1) E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations II. Stiff and differential-algebraic Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1996) 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_Radaua_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, Mode As Long, Neval As Long, Info As Long Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, YYpd(N - 1, N - 1) As Double, IRev As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) Mode = 2 T = 0: Y(0) = 1: Y(1) = 2 Tend = 10 Info = 0 Do Tout = T + 1 IRev = 0 Do Call Radaua_r(N, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Info, TT, YY(), YYp(), YYpd(), IRev, Neval) If IRev = 1 Then Call F2(N, TT, YY(), YYp()) Loop While IRev <> 0 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.367879441175671 0.908181747032157 2 0.135335283238606 -0.280811553265263 3 4.97870683691699E-02 -0.940205428205659 4 1.83156388907472E-02 -0.635327981830992 5 6.73794700079265E-03 0.290400132486203 6 2.47875217861059E-03 0.962649036304768 7 9.11881980026478E-04 0.75481411025709 8 3.3546262846512E-04 -0.145164571251474 9 1.23409804295512E-04 -0.91100685211013 10 4.53999298582475E-05 -0.839026129195407 270 0