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

Radaup_r() Sub Radaup_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Cont() As  Double, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, YYpd() As  Double, Irtrn As  Long, IRev As  Long, Optional Ns As  Long, Optional Iout 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 Neval As  Long, Optional Njac As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional M1 As  Long, Optional M2 As  Long, Optional Hinit As  Double, Optional Hmax As  Double, Optional MaxIter As  Long, Optional MaxNit As  Long, Optional StartN As  Long, Optional Nind1 As  Long, Optional Nind2 As  Long, Optional Nind3 As  Long, Optional Pred As  Long, Optional Hess As  Long, Optional Safe As  Double, Optional Thet As  Double, Optional FNewt As  Double, Optional Quot1 As  Double, Optional Quot2 As  Double, Optional Facl As  Double, Optional Facr As  Double, Optional Cnt As  Long  ) Initial value problem of ordinary differential equations (5, 9, 13-th order implicit Runge-Kutta method (Radau IIA)) (reverse communication version) NOTE - THIS PROGRAM IS DEPRECATED AND WILL BE REMOVED IN THE NEXT VERSION. Purpose This routine 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 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. M is the mass matrix. The system can be linearly implicit (M is not I (identity matrix)) or explicit (M = I). Radaup is the implicit Runge-Kutta code based on the 5, 9 or 13-th order Radau IIA method. It is provided with the step control algorithm and the continuous output feature. The order can be selected by setting Ns. In the case of 5th order (default), Radaup is equivalent to Radau5 but it will be slightly slower than Radau5. Radaup refers to the coefficient table during computation so that the order can be switched. On the other hand, radau5 is hard coded and dedicated to 5th order. See for details in the reference below. Radaup_r is the reverse communication version of Radaup. 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. [in] Initial value of the independent variable T. [out] Last value of the independent variable T of the final step (normally equals to Tout). The solution was successfully advanced to this point. It is possible to continue the integration to new point by recalling this routine with the new Tout value with setting Info = 1. [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 (normally equals to Tout). [in] Tout Set Tout to the point at which a solution is desired. Integration either forward in T (Tout > T) or backward in T (Tout < T) is permitted. The routine advances the solution from T to Tout using step sizes which are automatically selected so as to achieve the desired accuracy. [in] RTol() Array 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 a vector (LRTol = N). If LRTol = 2, ... or N-1, LRTol = 1 is assumed. If LRTol > N, LRTol = N is assumed. Even if LRTol = N, it is assumed to be 1 if LATol = 1. The tolerances are used by the code in a local error test at each step which requires roughly that   abs(local error of Y(i)) <= RTol(i)*abs(Y(i)) + ATol(i) for each component of Y() (i = 0 to LRTol-1). Setting RTol(i) = 0 results in a pure absolute error test on that component. RTol(i) and ATol(i) should not be zero at the same time (i = 0 to LRTol-1). [in] ATol() Array 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 a vector (LATol = N). If LATol = 2, ... or N-1, LATol = 1 is assumed. If LATol > N, LATol = N is assumed. Even if LATol = N, it is assumed to be 1 if LRTol = 1. The tolerances are used by the code in a local error test at each step which requires roughly that   abs(local error of Y(i)) <= RTol(i)*abs(Y(i)) + ATol(i) for each component of Y() (i = 0 to LATol-1). Setting ATol(i) = 0 results in a pure relative error test on that component. RTol(i) and ATol(i) should not be zero at the same time (i = 0 to LRTol-1). [in,out] Cont() Array Cont(LCont - 1) (LCont >= (Ns + 1)*N) Control information for dense output. (Required even if Iout = 0) [in,out] Info [in] = 0: Initialize and start computation (Solve new problem). = 1: Continue computation with new Tout value (Resume computation of previous call). [out] = -1: The argument N had an illegal value. (N < 1) = -3: The argument Y() is invalid. = -5: The argument RTol() had an illegal value. (RTol(i) < 0, RTol(i) = 0 and ATol(i) = 0) = -6: The argument ATol() had an illegal value. (ATol(i) < 0) = -7: The argument Cont() is invalid. = -8: the argument Info had an illegal value (Info <> 0 and Info <> 1) = -10: The argument YY() is invalid. = -11: The argument YYp() is invalid. = -12: The argument YYpd() is invalid. = -14: the argument IRev had an illegal value. (IRev <> 0, 1, 2, 4 nor 5) = -20: The argument Mljac or Mlmas had an illegal value. (Mlmas > Mljac) = -21: The argument Mujac or Mumas had an illegal value. (Mumas > Mujac) = -28: The argument M1 had an illegal value. (M1 < 0) = -29: The argument M1 or M2 had an illegal value. (M2 < 0 or M1 + M2 > N) = -35: The argument Nind1, Nind2 or Nind3 had an illegal value. (Nind1 + Nind2 + Nind3 <> N) = -42: The argument FNewt had an illegal value. (FNewt too small) = 1: Successful exit. = 2: Interrupted by Irtrn (normal return). = 11: Maximum number of steps exceeded. = 12: Step size becomes too small. = 13: Matrix is repeatedly singular. [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] YYpd() Array YYpd(LYYpd1 - 1, LYYpd2 - 1) (LYYpd1 >= max(Ljac, Lmas), LYYpd2 >= N) IRev = 2: The Jacobian (dfi/dyj) at TT and YY() should be provided in the next call.   Note - If Mljac = N, a matrix is stored as N x N full matrix (Ljac = N). If Mljac < N, a matrix is stored in band matrix form (Ljac = Mljac + Mujac + 1). If Ijac = 0, Ljac is 0. IRev = 4: The mass matrix M should be provided in the next call.   Note - If Mlmas = N, a matrix is stored as N x N full matrix (Lmas = N). If Mlmas < N, a matrix is stored in band matrix form (Lmas = Mlmas + Mumas + 1). If Imas = 0, Lmas is 0. [in,out] Irtrn [in] IRev = 5: Do not alter unless user want to interrupt the integration. If Irtrn is set to the negative value, the integration will be interrupted and exit with Info = 2. [out] IRev = 5: Returns 0, 1 or 2 in the first, intermediate or last return with IRev = 5, respectively. [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(). = 2: User should set the Jacobian (dfi/dxj) at TT and YY() in YYpd(). Do not alter other variables. = 4: User should set the mass matrix M in YYpd(). Do not alter other variables. = 5: User may output the intermediate result. Do not alter any variables. (See Iout) [in] Ns (Optional) The number of stages. (Ns = 3, 5 or 7) (default = 7) Ns values correspond to implicit Runge-Kutta method of order 5, 9 or 13, respectively. (For other values, the default value will be used) [in] Iout (Optional) Specifies if the intermediate result output is required. (default = 0) = 0: Output is not required. (Not return with IRev = 5) = 1: Returns after every successful step with IRev = 5 to output the intermediate results. This is same as calling Solout in the case of normal Radaup (non reverse communication version). The corresponding information is as follows. Nr = Naccept + 1, Told = previous T, T = current T, Y() = current Y(). Cont() is used in the same way with Solout for the dense output. Y(i) = Contrp_r(i, T2, Cont()) (If other value is specified for Iout, Iout = 0 will be assumed) [in] Ijac (Optional) Switch for calculation of Jacobian. (default = 0) = 0: Jacobian is calculated by finite differences. (Not return with IRev = 2) = 1: Jacobian is calculated by user externally. (Return with IRev = 2 when calculation is required) (For other values, the default value will be used) [in] Mljac (Optional) The lower bandwidth of Jacobian. (0 <= Mljac <= N) (default = N) If Mljac = N, Jacobian is stored as N x N full matrix. If Mljac < N, Jacobian is stored in band matrix form. (If Mljac < 0 or Mljac > N, the default value will be used) [in] Mujac (Optional) The upper bandwidth of Jacobian. (0 <= Mujac <= N) (default = 0) If Mljac = N, Mujac is ignored. (If Mujac < 0 or Mujac > N, the default value will be used) [in] Imas (Optional) Switch for calculation of mass matrix M. (default = 0) = 0: M is supposed to be the identity matrix. (Not return with IRev = 4) = 1: M is calculated by user externally. (Return with IRev = 4 when calculation is required). (For other values, the default value will be used) [in] Mlmas (Optional) The lower bandwidth of mass matrix M. (0 <= Mlmas <= N) (default = N) If Mlmas = N, M is stored as N x N full matrix. If Mlmas < N, M is stored in band matrix form. (If Mlmas < 0 or Mlmas > N, the default value will be used) [in] Mumas (Optional) The upper bandwidth of mass matrix M. (0 <= Mumas <= N) (default = 0) If Mlmas = N, Mumas is ignored. (If Mumas < 0 or Mumas > N, the default value will be used) [out] Neval (Optional) Number of function evaluations. (Those for Jacobian evaluations are not included) [out] Njac (Optional) Number of Jacobian evaluations. (Those by finite differences are included) [out] Nstep (Optional) Number of computed steps. [out] Naccept (Optional) Number of accepted steps. [out] Nreject (Optional) Number of rejected steps. (Step rejections in the first step are not counted) [in] M1,M2 (Optional) If the first M1 equations has the following form   y'(i) = y(i + M2) for i = 1 to M1, with M1 a multiple of M2, and the remaining equations do not explicitly depend on y'(M1), ..., y'(N-1), efficient computation can be achieved by setting parameters M1 and M2 to nonzero values. (M1 > 0, M2 > 0, M1 + M2 <= N) (default M1 = M2 = 0) When parameters are set to nonzero, only the elements of non-trivial part of the Jacobian (rows M1+1 to N) heve to be stored in (N - M1) x N array. Also only the elements of right lower block of order N - M1 of the mass matrix M have to be stored in (N - M1) x (N - M1) array. [in] Hinit (Optional) Initial step size. (default = 1.0e-6) H = 1/||f'||, usually 1.0e-3 or 1.0e-5 is good for stiff equations with initial transient. (If Hinit = 0, 1.0e-6 will be used) [in] Hmax (Optional) Maximal step size. (default = Tout - T) (If Hmax = 0, the default value will be used) [in] MaxIter (Optional) Maximum number of allowed steps. (default = 100000) (If MaxIter <= 0, the default value will be used) [in] MaxNit (Optional) Maximum number of Newton iterations for the solution of the implicit system in each step. (default = 7) (if MaxNit <= 0, the default value will be used) [in] StartN (Optional) Starting value for Newton's method. (default = 0) = 0: Extrapolated collocation solution. <> 0: Zero starting value. The latter is recommended if Newton's method has difficulties with convergence. [in] Nind1,Nind2,Nind3 (Optional) Dimension of the index 1, 2 and 3 variables, respectively. (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. For ordinary defferential equations (ODEs), Nind1 equals to the number of equations of the system. [in] Pred (Optional) Switch for step size strategy. (default = 1) = 1: Model predictive controller (Gustafsson). = 2: Classical step size control. (For other values, the default value will be assumed) [in] Hess (Optional) Switch to transform the Jacobian matrix to Hessenberg form. (default = 0) = 0: Do not transform. (This option is not supported. For other values, 0 will be assumed.) This is particularly advantageous for large systems with full Jacobian. It does not work for banded Jacobian and not for implicit systems (Mas <> NullPtr). Parameters [in] Safe (Optional) The safety factor in step size prediction. (default = 0.9) (If Safe <= 0.001 or Safe >= 1, the default value will be used) [in] Thet (Optional) Decides whether the Jacobian should be recomputed. (Thet < 1) (default = 0.001) Increase Thet (e.g. 0.1) when Jacobian evaluations are costly. For small systems, Thet should be smaller (e.g. 0.001). Negative Thet forces the code to compute the Jacobian after every accepted step. (If Thet = 0 or Thet >= 1, the default value will be used) [in] FNewt (Optional) Stopping criterion for Newton's method. (Usually FNewt < 1) (default = automatically computed from RTol(0)) Smaller values make the code slower, but safer. (If FNewt = 0, the default value will be used) [in] Quot1,Quot2 (Optional) If Quot1 < Hnew/Hold < Quot2, the step size is not chenged. (Quot1 <= 1, Quot2 >= 1) (default: Quot1 = 1, Quot2 = 1.2) Quot1 = 1 and Quot2 = 1.2 for small systems, Quot1 = 0.99 and Quot2 = 2 for large systems might be good. (If Quot1 = 0 or Quot1 > 1, Quot2 = 0 or Quot2 < 1, the default values will be used respectively) [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. (If Facl = 0 or Facl > 1, Facr = 0 or Facr < 1, the default values will be used respectively) [in] Cnt (Optional) Specifies when Neval, Njac, Nstep, Naccept and Nreject are reset to zero. (default = 0) = 0: Reset whenever this routine is called. <> 0: Reset only if this routine is called with Info = 0. Reference 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 (1) 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_Radaup_r() Const N = 2, Ns = 7 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 Cont((Ns + 1) * N) As Double Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, YYpd(0) As Double Dim Irtrn As Long, IRev As Long, Iout 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 Radaup_r(N, T, Y(), Tout, RTol(), ATol(), Cont(), Info, TT, YY(), YYp(), YYpd(), Irtrn, IRev) If IRev = 1 Then Call F2(N, TT, YY(), YYp()) Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Radaup_r: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub Example Results 1 0.367879441171439 0.90818174704269 2 0.135335283236607 -0.280811553308553 3 4.97870683678812E-02 -0.940205428272677 4 1.83156388887558E-02 -0.635327982020324 5 6.73794699908958E-03 0.29040013245329 6 2.47875217664834E-03 0.962649038862733 7 9.11881965530575E-04 0.754814136356455 8 3.35462627894578E-04 -0.145164571164979 9 1.23409804101964E-04 -0.911006852111188 10 4.53999297868879E-05 -0.839026129195484 Example Program (2) Solve the following initial value problem of ordinary differential equations (stiff problem) (using dense output). 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_Radaup_r_2() Const N = 2, Ns = 7 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 Cont((Ns + 1) * N) As Double Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, YYpd(0) As Double Dim Irtrn As Long, IRev As Long, Iout As Long Dim Y0 As Double, Y1 As Double RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) Iout = 1 T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Tout = 1 Info = 0 IRev = 0 Do Call Radaup_r(N, T, Y(), Tend, RTol(), ATol(), Cont(), Info, TT, YY(), YYp(), YYpd(), Irtrn, IRev, , Iout) If IRev = 1 Then Call F2(N, TT, YY(), YYp()) ElseIf IRev = 5 Then While T >= Tout Y0 = Contrp_r(0, Tout, Cont()) Y1 = Contrp_r(1, Tout, Cont()) Debug.Print Tout, Y0, Y1 Tout = Tout + 1 Wend End If Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Radaup_r: Info =", Info End Sub Example Results 1 0.367879434499153 0.908181731594218 2 0.135335282851716 -0.280811562130519 3 4.97870677636028E-02 -0.940205433642388 4 1.83156366415058E-02 -0.635327957933275 5 6.73794563215949E-03 0.290399971657046 6 2.47874989741329E-03 0.962649104070445 7 9.11882397527037E-04 0.754814137911431 8 3.35462430575541E-04 -0.145163891518891 9 1.23409474270825E-04 -0.911006209259385 10 4.53999295433717E-05 -0.839026129255171