diff --git a/CHANGELOG.md b/CHANGELOG.md index 53a51e0..287d59c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,7 +17,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added "Lint" GitHub Action to check other actions for security issues - Added new example files: `rkadj.kpp`, `sd4.kpp`, `sdtlm.kpp` - Added new C-I tests: `F90_rkadj`, `F90_sd4`, `F90_sdtlm` - +- Added `Rodas3_1` integration method to `rosenbrock_autoreduce.f90` and `rosenbrock.f90`, with updated coefficents by @msl3v +- Added `Rodas3_1` integration method to `rosenbrock_autoreduce.f90` and `rosenbrock.f90`, with updated coefficents by @msl3v +- Added `Rodas3_1` integration method to `rosenbrock_autoreduce.f90`, `rosenbrock.f90`, `rosenbrock.c`, and `rosenbrock.m`, with updated coefficents by @msl3v +- Added `Rodas3_1` integration method (with updated coefficents by @msl3v) to `rosenbrock_autoreduce.f90`, `rosenbrock.{c,f90,m}`, `rosenbrock_adj.{c,f90}` + ### Changed - Updated ReadTheDocs documentation to reflect that C-I tests are now done as a GitHub Action @@ -26,6 +30,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Updated rules to ignore files in `.gitignore` and updated comments accordingly - Fixed a bug that prevented `.ci-pipelines/ci-cleanup-script.sh` from removing KPP-generated files for MCM mechanisms - Fixed typo in error message in `int/rosenbrock_autoreduce.f90` +- Fixed incorrect comment in `int/rosenbrock_autoreduce.f90`: `10*H` should be `0.1*H` ### Removed - Removed C-I tests on Microsoft Azure Dev Pipelines @@ -521,3 +526,4 @@ Changes by Domenico Taraborrelli: - examples/mimi* deleted ## [2.1] - 2005-07-19 + diff --git a/int/rosenbrock.c b/int/rosenbrock.c index 27eb893..f54a7f3 100644 --- a/int/rosenbrock.c +++ b/int/rosenbrock.c @@ -1,7 +1,7 @@ #define MAX(a,b) ( ((a) >= (b)) ?(a):(b) ) - #define MIN(b,c) ( ((b) < (c)) ?(b):(c) ) - #define ABS(x) ( ((x) >= 0 ) ?(x):(-x) ) + #define MIN(b,c) ( ((b) < (c)) ?(b):(c) ) + #define ABS(x) ( ((x) >= 0 ) ?(x):(-x) ) #define SQRT(d) ( pow((d),0.5) ) #define SIGN(x) ( ((x) >= 0 ) ?[0]:(-1) ) @@ -9,47 +9,47 @@ #define ZERO (KPP_REAL)0.0 #define ONE (KPP_REAL)1.0 #define HALF (KPP_REAL)0.5 - #define DeltaMin (KPP_REAL)1.0e-6 - -/*~~~> Collect statistics: global variables */ + #define DeltaMin (KPP_REAL)1.0e-6 + +/*~~~> Collect statistics: global variables */ int Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng; -/*~~~> Function headers */ - void FunTemplate(KPP_REAL, KPP_REAL [], KPP_REAL []); +/*~~~> Function headers */ + void FunTemplate(KPP_REAL, KPP_REAL [], KPP_REAL []); void JacTemplate(KPP_REAL, KPP_REAL [], KPP_REAL []) ; int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL AbsTol[], KPP_REAL RelTol[], - void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), + void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []), KPP_REAL RPAR[], int IPAR[]); int RosenbrockIntegrator( - KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend , + KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend , KPP_REAL AbsTol[], KPP_REAL RelTol[], - void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), + void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []), int ros_S, - KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_M[], KPP_REAL ros_E[], KPP_REAL ros_A[], KPP_REAL ros_C[], KPP_REAL ros_Alpha[],KPP_REAL ros_Gamma[], KPP_REAL ros_ELO, char ros_NewF[], - char Autonomous, char VectorTol, int Max_no_steps, + char Autonomous, char VectorTol, int Max_no_steps, KPP_REAL Roundoff, KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart, - KPP_REAL FacMin, KPP_REAL FacMax, KPP_REAL FacRej, KPP_REAL FacSafe, - KPP_REAL *Texit, KPP_REAL *Hexit ); + KPP_REAL FacMin, KPP_REAL FacMax, KPP_REAL FacRej, KPP_REAL FacSafe, + KPP_REAL *Texit, KPP_REAL *Hexit ); char ros_PrepareMatrix ( - KPP_REAL* H, - int Direction, KPP_REAL gam, KPP_REAL Jac0[], + KPP_REAL* H, + int Direction, KPP_REAL gam, KPP_REAL Jac0[], KPP_REAL Ghimj[], int Pivot[] ); - KPP_REAL ros_ErrorNorm ( - KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[], - KPP_REAL AbsTol[], KPP_REAL RelTol[], + KPP_REAL ros_ErrorNorm ( + KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[], + KPP_REAL AbsTol[], KPP_REAL RelTol[], char VectorTol ); int ros_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H); - void ros_FunTimeDerivative ( - KPP_REAL T, KPP_REAL Roundoff, - KPP_REAL Y[], KPP_REAL Fcn0[], - void ode_Fun(KPP_REAL, KPP_REAL [], KPP_REAL []), + void ros_FunTimeDerivative ( + KPP_REAL T, KPP_REAL Roundoff, + KPP_REAL Y[], KPP_REAL Fcn0[], + void ode_Fun(KPP_REAL, KPP_REAL [], KPP_REAL []), KPP_REAL dFdT[] ); void Fun( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] ); void Jac_SP( KPP_REAL Y[], KPP_REAL FIX[], KPP_REAL RCONST[], KPP_REAL Ydot[] ); @@ -61,31 +61,35 @@ void WAXPY(int N, KPP_REAL Alpha, KPP_REAL X[], int incX, KPP_REAL Y[], int incY ); void WSCAL(int N, KPP_REAL Alpha, KPP_REAL X[], int incX); KPP_REAL WLAMCH( char C ); - void Ros2 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + void Ros2 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); - void Ros3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + void Ros3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); - void Ros4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + void Ros4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); - void Rodas3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + void Rodas3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); - void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + void Rodas3_1 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); + void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ); int KppDecomp( KPP_REAL A[] ); void KppSolve ( KPP_REAL A[], KPP_REAL b[] ); void Update_SUN(); void Update_RCONST(); - + /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ void INTEGRATE( KPP_REAL TIN, KPP_REAL TOUT ) /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ @@ -98,8 +102,8 @@ void INTEGRATE( KPP_REAL TIN, KPP_REAL TOUT ) IPAR[i] = 0; RPAR[i] = ZERO; } /* for */ - - + + IPAR[0] = 0; /* non-autonomous */ IPAR[1] = 1; /* vector tolerances */ RPAR[2] = STEPMIN; /* starting step */ @@ -110,7 +114,7 @@ void INTEGRATE( KPP_REAL TIN, KPP_REAL TOUT ) &FunTemplate, &JacTemplate, RPAR, IPAR); - + Ns=Ns+IPAR[12]; Na=Na+IPAR[13]; Nr=Nr+IPAR[14]; @@ -122,21 +126,21 @@ void INTEGRATE( KPP_REAL TIN, KPP_REAL TOUT ) if (IERR < 0) printf("\n Rosenbrock: Unsucessful step at T=%g: IERR=%d\n", TIN,IERR); - + TIN = RPAR[10]; /* Exit time */ STEPMIN = RPAR[11]; /* Last step */ - + } /* INTEGRATE */ /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, KPP_REAL AbsTol[], KPP_REAL RelTol[], - void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), + void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []), KPP_REAL RPAR[], int IPAR[]) /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - + Solves the system y'=F(t,y) using a Rosenbrock method defined by: G = 1/(H*gamma[0]) - ode_Jac(t0,Y0) @@ -144,42 +148,42 @@ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j G * K_i = ode_Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j + gamma(i)*dF/dT(t0, Y0) - Y1 = Y0 + \sum_{j=1}^S M(j)*K_j + Y1 = Y0 + \sum_{j=1}^S M(j)*K_j For details on Rosenbrock methods and their implementation consult: E. Hairer and G. Wanner "Solving ODEs II. Stiff and differential-algebraic problems". - Springer series in computational mathematics, Springer-Verlag, 1996. - The codes contained in the book inspired this implementation. + Springer series in computational mathematics, Springer-Verlag, 1996. + The codes contained in the book inspired this implementation. (C) Adrian Sandu, August 2004 - Virginia Polytechnic Institute and State University + Virginia Polytechnic Institute and State University Contact: sandu@cs.vt.edu This implementation is part of KPP - the Kinetic PreProcessor ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - *~~~> INPUT ARGUMENTS: - + + *~~~> INPUT ARGUMENTS: + - Y(NVAR) = vector of initial conditions (at T=Tstart) - [Tstart,Tend] = time range of integration - (if Tstart>Tend the integration is performed backwards in time) + (if Tstart>Tend the integration is performed backwards in time) - RelTol, AbsTol = user precribed accuracy -- void ode_Fun( T, Y, Ydot ) = ODE function, - returns Ydot = Y' = F(T,Y) +- void ode_Fun( T, Y, Ydot ) = ODE function, + returns Ydot = Y' = F(T,Y) - void ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function, - returns Jcb = dF/dY + returns Jcb = dF/dY - IPAR(1:10) = int inputs parameters - RPAR(1:10) = real inputs parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - *~~~> OUTPUT ARGUMENTS: - + *~~~> OUTPUT ARGUMENTS: + - Y(NVAR) -> vector of final states (at T->Tend) - IPAR(11:20) -> int output parameters - RPAR(11:20) -> real output parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - *~~~> RETURN VALUE (int): + *~~~> RETURN VALUE (int): - IERR -> job status upon return - succes (positive value) or failure (negative value) - @@ -193,7 +197,7 @@ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, = -7 : Step size too small = -8 : Matrix is repeatedly singular ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - + *~~~> INPUT PARAMETERS: Note: For input parameters equal to zero the default values of the @@ -209,24 +213,25 @@ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, IPAR[3] -> selection of a particular Rosenbrock method = 0 : default method is Rodas3 = 1 : method is Ros2 - = 2 : method is Ros3 - = 3 : method is Ros4 + = 2 : method is Ros3 + = 3 : method is Ros4 = 4 : method is Rodas3 - = 5: method is Rodas4 + = 5 : method is Rodas4 + = 7 : method is Rodas3.1 RPAR[0] -> Hmin, lower bound for the integration step size - It is strongly recommended to keep Hmin = ZERO + It is strongly recommended to keep Hmin = ZERO RPAR[1] -> Hmax, upper bound for the integration step size RPAR[2] -> Hstart, starting value for the integration step size - + RPAR[3] -> FacMin, lower bound on step decrease factor (default=0.2) RPAR[4] -> FacMin,upper bound on step increase factor (default=6) RPAR[5] -> FacRej, step decrease factor after multiple rejections (default=0.1) - RPAR[6] -> FacSafe, by which the new step is slightly smaller + RPAR[6] -> FacSafe, by which the new step is slightly smaller than the predicted value (default=0.9) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - + *~~~> OUTPUT PARAMETERS: Note: each call to Rosenbrock adds the corrent no. of fcn calls @@ -242,21 +247,21 @@ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, IPAR[16] = No. of forward/backward substitutions IPAR[17] = No. of singular matrix decompositions - RPAR[10] -> Texit, the time corresponding to the + RPAR[10] -> Texit, the time corresponding to the computed Y upon return RPAR[11] -> Hexit, last accepted step before exit - For multiple restarts, use Hexit as Hstart in the following run + For multiple restarts, use Hexit as Hstart in the following run ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ -{ +{ - /*~~~> The method parameters */ + /*~~~> The method parameters */ #define Smax 6 int Method, ros_S; KPP_REAL ros_M[Smax], ros_E[Smax]; KPP_REAL ros_A[Smax*(Smax-1)/2], ros_C[Smax*(Smax-1)/2]; KPP_REAL ros_Alpha[Smax], ros_Gamma[Smax], ros_ELO; char ros_NewF[Smax], ros_Name[12]; - /*~~~> Local variables */ + /*~~~> Local variables */ int Max_no_steps, IERR, i, UplimTol; char Autonomous, VectorTol; KPP_REAL Roundoff,FacMin,FacMax,FacRej,FacSafe; @@ -272,7 +277,7 @@ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, Ndec = IPAR[15]; Nsol = IPAR[16]; Nsng = IPAR[17]; - + /*~~~> Autonomous or time dependent ODE. Default is time dependent. */ Autonomous = !(IPAR[0] == 0); @@ -280,54 +285,54 @@ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, ! For Vector tolerances (IPAR[1] == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR) */ if (IPAR[1] == 0) { VectorTol = 1; UplimTol = KPP_NVAR; - } else { + } else { VectorTol = 0; UplimTol = 1; } /* end if */ - + /*~~~> The maximum number of steps admitted */ - if (IPAR[2] == 0) + if (IPAR[2] == 0) Max_no_steps = 100000; - else + else Max_no_steps=IPAR[2]; - if (Max_no_steps < 0) { + if (Max_no_steps < 0) { printf("\n User-selected max no. of steps: IPAR[2]=%d\n",IPAR[2]); return ros_ErrorMsg(-1,Tstart,ZERO); } /* end if */ /*~~~> The particular Rosenbrock method chosen */ - if (IPAR[3] == 0) + if (IPAR[3] == 0) Method = 3; - else + else Method = IPAR[3]; - if ( (IPAR[3] < 1) || (IPAR[3] > 5) ){ + if ( (IPAR[3] < 1) || (IPAR[3] > 5) ){ printf("\n User-selected Rosenbrock method: IPAR[3]=%d\n",IPAR[3]); return ros_ErrorMsg(-2,Tstart,ZERO); } /* end if */ - + /*~~~> Unit Roundoff (1+Roundoff>1) */ Roundoff = WLAMCH('E'); /*~~~> Lower bound on the step size: (positive value) */ Hmin = RPAR[0]; - if (RPAR[0] < ZERO) { + if (RPAR[0] < ZERO) { printf("\n User-selected Hmin: RPAR[0]=%e\n", RPAR[0]); return ros_ErrorMsg(-3,Tstart,ZERO); } /* end if */ /*~~~> Upper bound on the step size: (positive value) */ - if (RPAR[1] == ZERO) + if (RPAR[1] == ZERO) Hmax = ABS(Tend-Tstart); - else + else Hmax = MIN(ABS(RPAR[1]),ABS(Tend-Tstart)); - if (RPAR[1] < ZERO) { + if (RPAR[1] < ZERO) { printf("\n User-selected Hmax: RPAR[1]=%e\n", RPAR[1]); return ros_ErrorMsg(-3,Tstart,ZERO); } /* end if */ /*~~~> Starting step size: (positive value) */ - if (RPAR[2] == ZERO) + if (RPAR[2] == ZERO) Hstart = MAX(Hmin,DeltaMin); else Hstart = MIN(ABS(RPAR[2]),ABS(Tend-Tstart)); - if (RPAR[2] < ZERO) { + if (RPAR[2] < ZERO) { printf("\n User-selected Hstart: RPAR[2]=%e\n", RPAR[2]); return ros_ErrorMsg(-3,Tstart,ZERO); } /* end if */ @@ -336,33 +341,33 @@ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, FacMin = (KPP_REAL)0.2; else FacMin = RPAR[3]; - if (RPAR[3] < ZERO) { + if (RPAR[3] < ZERO) { printf("\n User-selected FacMin: RPAR[3]=%e\n", RPAR[3]); return ros_ErrorMsg(-4,Tstart,ZERO); } /* end if */ - if (RPAR[4] == ZERO) + if (RPAR[4] == ZERO) FacMax = (KPP_REAL)6.0; else FacMax = RPAR[4]; - if (RPAR[4] < ZERO) { + if (RPAR[4] < ZERO) { printf("\n User-selected FacMax: RPAR[4]=%e\n", RPAR[4]); return ros_ErrorMsg(-4,Tstart,ZERO); } /* end if */ /*~~~> FacRej: Factor to decrease step after 2 succesive rejections */ - if (RPAR[5] == ZERO) + if (RPAR[5] == ZERO) FacRej = (KPP_REAL)0.1; else FacRej = RPAR[5]; - if (RPAR[5] < ZERO) { + if (RPAR[5] < ZERO) { printf("\n User-selected FacRej: RPAR[5]=%e\n", RPAR[5]); return ros_ErrorMsg(-4,Tstart,ZERO); } /* end if */ /*~~~> FacSafe: Safety Factor in the computation of new step size */ - if (RPAR[6] == ZERO) + if (RPAR[6] == ZERO) FacSafe = (KPP_REAL)0.9; else FacSafe = RPAR[6]; - if (RPAR[6] < ZERO) { + if (RPAR[6] < ZERO) { printf("\n User-selected FacSafe: RPAR[6]=%e\n", RPAR[6]); return ros_ErrorMsg(-4,Tstart,ZERO); } /* end if */ @@ -375,47 +380,51 @@ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, return ros_ErrorMsg(-5,Tstart,ZERO); } /* end if */ } /* for */ - - + + /*~~~> Initialize the particular Rosenbrock method */ switch (Method) { case 1: - Ros2(&ros_S, ros_A, ros_C, ros_M, ros_E, + Ros2(&ros_S, ros_A, ros_C, ros_M, ros_E, ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); - break; + break; case 2: - Ros3(&ros_S, ros_A, ros_C, ros_M, ros_E, + Ros3(&ros_S, ros_A, ros_C, ros_M, ros_E, ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); - break; + break; case 3: - Ros4(&ros_S, ros_A, ros_C, ros_M, ros_E, + Ros4(&ros_S, ros_A, ros_C, ros_M, ros_E, ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); - break; + break; case 4: - Rodas3(&ros_S, ros_A, ros_C, ros_M, ros_E, + Rodas3(&ros_S, ros_A, ros_C, ros_M, ros_E, ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); - break; + break; case 5: - Rodas4(&ros_S, ros_A, ros_C, ros_M, ros_E, + Rodas4(&ros_S, ros_A, ros_C, ros_M, ros_E, + ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); + break; + case 7: + Rodas3_1(&ros_S, ros_A, ros_C, ros_M, ros_E, ros_Alpha, ros_Gamma, ros_NewF, &ros_ELO, ros_Name); - break; + break; default: printf("\n Unknown Rosenbrock method: IPAR[3]= %d", Method); - return ros_ErrorMsg(-2,Tstart,ZERO); + return ros_ErrorMsg(-2,Tstart,ZERO); } /* end switch */ /*~~~> Rosenbrock method */ IERR = RosenbrockIntegrator( Y,Tstart,Tend, AbsTol, RelTol, ode_Fun,ode_Jac , - /* Rosenbrock method coefficients */ - ros_S, ros_M, ros_E, ros_A, ros_C, + /* Rosenbrock method coefficients */ + ros_S, ros_M, ros_E, ros_A, ros_C, ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, - /* Integration parameters */ + /* Integration parameters */ Autonomous, VectorTol, Max_no_steps, Roundoff, Hmin, Hmax, Hstart, - FacMin, FacMax, FacRej, FacSafe, - /* Output parameters */ + FacMin, FacMax, FacRej, FacSafe, + /* Output parameters */ &Texit, &Hexit ); @@ -430,87 +439,87 @@ int Rosenbrock(KPP_REAL Y[], KPP_REAL Tstart, KPP_REAL Tend, IPAR[17] = Nsng; /*~~~> Last T and H */ RPAR[10] = Texit; - RPAR[11] = Hexit; - + RPAR[11] = Hexit; + return IERR; - + } /* Rosenbrock */ - + /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ int RosenbrockIntegrator( - /*~~~> Input: the initial condition at Tstart; Output: the solution at T */ + /*~~~> Input: the initial condition at Tstart; Output: the solution at T */ KPP_REAL Y[], - /*~~~> Input: integration interval */ - KPP_REAL Tstart, KPP_REAL Tend , - /*~~~> Input: tolerances */ + /*~~~> Input: integration interval */ + KPP_REAL Tstart, KPP_REAL Tend , + /*~~~> Input: tolerances */ KPP_REAL AbsTol[], KPP_REAL RelTol[], - /*~~~> Input: ode function and its Jacobian */ - void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), + /*~~~> Input: ode function and its Jacobian */ + void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), void (*ode_Jac)(KPP_REAL, KPP_REAL [], KPP_REAL []) , - /*~~~> Input: The Rosenbrock method parameters */ + /*~~~> Input: The Rosenbrock method parameters */ int ros_S, - KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_M[], KPP_REAL ros_E[], KPP_REAL ros_A[], KPP_REAL ros_C[], KPP_REAL ros_Alpha[],KPP_REAL ros_Gamma[], KPP_REAL ros_ELO, char ros_NewF[], - /*~~~> Input: integration parameters */ + /*~~~> Input: integration parameters */ char Autonomous, char VectorTol, - int Max_no_steps, + int Max_no_steps, KPP_REAL Roundoff, KPP_REAL Hmin, KPP_REAL Hmax, KPP_REAL Hstart, - KPP_REAL FacMin, KPP_REAL FacMax, KPP_REAL FacRej, KPP_REAL FacSafe, - /*~~~> Output: time at which the solution is returned (T=Tend if success) - and last accepted step */ - KPP_REAL *Texit, KPP_REAL *Hexit ) + KPP_REAL FacMin, KPP_REAL FacMax, KPP_REAL FacRej, KPP_REAL FacSafe, + /*~~~> Output: time at which the solution is returned (T=Tend if success) + and last accepted step */ + KPP_REAL *Texit, KPP_REAL *Hexit ) /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - Template for the implementation of a generic Rosenbrock method + Template for the implementation of a generic Rosenbrock method defined by ros_S (no of stages) and coefficients ros_{A,C,M,E,Alpha,Gamma} - - returned value: IERR, indicator of success (if positive) + + returned value: IERR, indicator of success (if positive) or failure (if negative) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ +{ KPP_REAL Ynew[KPP_NVAR], Fcn0[KPP_NVAR], Fcn[KPP_NVAR], dFdT[KPP_NVAR], Jac0[KPP_LU_NONZERO], Ghimj[KPP_LU_NONZERO]; - KPP_REAL K[KPP_NVAR*ros_S]; - KPP_REAL H, T, Hnew, HC, HG, Fac, Tau; + KPP_REAL K[KPP_NVAR*ros_S]; + KPP_REAL H, T, Hnew, HC, HG, Fac, Tau; KPP_REAL Err, Yerr[KPP_NVAR]; int Pivot[KPP_NVAR], Direction, ioffset, j, istage; char RejectLastH, RejectMoreH; /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ - + /*~~~> INITIAL PREPARATIONS */ T = Tstart; *Hexit = 0.0; - H = MIN(Hstart,Hmax); - if (ABS(H) <= 10.0*Roundoff) + H = MIN(Hstart,Hmax); + if (ABS(H) <= 10.0*Roundoff) H = DeltaMin; - + if (Tend >= Tstart) { Direction = +1; } else { Direction = -1; - } /* end if */ + } /* end if */ RejectLastH=0; RejectMoreH=0; - - /*~~~> Time loop begins below */ + + /*~~~> Time loop begins below */ while ( ( (Direction > 0) && ((T-Tend)+Roundoff <= ZERO) ) - || ( (Direction < 0) && ((Tend-T)+Roundoff <= ZERO) ) ) { - + || ( (Direction < 0) && ((Tend-T)+Roundoff <= ZERO) ) ) { + if ( Nstp > Max_no_steps ) { /* Too many steps */ *Texit = T; return ros_ErrorMsg(-6,T,H); - } + } if ( ((T+0.1*H) == T) || (H <= Roundoff) ) { /* Step size too small */ *Texit = T; return ros_ErrorMsg(-7,T,H); - } - - /*~~~> Limit H if necessary to avoid going beyond Tend */ + } + + /*~~~> Limit H if necessary to avoid going beyond Tend */ *Hexit = H; H = MIN(H,ABS(Tend-T)); @@ -518,16 +527,16 @@ int RosenbrockIntegrator( (*ode_Fun)(T,Y,Fcn0); /*~~~> Compute the function derivative with respect to T */ - if (!Autonomous) + if (!Autonomous) ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, ode_Fun, dFdT ); - + /*~~~> Compute the Jacobian at current time */ (*ode_Jac)(T,Y,Jac0); - + /*~~~> Repeat step calculation until current step accepted */ while (1) { /* WHILE STEP NOT ACCEPTED */ - + if( ros_PrepareMatrix( &H, Direction, ros_Gamma[0], Jac0, Ghimj, Pivot) ) { /* More than 5 consecutive failed decompositions */ *Texit = T; @@ -536,10 +545,10 @@ int RosenbrockIntegrator( /*~~~> Compute the stages */ for (istage = 1; istage <= ros_S; istage++) { - + /* Current istage offset. Current istage vector is K[ioffset:ioffset+KPP_NVAR-1] */ ioffset = KPP_NVAR*(istage-1); - + /* For the 1st istage the function has been computed previously */ if ( istage == 1 ) WCOPY(KPP_NVAR,Fcn0,1,Fcn,1); @@ -548,27 +557,27 @@ int RosenbrockIntegrator( WCOPY(KPP_NVAR,Y,1,Ynew,1); for (j = 1; j <= istage-1; j++) WAXPY(KPP_NVAR,ros_A[(istage-1)*(istage-2)/2+j-1], - &K[KPP_NVAR*(j-1)],1,Ynew,1); + &K[KPP_NVAR*(j-1)],1,Ynew,1); Tau = T + ros_Alpha[istage-1]*Direction*H; (*ode_Fun)(Tau,Ynew,Fcn); } /*end if ros_NewF(istage)*/ } /* end if istage */ - + WCOPY(KPP_NVAR,Fcn,1,&K[ioffset],1); for (j = 1; j <= istage-1; j++) { HC = ros_C[(istage-1)*(istage-2)/2+j-1]/(Direction*H); WAXPY(KPP_NVAR,HC,&K[KPP_NVAR*(j-1)],1,&K[ioffset],1); } /* for j */ - - if ((!Autonomous) && (ros_Gamma[istage-1])) { + + if ((!Autonomous) && (ros_Gamma[istage-1])) { HG = Direction*H*ros_Gamma[istage-1]; WAXPY(KPP_NVAR,HG,dFdT,1,&K[ioffset],1); } /* end if !Autonomous */ - + SolveTemplate(Ghimj, Pivot, &K[ioffset]); - - } /* for istage */ - + + } /* for istage */ + /*~~~> Compute the new solution */ WCOPY(KPP_NVAR,Y,1,Ynew,1); @@ -577,13 +586,13 @@ int RosenbrockIntegrator( /*~~~> Compute the error estimation */ WSCAL(KPP_NVAR,ZERO,Yerr,1); - for (j=1; j<=ros_S; j++) + for (j=1; j<=ros_S; j++) WAXPY(KPP_NVAR,ros_E[j-1],&K[KPP_NVAR*(j-1)],1,Yerr,1); Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ); /*~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax */ Fac = MIN(FacMax,MAX(FacMin,FacSafe/pow(Err,ONE/ros_ELO))); - Hnew = H*Fac; + Hnew = H*Fac; /*~~~> Check the error magnitude and adjust step size */ Nstp++; @@ -593,46 +602,46 @@ int RosenbrockIntegrator( T += Direction*H; Hnew = MAX(Hmin,MIN(Hnew,Hmax)); /* No step size increase after a rejected step */ - if (RejectLastH) - Hnew = MIN(Hnew,H); + if (RejectLastH) + Hnew = MIN(Hnew,H); RejectLastH = 0; RejectMoreH = 0; H = Hnew; break; /* EXIT THE LOOP: WHILE STEP NOT ACCEPTED */ } else { /*~~~> Reject step */ - if (Nacc >= 1) - Nrej++; - if (RejectMoreH) - Hnew=H*FacRej; + if (Nacc >= 1) + Nrej++; + if (RejectMoreH) + Hnew=H*FacRej; RejectMoreH = RejectLastH; RejectLastH = 1; H = Hnew; } /* end if Err <= 1 */ } /* while LOOP: WHILE STEP NOT ACCEPTED */ - } /* while: time loop */ - + } /* while: time loop */ + /*~~~> The integration was successful */ *Texit = T; - return 1; + return 1; } /* RosenbrockIntegrator */ - - + + /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -KPP_REAL ros_ErrorNorm ( - /*~~~> Input arguments */ - KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[], - KPP_REAL AbsTol[], KPP_REAL RelTol[], +KPP_REAL ros_ErrorNorm ( + /*~~~> Input arguments */ + KPP_REAL Y[], KPP_REAL Ynew[], KPP_REAL Yerr[], + KPP_REAL AbsTol[], KPP_REAL RelTol[], char VectorTol ) /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Computes and returns the "scaled norm" of the error vector Yerr -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ - /*~~~> Local variables */ - KPP_REAL Err, Scale, Ymax; +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ + /*~~~> Local variables */ + KPP_REAL Err, Scale, Ymax; int i; - + Err = ZERO; for (i=0; i Input arguments: */ - KPP_REAL T, KPP_REAL Roundoff, - KPP_REAL Y[], KPP_REAL Fcn0[], - void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), - /*~~~> Output arguments: */ +void ros_FunTimeDerivative ( + /*~~~> Input arguments: */ + KPP_REAL T, KPP_REAL Roundoff, + KPP_REAL Y[], KPP_REAL Fcn0[], + void (*ode_Fun)(KPP_REAL, KPP_REAL [], KPP_REAL []), + /*~~~> Output arguments: */ KPP_REAL dFdT[] ) /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The time partial derivative of the function by finite differences -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ { - /*~~~> Local variables */ - KPP_REAL Delta; - + /*~~~> Local variables */ + KPP_REAL Delta; + Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)); (*ode_Fun)(T+Delta,Y,dFdT); WAXPY(KPP_NVAR,(-ONE),Fcn0,1,dFdT,1); @@ -673,13 +682,13 @@ void ros_FunTimeDerivative ( } /* ros_FunTimeDerivative */ -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ char ros_PrepareMatrix ( - /* Inout argument: (step size is decreased when LU fails) */ - KPP_REAL* H, - /* Input arguments: */ - int Direction, KPP_REAL gam, KPP_REAL Jac0[], - /* Output arguments: */ + /* Inout argument: (step size is decreased when LU fails) */ + KPP_REAL* H, + /* Input arguments: */ + int Direction, KPP_REAL gam, KPP_REAL Jac0[], + /* Output arguments: */ KPP_REAL Ghimj[], int Pivot[] ) /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Prepares the LHS matrix for stage calculations @@ -690,16 +699,16 @@ char ros_PrepareMatrix ( -exit after 5 consecutive fails Return value: Singular (true=1=failed_LU or false=0=successful_LU) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ - /*~~~> Local variables */ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ + /*~~~> Local variables */ int i, ising, Nconsecutive; KPP_REAL ghinv; - + Nconsecutive = 0; - + while (1) { /* while Singular */ - + /*~~~> Construct Ghimj = 1/(H*ham) - Jac0 */ WCOPY(KPP_LU_NONZERO,Jac0,1,Ghimj,1); WSCAL(KPP_LU_NONZERO,(-ONE),Ghimj,1); @@ -722,29 +731,29 @@ char ros_PrepareMatrix ( return 1; /* Singular = true */ } /* end if Nconsecutive */ } /* end if ising */ - + } /* while Singular */ } /* ros_PrepareMatrix */ -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ int ros_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H) /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Handles all error messages and returns IERR = error Code -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ - printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); - printf("\nForced exit from Rosenbrock due to the following error:\n"); - +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ + printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); + printf("\nForced exit from Rosenbrock due to the following error:\n"); + switch (Code) { - case -1: + case -1: printf("--> Improper value for maximal no of steps"); break; - case -2: + case -2: printf("--> Selected Rosenbrock method not implemented"); break; - case -3: + case -3: printf("--> Hmin/Hmax/Hstart must be positive"); break; - case -4: + case -4: printf("--> FacMin/FacMax/FacRej must be positive"); break; case -5: printf("--> Improper tolerance values"); break; @@ -752,90 +761,90 @@ int ros_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H) printf("--> No of steps exceeds maximum bound"); break; case -7: printf("--> Step size too small (T + H/10 = T) or H < Roundoff"); break; - case -8: + case -8: printf("--> Matrix is repeatedly singular"); break; default: - printf("Unknown Error code: %d ",Code); + printf("Unknown Error code: %d ",Code); } /* end switch */ - + printf("\n Time = %15.7e, H = %15.7e",T,H); - printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"); - - return Code; - + printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"); + + return Code; + } /* ros_ErrorMsg */ - -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -void Ros2 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +void Ros2 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ AN L-STABLE METHOD, 2 stages, order 2 -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ double g = (KPP_REAL)1.70710678118655; /* 1.0 + 1.0/SQRT(2.0) */ - + /*~~~> Name of the method */ - strcpy(ros_Name, "ROS-2"); - + strcpy(ros_Name, "ROS-2"); + /*~~~> Number of stages */ *ros_S = 2; - + /*~~~> The coefficient matrices A and C are strictly lower triangular. The lower triangular (subdiagonal) elements are stored in row-wise order: A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. The general mapping formula is: A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ ros_A[0] = 1.0/g; - + /*~~~> C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1] */ ros_C[0] = (-2.0)/g; - + /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ ros_NewF[0] = 1; ros_NewF[1] = 1; - + /*~~~> M_i = Coefficients for new step solution */ ros_M[0]= (3.0)/(2.0*g); ros_M[1]= (1.0)/(2.0*g); - - /*~~~> E_i = Coefficients for error estimator */ + + /*~~~> E_i = Coefficients for error estimator */ ros_E[0] = 1.0/(2.0*g); ros_E[1] = 1.0/(2.0*g); - + /*~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus one */ - *ros_ELO = (KPP_REAL)2.0; - + *ros_ELO = (KPP_REAL)2.0; + /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ ros_Alpha[0] = (KPP_REAL)0.0; - ros_Alpha[1] = (KPP_REAL)1.0; - - /*~~~> Gamma_i = \sum_j gamma_{i,j} */ + ros_Alpha[1] = (KPP_REAL)1.0; + + /*~~~> Gamma_i = \sum_j gamma_{i,j} */ ros_Gamma[0] = g; ros_Gamma[1] = -g; - + } /* Ros2 */ -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -void Ros3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +void Ros3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ /*~~~> Name of the method */ - strcpy(ros_Name, "ROS-3"); + strcpy(ros_Name, "ROS-3"); /*~~~> Number of stages */ *ros_S = 3; - + /*~~~> The coefficient matrices A and C are strictly lower triangular. The lower triangular (subdiagonal) elements are stored in row-wise order: A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. @@ -849,63 +858,63 @@ void Ros3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], ros_C[0] = (KPP_REAL)(-1.0156171083877702091975600115545); ros_C[1] = (KPP_REAL)4.0759956452537699824805835358067; ros_C[2] = (KPP_REAL)9.2076794298330791242156818474003; - + /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ ros_NewF[0] = 1; ros_NewF[1] = 1; ros_NewF[2] = 0; - + /*~~~> M_i = Coefficients for new step solution */ ros_M[0] = (KPP_REAL)1.0; ros_M[1] = (KPP_REAL)6.1697947043828245592553615689730; ros_M[2] = (KPP_REAL)(-0.4277225654321857332623837380651); - - /*~~~> E_i = Coefficients for error estimator */ + + /*~~~> E_i = Coefficients for error estimator */ ros_E[0] = (KPP_REAL)0.5; ros_E[1] = (KPP_REAL)(-2.9079558716805469821718236208017); ros_E[2] = (KPP_REAL)0.2235406989781156962736090927619; - + /*~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus 1 */ - *ros_ELO = (KPP_REAL)3.0; - + *ros_ELO = (KPP_REAL)3.0; + /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ ros_Alpha[0]= (KPP_REAL)0.0; ros_Alpha[1]= (KPP_REAL)0.43586652150845899941601945119356; ros_Alpha[2]= (KPP_REAL)0.43586652150845899941601945119356; - - /*~~~> Gamma_i = \sum_j gamma_{i,j} */ + + /*~~~> Gamma_i = \sum_j gamma_{i,j} */ ros_Gamma[0]= (KPP_REAL)0.43586652150845899941601945119356; ros_Gamma[1]= (KPP_REAL)0.24291996454816804366592249683314; ros_Gamma[2]= (KPP_REAL)2.1851380027664058511513169485832; } /* Ros3 */ -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -void Ros4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +void Ros4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES - L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3 + L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3 E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, - SPRINGER-VERLAG (1990) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ + SPRINGER-VERLAG (1990) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ /*~~~> Name of the method */ - strcpy(ros_Name, "ROS-4"); - + strcpy(ros_Name, "ROS-4"); + /*~~~> Number of stages */ *ros_S = 4; - + /*~~~> The coefficient matrices A and C are strictly lower triangular. The lower triangular (subdiagonal) elements are stored in row-wise order: A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. @@ -925,64 +934,64 @@ void Ros4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], ros_C[3] = (KPP_REAL)(-0.2137148994382534e+01); ros_C[4] = (KPP_REAL)(-0.3214669691237626); ros_C[5] = (KPP_REAL)(-0.6949742501781779); - + /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ ros_NewF[0] = 1; ros_NewF[1] = 1; ros_NewF[2] = 1; ros_NewF[3] = 0; - + /*~~~> M_i = Coefficients for new step solution */ ros_M[0] = (KPP_REAL)0.2255570073418735e+01; ros_M[1] = (KPP_REAL)0.2870493262186792; ros_M[2] = (KPP_REAL)0.4353179431840180; ros_M[3] = (KPP_REAL)0.1093502252409163e+01; - - /*~~~> E_i = Coefficients for error estimator */ + + /*~~~> E_i = Coefficients for error estimator */ ros_E[0] = (KPP_REAL)(-0.2815431932141155); ros_E[1] = (KPP_REAL)(-0.7276199124938920e-01); ros_E[2] = (KPP_REAL)(-0.1082196201495311); ros_E[3] = (KPP_REAL)(-0.1093502252409163e+01); - + /*~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus 1 */ - *ros_ELO = (KPP_REAL)4.0; - + *ros_ELO = (KPP_REAL)4.0; + /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ ros_Alpha[0] = (KPP_REAL)0.0; ros_Alpha[1] = (KPP_REAL)0.1145640000000000e+01; ros_Alpha[2] = (KPP_REAL)0.6552168638155900; ros_Alpha[3] = (KPP_REAL)ros_Alpha[2]; - - /*~~~> Gamma_i = \sum_j gamma_{i,j} */ + + /*~~~> Gamma_i = \sum_j gamma_{i,j} */ ros_Gamma[0] = (KPP_REAL)( 0.5728200000000000); ros_Gamma[1] = (KPP_REAL)(-0.1769193891319233e+01); ros_Gamma[2] = (KPP_REAL)( 0.7592633437920482); ros_Gamma[3] = (KPP_REAL)(-0.1049021087100450); } /* Ros4 */ - -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -void Rodas3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +void Rodas3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ --- A STIFFLY-STABLE METHOD, 4 stages, order 3 -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ /*~~~> Name of the method */ - strcpy(ros_Name, "RODAS-3"); - + strcpy(ros_Name, "RODAS-3"); + /*~~~> Number of stages */ *ros_S = 4; - + /*~~~> The coefficient matrices A and C are strictly lower triangular. The lower triangular (subdiagonal) elements are stored in row-wise order: A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. The general mapping formula is: - A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ + A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ ros_A[0] = (KPP_REAL)0.0; ros_A[1] = (KPP_REAL)2.0; ros_A[2] = (KPP_REAL)0.0; @@ -995,75 +1004,148 @@ void Rodas3 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], ros_C[1] = (KPP_REAL)1.0; ros_C[2] = (KPP_REAL)(-1.0); ros_C[3] = (KPP_REAL)1.0; - ros_C[4] = (KPP_REAL)(-1.0); - ros_C[5] = (KPP_REAL)(-2.66666666666667); /* -8/3 */ - + ros_C[4] = (KPP_REAL)(-1.0); + ros_C[5] = (KPP_REAL)(-2.66666666666667); /* -8/3 */ + /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ ros_NewF[0] = 1; ros_NewF[1] = 0; ros_NewF[2] = 1; ros_NewF[3] = 1; - + /*~~~> M_i = Coefficients for new step solution */ ros_M[0] = (KPP_REAL)2.0; ros_M[1] = (KPP_REAL)0.0; ros_M[2] = (KPP_REAL)1.0; ros_M[3] = (KPP_REAL)1.0; - - /*~~~> E_i = Coefficients for error estimator */ + + /*~~~> E_i = Coefficients for error estimator */ ros_E[0] = (KPP_REAL)0.0; ros_E[1] = (KPP_REAL)0.0; ros_E[2] = (KPP_REAL)0.0; ros_E[3] = (KPP_REAL)1.0; - + /*~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus 1 */ *ros_ELO = (KPP_REAL)3.0; - + /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ ros_Alpha[0] = (KPP_REAL)0.0; ros_Alpha[1] = (KPP_REAL)0.0; ros_Alpha[2] = (KPP_REAL)1.0; ros_Alpha[3] = (KPP_REAL)1.0; - - /*~~~> Gamma_i = \sum_j gamma_{i,j} */ + + /*~~~> Gamma_i = \sum_j gamma_{i,j} */ ros_Gamma[0] = (KPP_REAL)0.5; ros_Gamma[1] = (KPP_REAL)1.5; ros_Gamma[2] = (KPP_REAL)0.0; ros_Gamma[3] = (KPP_REAL)0.0; - + } /* Rodas3 */ - -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], - KPP_REAL ros_M[], KPP_REAL ros_E[], - KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +void Rodas3_1 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], + char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + --- A STIFFLY-STABLE METHOD, 4 stages, order 3 + --- Updated coefficients by Mike Long (17 Jul 2025) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ + /*~~~> Name of the method */ + strcpy(ros_Name, "RODAS-3.1"); + + /*~~~> Number of stages */ + *ros_S = 4; + + /*~~~> The coefficient matrices A and C are strictly lower triangular. + The lower triangular (subdiagonal) elements are stored in row-wise order: + A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. + The general mapping formula is: + A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ + ros_A[0] = (KPP_REAL)0.000000000000000; + ros_A[1] = (KPP_REAL)0.646601929740551; + ros_A[2] = (KPP_REAL)0.409567801987914; + ros_A[3] = (KPP_REAL)0.646601929740551; + ros_A[4] = (KPP_REAL)0.409567801987914; + ros_A[5] = (KPP_REAL)1.000000000000000; + + /*~~~> C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1] */ + ros_C[0] = (KPP_REAL)4.198495621784201; + ros_C[1] = (KPP_REAL)3.711590161613010; + ros_C[2] = (KPP_REAL)(-1.787771994729384); + ros_C[3] = (KPP_REAL)4.458898153216104; + ros_C[4] = (KPP_REAL)(-2.024095448516552); + ros_C[5] = (KPP_REAL)(-2.626700600119396); + + /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) + or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ + ros_NewF[0] = 1; + ros_NewF[1] = 0; + ros_NewF[2] = 1; + ros_NewF[3] = 1; + + /*~~~> M_i = Coefficients for new step solution */ + ros_M[0] = (KPP_REAL)0.646601929740551; + ros_M[1] = (KPP_REAL)0.409567801987914; + ros_M[2] = (KPP_REAL)1.000000000000000; + ros_M[3] = (KPP_REAL)1.000000000000000; + + /*~~~> E_i = Coefficients for error estimator */ + ros_E[0] = (KPP_REAL)0.000000000000000; + ros_E[1] = (KPP_REAL)0.000000000000000; + ros_E[2] = (KPP_REAL)0.000000000000000; + ros_E[3] = (KPP_REAL)1.000000000000000; + + /*~~~> ros_ELO = estimator of local order - the minimum between the +! main and the embedded scheme orders plus 1 */ + *ros_ELO = (KPP_REAL)3.0000000000000000; + + /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ + ros_Alpha[0] = (KPP_REAL)0.000000000000000; + ros_Alpha[1] = (KPP_REAL)0.000000000000000; + ros_Alpha[2] = (KPP_REAL)1.000000000000000; + ros_Alpha[3] = (KPP_REAL)1.000000000000000; + + /*~~~> Gamma_i = \sum_j gamma_{i,j} */ + ros_Gamma[0] = (KPP_REAL)0.515000000000000; + ros_Gamma[1] = (KPP_REAL)1.628546001287715; + ros_Gamma[2] = (KPP_REAL)0.000000000000000; + ros_Gamma[3] = (KPP_REAL)0.000000000000000; + +} /* Rodas3.1 */ + +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], + KPP_REAL ros_M[], KPP_REAL ros_E[], + KPP_REAL ros_Alpha[], KPP_REAL ros_Gamma[], char ros_NewF[], KPP_REAL *ros_ELO, char* ros_Name ) -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, - SPRINGER-VERLAG (1996) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ + SPRINGER-VERLAG (1996) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ { /*~~~> Name of the method */ - strcpy(ros_Name, "RODAS-4"); - + strcpy(ros_Name, "RODAS-4"); + /*~~~> Number of stages */ *ros_S = 6; /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ ros_Alpha[0] = (KPP_REAL)0.000; ros_Alpha[1] = (KPP_REAL)0.386; - ros_Alpha[2] = (KPP_REAL)0.210; + ros_Alpha[2] = (KPP_REAL)0.210; ros_Alpha[3] = (KPP_REAL)0.630; ros_Alpha[4] = (KPP_REAL)1.000; ros_Alpha[5] = (KPP_REAL)1.000; - - /*~~~> Gamma_i = \sum_j gamma_{i,j} */ + + /*~~~> Gamma_i = \sum_j gamma_{i,j} */ ros_Gamma[0] = (KPP_REAL)0.2500000000000000; ros_Gamma[1] = (KPP_REAL)(-0.1043000000000000); ros_Gamma[2] = (KPP_REAL)0.1035000000000000; @@ -1116,7 +1198,7 @@ void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], ros_M[4] = (KPP_REAL)1.0; ros_M[5] = (KPP_REAL)1.0; - /*~~~> E_i = Coefficients for error estimator */ + /*~~~> E_i = Coefficients for error estimator */ ros_E[0] = (KPP_REAL)0.0; ros_E[1] = (KPP_REAL)0.0; ros_E[2] = (KPP_REAL)0.0; @@ -1132,53 +1214,53 @@ void Rodas4 ( int *ros_S, KPP_REAL ros_A[], KPP_REAL ros_C[], ros_NewF[3] = 1; ros_NewF[4] = 1; ros_NewF[5] = 1; - + /*~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus 1 */ *ros_ELO = (KPP_REAL)4.0; - + } /* Rodas4 */ - -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ + +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ void DecompTemplate( KPP_REAL A[], int Pivot[], int* ising ) -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - Template for the LU decomposition -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Template for the LU decomposition +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ *ising = KppDecomp ( A ); /*~~~> Note: for a full matrix use Lapack: DGETRF( KPP_NVAR, KPP_NVAR, A, KPP_NVAR, Pivot, ising ) */ - + Ndec++; } /* DecompTemplate */ - -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ + +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ void SolveTemplate( KPP_REAL A[], int Pivot[], KPP_REAL b[] ) -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - Template for the forward/backward substitution (using pre-computed LU decomposition) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -{ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Template for the forward/backward substitution (using pre-computed LU decomposition) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ KppSolve( A, b ); /*~~~> Note: for a full matrix use Lapack: NRHS = 1 DGETRS( 'N', KPP_NVAR , NRHS, A, KPP_NVAR, Pivot, b, KPP_NVAR, INFO ) */ - + Nsol++; } /* SolveTemplate */ -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] ) -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Template for the ODE function call. - Updates the rate coefficients (and possibly the fixed species) at each call -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ + Updates the rate coefficients (and possibly the fixed species) at each call +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ { - KPP_REAL Told; + KPP_REAL Told; Told = TIME; TIME = T; @@ -1186,31 +1268,29 @@ void FunTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Ydot[] ) Update_RCONST(); Fun( Y, FIX, RCONST, Ydot ); TIME = Told; - + Nfun++; - + } /* FunTemplate */ - -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ + +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Jcb[] ) -/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Template for the ODE Jacobian call. - Updates the rate coefficients (and possibly the fixed species) at each call -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ + Updates the rate coefficients (and possibly the fixed species) at each call +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ { /*~~~> Local variables */ - KPP_REAL Told; + KPP_REAL Told; Told = TIME; - TIME = T ; + TIME = T ; Update_SUN(); Update_RCONST(); Jac_SP( Y, FIX, RCONST, Jcb ); TIME = Told; - - Njac++; - -} /* JacTemplate */ + Njac++; +} /* JacTemplate */ diff --git a/int/rosenbrock.f90 b/int/rosenbrock.f90 index 2460c8d..477aac7 100644 --- a/int/rosenbrock.f90 +++ b/int/rosenbrock.f90 @@ -4,7 +4,9 @@ ! * Ros3 ! ! * Ros4 ! ! * Rodas3 ! +! * Rodas3.1 ! ! * Rodas4 ! +! * Rang3 ! ! By default the code employs the KPP sparse linear algebra routines ! ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! ! ! @@ -192,6 +194,8 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ! = 3 : Ros4 ! = 4 : Rodas3 ! = 5 : Rodas4 +! = 6 : Rang3 +! = 7 : Rodas3.1 ! ! ICNTRL(4) -> maximum number of integration steps ! For ICNTRL(4)=0) the default value of 200000 is used @@ -207,7 +211,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ! = 6 : Call Update_SUN and Update_PHOTO from within the int. ! = 7 : Call Update_SUN, Update_PHOTO, Update_RCONST w/in the int. ! -! ICNTRL(16) -> +! ICNTRL(16) -> ! = 0 : allow negative concentrations (default) ! = 1 : set negative concentrations to zero ! @@ -305,18 +309,20 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & !~~~> Initialize the particular Rosenbrock method selected SELECT CASE (ICNTRL(3)) + CASE (0,4) + CALL Rodas3 CASE (1) CALL Ros2 CASE (2) CALL Ros3 CASE (3) CALL Ros4 - CASE (0,4) - CALL Rodas3 CASE (5) CALL Rodas4 CASE (6) CALL Rang3 + CASE (7) + CALL Rodas3_1 CASE DEFAULT PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) @@ -637,7 +643,7 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & Y = MAX(Ynew,ZERO) ELSE Y(1:N) = Ynew(1:N) - ENDIF + ENDIF T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1109,6 +1115,73 @@ SUBROUTINE Rodas3 END SUBROUTINE Rodas3 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + SUBROUTINE Rodas3_1 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! --- A STIFFLY-STABLE METHOD, 4 stages, order 3 +! --- Updated coefficients by Mike Long (17 Jul 2025) +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + IMPLICIT NONE + + rosMethod = RD3 +!~~~> Name of the method + ros_Name = 'RODAS-3.1' +!~~~> Number of stages + ros_S = 4 + +!~~~> The coefficient matrices A and C are strictly lower triangular. +! The lower triangular (subdiagonal) elements are stored +! in row-wise order: +! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. +! The general mapping formula is: +! A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) +! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) + ros_A(1) = 0.000000000000000_dp + ros_A(2) = 0.646601929740551_dp + ros_A(3) = 0.409567801987914_dp + ros_A(4) = 0.646601929740551_dp + ros_A(5) = 0.409567801987914_dp + ros_A(6) = 1.000000000000000_dp + ros_C(1) = 4.198495621784201_dp + ros_C(2) = 3.711590161613010_dp + ros_C(3) = -1.787771994729384_dp + ros_C(4) = 4.458898153216104_dp + ros_C(5) = -2.024095448516552_dp + ros_C(6) = -2.626700600119396_dp +!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) +! or does it re-use the function evaluation from stage i-1 +! (ros_NewF(i)=FALSE) + ros_NewF(1) = .TRUE. + ros_NewF(2) = .FALSE. + ros_NewF(3) = .TRUE. + ros_NewF(4) = .TRUE. +!~~~> M_i = Coefficients for new step solution + ros_M(1) = 0.646601929740551_dp + ros_M(2) = 0.409567801987914_dp + ros_M(3) = 1.000000000000000_dp + ros_M(4) = 1.000000000000000_dp +!~~~> E_i = Coefficients for error estimator + ros_E(1) = 0.000000000000000_dp + ros_E(2) = 0.000000000000000_dp + ros_E(3) = 0.000000000000000_dp + ros_E(4) = 1.000000000000000_dp +!~~~> ros_ELO = estimator of local order - the minimum between the +! main and the embedded scheme orders plus 1 + ros_ELO = 3.000000000000000_dp +! ~~~> Y_stage_i ~ Y( T + H*Alpha_i ) + ros_Alpha(1) = 0.000000000000000_dp + ros_Alpha(2) = 0.000000000000000_dp + ros_Alpha(3) = 1.000000000000000_dp + ros_Alpha(4) = 1.000000000000000_dp +!~~~> Gamma_i = \sum_j gamma_{i,j} + ros_Gamma(1) = 0.515000000000000_dp + ros_Gamma(2) = 1.628546001287715_dp + ros_Gamma(3) = 0.000000000000000_dp + ros_Gamma(4) = 0.000000000000000_dp + + END SUBROUTINE Rodas3_1 + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Rodas4 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/int/rosenbrock.m b/int/rosenbrock.m index 94ec63d..25e73e6 100644 --- a/int/rosenbrock.m +++ b/int/rosenbrock.m @@ -2,16 +2,17 @@ Rosenbrock(Function, Tspan, Y0, Options, RCNTRL, ICNTRL) %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % Rosenbrock - Implementation of several Rosenbrock methods: -% * Ros2 -% * Ros3 -% * Ros4 -% * Rodas3 -% * Rodas4 -% +% * Ros2 +% * Ros3 +% * Ros4 +% * Rodas3 +% * Rodas4 +% * Rodas3.1 +% % Solves the system y'=F(t,y) using a Rosenbrock method defined by: -% -% G = 1/(H*gamma(1)) - Jac(t0,Y0) -% T_i = t0 + Alpha(i)*H +% +% G = 1/(H*gamma(1)) - Jac(t0,Y0) +% T_i = t0 + Alpha(i)*H % Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j % G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j + % gamma(i)*dF/dT(t0, Y0) @@ -23,9 +24,9 @@ % Springer series in computational mathematics, Springer-Verlag, 1996. % The codes contained in the book inspired this implementation. % -% MATLAB implementation (C) John C. Linford (jlinford@vt.edu). -% Virginia Polytechnic Institute and State University -% November, 2009 +% MATLAB implementation (C) John C. Linford (jlinford@vt.edu). +% Virginia Polytechnic Institute and State University +% November, 2009 % % Based on the Fortran90 implementation (C) Adrian Sandu, August 2004 % and revised by Philipp Miehe and Adrian Sandu, May 2006. @@ -34,15 +35,15 @@ % %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % Input Arguments : -% The first four arguments are similar to the input arguments of +% The first four arguments are similar to the input arguments of % MATLAB's ODE solvers % Function - A function handle for the ODE function % Tspan - The time space to integrate % Y0 - Initial value -% Options - ODE solver options created by odeset(): +% Options - ODE solver options created by odeset(): % AbsTol, InitialStep, Jacobian, MaxStep, and RelTol -% are considered. Other options are ignored. -% 'Jacobian' must be set. +% are considered. Other options are ignored. +% 'Jacobian' must be set. % RCNTRL - real value input parameters (explained below) % ICNTRL - integer input parameters (explained below) %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -76,6 +77,7 @@ % = 3 : Ros4 % = 4 : Rodas3 % = 5 : Rodas4 +% = 7 : Rodas3.1 % % ICNTRL(4) -> maximum number of integration steps % For ICNTRL(4)=0) the default value of 100000 is used @@ -222,6 +224,8 @@ ros_Param = Rodas3; case (5) ros_Param = Rodas4; + case (7) + ros_Param = Rodas3_1; otherwise disp(['Unknown Rosenbrock method: ICNTRL(3)=',num2str(ICNTRL(3))]); IERR = ros_ErrorMsg(-2,Tstart,ZERO); @@ -409,7 +413,7 @@ [H, Ghimj, Pivot, Singular] = ... ros_PrepareMatrix(N, H, Direction, ros_Gamma(1), Jac0); - % Not calculating LU decomposition in ros_PrepareMatrix anymore + % Not calculating LU decomposition in ros_PrepareMatrix anymore % so don't need to check if the matrix is singular % %if Singular % More than 5 consecutive failed decompositions @@ -425,7 +429,7 @@ K(:,1) = K(:,1) + HG * dFdT; end K(:,1) = ros_Solve(Ghimj, Pivot, K(:,1)); - + %~~~> Compute the remaining stages for istage=2:ros_S @@ -615,7 +619,7 @@ %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [X] = ros_Solve(JVS, Pivot, X) %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -% Template for the forward/backward substitution +% Template for the forward/backward substitution %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ global Nsol ISTATUS @@ -875,6 +879,75 @@ return +%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +function [ params ] = Rodas3_1() +%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +% --- A STIFFLY-STABLE METHOD, 4 stages, order 3 +% --- Updated coefficients by Mike Long (17 Jul 2025) +%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +rosMethod = 4; +%~~~> Name of the method +ros_Name = 'RODAS-3.1'; +%~~~> Number of stages +ros_S = 4; + +%~~~> The coefficient matrices A and C are strictly lower triangular. +% The lower triangular (subdiagonal) elements are stored in row-wise order: +% A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. +% The general mapping formula is: +% A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) +% C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) + +ros_A(1) = 0.000000000000000; +ros_A(2) = 0.646601929740551; +ros_A(3) = 0.409567801987914; +ros_A(4) = 0.646601929740551; +ros_A(5) = 0.409567801987914; +ros_A(6) = 1.000000000000000; + +ros_C(1) = 4.198495621784201; +ros_C(2) = 3.711590161613010; +ros_C(3) = -1.787771994729384; +ros_C(4) = 4.458898153216104; +ros_C(5) = -2.024095448516552; +ros_C(6) = -2.626700600119396; + +%~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) +% or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) +ros_NewF(1) = true; +ros_NewF(2) = false; +ros_NewF(3) = true; +ros_NewF(4) = true; +%~~~> M_i = Coefficients for new step solution +ros_M(1) = 0.646601929740551; +ros_M(2) = 0.409567801987914; +ros_M(3) = 1.000000000000000; +ros_M(4) = 1.000000000000000; +%~~~> E_i = Coefficients for error estimator +ros_E(1) = 0.0000000000000000; +ros_E(2) = 0.0000000000000000; +ros_E(3) = 0.0000000000000000; +ros_E(4) = 1.0000000000000000; +%~~~> ros_ELO = estimator of local order - the minimum between the +% main and the embedded scheme orders plus 1 +ros_ELO = 3.0; +%~~~> Y_stage_i ~ Y( T + H*Alpha_i ) +ros_Alpha(1) = 0.0000000000000000; +ros_Alpha(2) = 0.0000000000000000; +ros_Alpha(3) = 1.0000000000000000; +ros_Alpha(4) = 1.0000000000000000; +%~~~> Gamma_i = \sum_j gamma_{i,j} +ros_Gamma(1) = 0.515000000000000; +ros_Gamma(2) = 1.628546001287715; +ros_Gamma(3) = 0.0000000000000000; +ros_Gamma(4) = 0.0000000000000000; + +params = { ros_S, rosMethod, ros_A, ros_C, ros_M, ros_E, ... + ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, ros_Name }; + +return + %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [ params ] = Rodas4() %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -982,5 +1055,3 @@ % End of INTEGRATE function % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - diff --git a/int/rosenbrock_adj.c b/int/rosenbrock_adj.c index eb04fdc..cf0db61 100644 --- a/int/rosenbrock_adj.c +++ b/int/rosenbrock_adj.c @@ -98,6 +98,7 @@ void Ros2(); void Ros3(); void Ros4(); void Rodas3(); +void Rodas3_1(); void Rodas4(); void JacTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Jcb[] ); void HessTemplate( KPP_REAL T, KPP_REAL Y[], KPP_REAL Hes[] ); @@ -395,6 +396,9 @@ int RosenbrockADJ( KPP_REAL Y[], int NADJ, KPP_REAL Lambda[][NVAR], case 5: Rodas4(); break; + case 7: + Rodas3_1(); + break; default: printf( "Unknown Rosenbrock method: ICNTRL[2]=%d", ICNTRL[2] ); return ros_ErrorMsg(-2, Tstart, ZERO); @@ -555,6 +559,8 @@ int RosenbrockADJ( KPP_REAL Y[], int NADJ, KPP_REAL Lambda[][NVAR], break; case 5: Rodas4(); + case 7: + Rodas3_1(); break; default: printf( "Unknown Rosenbrock method: ICNTRL[2]=%d", ICNTRL[2] ); @@ -2293,11 +2299,79 @@ void Rodas3() { } /* End of Rodas3 */ +void Rodas3_1() +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + --- A STIFFLY-STABLE METHOD, 4 stages, order 3 + --- Updated coefficients by Mike Long (22 May 2025) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ +{ + /*~~~> Name of the method */ + strcpy(ros_Name, "RODAS-3.1"); + + /*~~~> Number of stages */ + ros_S = 4; + + /*~~~> The coefficient matrices A and C are strictly lower triangular. + The lower triangular (subdiagonal) elements are stored in row-wise order: + A(2,1) = ros_A[0], A(3,1)=ros_A[1], A(3,2)=ros_A[2], etc. + The general mapping formula is: + A_{i,j} = ros_A[ (i-1)*(i-2)/2 + j -1 ] */ + ros_A[0] = (KPP_REAL)0.00000000000000000; + ros_A[1] = (KPP_REAL)1.5382237953138116; + ros_A[2] = (KPP_REAL)(-0.36440683885434433); + ros_A[3] = (KPP_REAL)1.538223795313811; + ros_A[4] = (KPP_REAL)(-0.36440683885434433); + ros_A[5] = (KPP_REAL)1.0000000000000000; + + /*~~~> C_{i,j} = ros_C[ (i-1)*(i-2)/2 + j -1] */ + ros_C[0] = (KPP_REAL)(-4.0919303685081028); + ros_C[1] = (KPP_REAL)(-3.0551174378039538e-002); + ros_C[2] = (KPP_REAL)1.7259281281917580; + ros_C[3] = (KPP_REAL)0.19561160936073679; + ros_C[4] = (KPP_REAL)1.9301670595355112; + ros_C[5] = (KPP_REAL)(-2.6267006001193960); + + /*~~~> does the stage i require a new function evaluation (ros_NewF(i)=TRUE) + or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) */ + ros_NewF[0] = 1; + ros_NewF[1] = 0; + ros_NewF[2] = 1; + ros_NewF[3] = 1; + + /*~~~> M_i = Coefficients for new step solution */ + ros_M[0] = (KPP_REAL)1.5382237953138116; + ros_M[1] = (KPP_REAL)(-0.36440683885434444); + ros_M[2] = (KPP_REAL)1.0000000000000002; + ros_M[3] = (KPP_REAL)1.0000000000000000; + + /*~~~> E_i = Coefficients for error estimator */ + ros_E[0] = (KPP_REAL)0.0000000000000000; + ros_E[1] = (KPP_REAL)0.0000000000000000; + ros_E[2] = (KPP_REAL)0.0000000000000000; + ros_E[3] = (KPP_REAL)1.0000000000000000; + + /*~~~> ros_ELO = estimator of local order - the minimum between the +! main and the embedded scheme orders plus 1 */ + ros_ELO = (KPP_REAL)3.0000000000000000; + + /*~~~> Y_stage_i ~ Y( T + H*Alpha_i ) */ + ros_Alpha[0] = (KPP_REAL)0.0000000000000000; + ros_Alpha[1] = (KPP_REAL)0.0000000000000000; + ros_Alpha[2] = (KPP_REAL)1.0000000000000000; + ros_Alpha[3] = (KPP_REAL)1.0000000000000000; + + /*~~~> Gamma_i = \sum_j gamma_{i,j} */ + ros_Gamma[0] = (KPP_REAL)0.5150000000000000; + ros_Gamma[1] = (KPP_REAL)(-0.57028223198756145); + ros_Gamma[2] = (KPP_REAL)0.0000000000000000; + ros_Gamma[3] = (KPP_REAL)0.0000000000000000; + +} /* End of Rodas3.1 */ + /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ void Rodas4() { /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES - + STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, diff --git a/int/rosenbrock_adj.f90 b/int/rosenbrock_adj.f90 index aa644f9..753a377 100644 --- a/int/rosenbrock_adj.f90 +++ b/int/rosenbrock_adj.f90 @@ -6,6 +6,7 @@ ! * Ros4 ! ! * Rodas3 ! ! * Rodas4 ! +! * Rodas3.1 ! By default the code employs the KPP sparse linear algebra routines ! ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! ! ! @@ -245,7 +246,8 @@ SUBROUTINE RosenbrockADJ( Y, NADJ, Lambda, & ! = 2 : method is Ros3 ! = 3 : method is Ros4 ! = 4 : method is Rodas3 -! = 5: method is Rodas4 +! = 5 : method is Rodas4 +! = 7 : method is Rodas3.1 ! ! ICNTRL(4) -> maximum number of integration steps ! For ICNTRL(4)=0) the default value of BUFSIZE is used @@ -394,6 +396,8 @@ SUBROUTINE RosenbrockADJ( Y, NADJ, Lambda, & CALL Rodas3 CASE (5) CALL Rodas4 + CASE (7) + CALL Rodas3_1 CASE DEFAULT PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) @@ -555,6 +559,9 @@ SUBROUTINE RosenbrockADJ( Y, NADJ, Lambda, & CALL Rodas3 CASE (5) CALL Rodas4 + CASE (7) + CALL Rodas4 + CASE DEFAULT PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=', ICNTRL(3) CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) @@ -2383,6 +2390,73 @@ SUBROUTINE Rodas3 END SUBROUTINE Rodas3 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + SUBROUTINE Rodas3_1 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! --- A STIFFLY-STABLE METHOD, 4 stages, order 3 +! --- Updated coefficients by Mike Long (17 Jul 2025) +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + IMPLICIT NONE + + rosMethod = RD3 +!~~~> Name of the method + ros_Name = 'RODAS-3.1' +!~~~> Number of stages + ros_S = 4 + +!~~~> The coefficient matrices A and C are strictly lower triangular. +! The lower triangular (subdiagonal) elements are stored +! in row-wise order: +! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. +! The general mapping formula is: +! A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) +! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) + ros_A(1) = 0.000000000000000_dp + ros_A(2) = 0.646601929740551_dp + ros_A(3) = 0.409567801987914_dp + ros_A(4) = 0.646601929740551_dp + ros_A(5) = 0.409567801987914_dp + ros_A(6) = 1.000000000000000_dp + ros_C(1) = 4.198495621784201_dp + ros_C(2) = 3.711590161613010_dp + ros_C(3) = -1.787771994729384_dp + ros_C(4) = 4.458898153216104_dp + ros_C(5) = -2.024095448516552_dp + ros_C(6) = -2.626700600119396_dp +!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) +! or does it re-use the function evaluation from stage i-1 +! (ros_NewF(i)=FALSE) + ros_NewF(1) = .TRUE. + ros_NewF(2) = .FALSE. + ros_NewF(3) = .TRUE. + ros_NewF(4) = .TRUE. +!~~~> M_i = Coefficients for new step solution + ros_M(1) = 0.646601929740551_dp + ros_M(2) = 0.409567801987914_dp + ros_M(3) = 1.000000000000000_dp + ros_M(4) = 1.000000000000000_dp +!~~~> E_i = Coefficients for error estimator + ros_E(1) = 0.000000000000000_dp + ros_E(2) = 0.000000000000000_dp + ros_E(3) = 0.000000000000000_dp + ros_E(4) = 1.000000000000000_dp +!~~~> ros_ELO = estimator of local order - the minimum between the +! main and the embedded scheme orders plus 1 + ros_ELO = 3.000000000000000_dp +! ~~~> Y_stage_i ~ Y( T + H*Alpha_i ) + ros_Alpha(1) = 0.000000000000000_dp + ros_Alpha(2) = 0.000000000000000_dp + ros_Alpha(3) = 1.000000000000000_dp + ros_Alpha(4) = 1.000000000000000_dp +!~~~> Gamma_i = \sum_j gamma_{i,j} + ros_Gamma(1) = 0.515000000000000_dp + ros_Gamma(2) = 1.628546001287715_dp + ros_Gamma(3) = 0.000000000000000_dp + ros_Gamma(4) = 0.000000000000000_dp + + END SUBROUTINE Rodas3_1 + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Rodas4 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/int/rosenbrock_autoreduce.f90 b/int/rosenbrock_autoreduce.f90 index b6e2821..98d197d 100644 --- a/int/rosenbrock_autoreduce.f90 +++ b/int/rosenbrock_autoreduce.f90 @@ -4,7 +4,9 @@ ! * Ros3 ! ! * Ros4 ! ! * Rodas3 ! +! * Rodas3.1 ! ! * Rodas4 ! +! * Rang ! ! By default the code employs the KPP sparse linear algebra routines ! ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! ! ! @@ -197,6 +199,8 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ! = 3 : Ros4 ! = 4 : Rodas3 ! = 5 : Rodas4 +! = 6 : Rang +! = 7 : Rodas3.1 ! ! ICNTRL(4) -> maximum number of integration steps ! For ICNTRL(4)=0) the default value of 200000 is used @@ -217,7 +221,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ! = 6 : Call Update_SUN and Update_PHOTO from within the int. ! = 7 : Call Update_SUN, Update_PHOTO, Update_RCONST w/in the int. ! -! ICNTRL(16) -> +! ICNTRL(16) -> ! = 0 : allow negative concentrations (default) ! = 1 : set negative concentrations to zero ! @@ -322,18 +326,20 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & !~~~> Initialize the particular Rosenbrock method selected SELECT CASE (ICNTRL(3)) + CASE (0,4) + CALL Rodas3 CASE (1) CALL Ros2 CASE (2) CALL Ros3 CASE (3) CALL Ros4 - CASE (0,4) - CALL Rodas3 CASE (5) CALL Rodas4 CASE (6) CALL Rang3 + CASE (7) + CALL Rodas3_1 CASE DEFAULT PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) @@ -484,7 +490,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ! Error indicator IERR) ENDIF - + !~~~> CALL Normal Rosenbrock method IF ( .not. Autoreduce .or. IERR .eq. -99 ) & CALL ros_Integrator(Y, Tstart, Tend, Texit, & @@ -594,7 +600,7 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & DO_SLV = .true. DO_FUN = .true. DO_JVS = .true. - + T = Tstart RSTATUS(Nhexit) = ZERO H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) ) @@ -712,7 +718,7 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & Y = MAX(Ynew,ZERO) ELSE Y(1:N) = Ynew(1:N) - ENDIF + ENDIF T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1181,7 +1187,7 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, & ! 1st order calculation for removed species per Shen et al. (2020) Eq. 4 ! -- currently, DO_FUN() selects ! -- DO_FUN loops over 1,NVAR. Only needs to loop over NVAR-rNVAR - ! but the structure doesn't exist. Maybe worth considering + ! but the structure doesn't exist. Maybe worth considering ! for efficiency purposes. DO i=1,N IF (.not. DO_SLV(i)) THEN @@ -1643,7 +1649,7 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, & ! 1st order calculation for removed species per Shen et al. (2020) Eq. 4 ! -- currently, DO_FUN() selects ! -- DO_FUN loops over 1,NVAR. Only needs to loop over NVAR-rNVAR - ! but the structure doesn't exist. Maybe worth considering + ! but the structure doesn't exist. Maybe worth considering ! for efficiency purposes. DO i=1,N IF (.not. DO_SLV(i)) THEN @@ -1744,7 +1750,7 @@ SUBROUTINE ros_cPrepareMatrix ( H, Direction, gam, & KPP_REAL, INTENT(IN) :: Jac0(N,N) #else KPP_REAL, INTENT(IN) :: Jac0(LU_NONZERO) -#endif +#endif KPP_REAL, INTENT(IN) :: gam INTEGER, INTENT(IN) :: Direction !~~~> Output arguments @@ -1752,7 +1758,7 @@ SUBROUTINE ros_cPrepareMatrix ( H, Direction, gam, & KPP_REAL, INTENT(OUT) :: Ghimj(N,N) #else KPP_REAL, INTENT(OUT) :: Ghimj(cNONZERO) -#endif +#endif LOGICAL, INTENT(OUT) :: Singular INTEGER, INTENT(OUT) :: Pivot(N) !~~~> Inout arguments @@ -1780,7 +1786,7 @@ SUBROUTINE ros_cPrepareMatrix ( H, Direction, gam, & DO i=1,rNVAR Ghimj(cLU_DIAG(i)) = Ghimj(cLU_DIAG(i))+ghinv END DO -#endif +#endif !~~~> Compute LU decomposition CALL ros_cDecomp( Ghimj, Pivot, ISING ) IF (ISING == 0) THEN @@ -1813,7 +1819,7 @@ SUBROUTINE ros_cDecomp( A, Pivot, ISING ) !~~~> Inout variables #ifdef FULL_ALGEBRA KPP_REAL, INTENT(INOUT) :: A(N,N) -#else +#else KPP_REAL, INTENT(INOUT) :: A(cNONZERO) #endif !~~~> Output variables @@ -1821,14 +1827,14 @@ SUBROUTINE ros_cDecomp( A, Pivot, ISING ) #ifdef FULL_ALGEBRA CALL DGETRF( N, N, A, N, Pivot, ISING ) -#else +#else CALL cKppDecomp ( A, ISING ) Pivot(1) = 1 #endif ISTATUS(Ndec) = ISTATUS(Ndec) + 1 END SUBROUTINE ros_cDecomp - + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_cSolve( A, Pivot, b, map1, map2 ) @@ -1840,7 +1846,7 @@ SUBROUTINE ros_cSolve( A, Pivot, b, map1, map2 ) #ifdef FULL_ALGEBRA KPP_REAL, INTENT(IN) :: A(N,N) INTEGER :: ISING -#else +#else KPP_REAL, INTENT(IN) :: A(cNONZERO) #endif INTEGER, INTENT(IN) :: Pivot(N) @@ -1853,8 +1859,8 @@ SUBROUTINE ros_cSolve( A, Pivot, b, map1, map2 ) CALL DGETRS( 'N', N , 1, A, N, Pivot, b, N, ISING ) IF ( Info < 0 ) THEN PRINT*,"Error in DGETRS. ISING=",ISING - END IF -#else + END IF +#else Atmp = 0.d0 Btmp = 0.d0 @@ -2249,6 +2255,73 @@ SUBROUTINE Rodas3 END SUBROUTINE Rodas3 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + SUBROUTINE Rodas3_1 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! --- A STIFFLY-STABLE METHOD, 4 stages, order 3 +! --- Updated coefficients by Mike Long (17 Jul 2025) +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + IMPLICIT NONE + + rosMethod = RD3 +!~~~> Name of the method + ros_Name = 'RODAS-3.1' +!~~~> Number of stages + ros_S = 4 + +!~~~> The coefficient matrices A and C are strictly lower triangular. +! The lower triangular (subdiagonal) elements are stored +! in row-wise order: +! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. +! The general mapping formula is: +! A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) +! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) + ros_A(1) = 0.000000000000000_dp + ros_A(2) = 0.646601929740551_dp + ros_A(3) = 0.409567801987914_dp + ros_A(4) = 0.646601929740551_dp + ros_A(5) = 0.409567801987914_dp + ros_A(6) = 1.000000000000000_dp + ros_C(1) = 4.198495621784201_dp + ros_C(2) = 3.711590161613010_dp + ros_C(3) = -1.787771994729384_dp + ros_C(4) = 4.458898153216104_dp + ros_C(5) = -2.024095448516552_dp + ros_C(6) = -2.626700600119396_dp +!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) +! or does it re-use the function evaluation from stage i-1 +! (ros_NewF(i)=FALSE) + ros_NewF(1) = .TRUE. + ros_NewF(2) = .FALSE. + ros_NewF(3) = .TRUE. + ros_NewF(4) = .TRUE. +!~~~> M_i = Coefficients for new step solution + ros_M(1) = 0.646601929740551_dp + ros_M(2) = 0.409567801987914_dp + ros_M(3) = 1.000000000000000_dp + ros_M(4) = 1.000000000000000_dp +!~~~> E_i = Coefficients for error estimator + ros_E(1) = 0.000000000000000_dp + ros_E(2) = 0.000000000000000_dp + ros_E(3) = 0.000000000000000_dp + ros_E(4) = 1.000000000000000_dp +!~~~> ros_ELO = estimator of local order - the minimum between the +! main and the embedded scheme orders plus 1 + ros_ELO = 3.000000000000000_dp +! ~~~> Y_stage_i ~ Y( T + H*Alpha_i ) + ros_Alpha(1) = 0.000000000000000_dp + ros_Alpha(2) = 0.000000000000000_dp + ros_Alpha(3) = 1.000000000000000_dp + ros_Alpha(4) = 1.000000000000000_dp +!~~~> Gamma_i = \sum_j gamma_{i,j} + ros_Gamma(1) = 0.515000000000000_dp + ros_Gamma(2) = 1.628546001287715_dp + ros_Gamma(3) = 0.000000000000000_dp + ros_Gamma(4) = 0.000000000000000_dp + + END SUBROUTINE Rodas3_1 + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Rodas4 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -2604,7 +2677,7 @@ SUBROUTINE cKppDecomp( JVS, IER ) JVS(kk) = W( cLU_ICOL(kk) ) END DO END DO - + END SUBROUTINE cKppDecomp SUBROUTINE APPEND(IDX) @@ -2616,7 +2689,7 @@ SUBROUTINE APPEND(IDX) DO_SLV(IDX) = .true. DO_FUN(IDX) = .true. ! increment rNVAR - rNVAR = rNVAR+1 + rNVAR = rNVAR+1 ! append SPC_MAP & iSPC_MAP SPC_MAP(rNVAR) = IDX ! From AR to full species iSPC_MAP(IDX) = rNVAR ! From full to AR species diff --git a/int/rosenbrock_h211b_qssa.f90 b/int/rosenbrock_h211b_qssa.f90 index 2d209ca..48ae8fe 100644 --- a/int/rosenbrock_h211b_qssa.f90 +++ b/int/rosenbrock_h211b_qssa.f90 @@ -5,6 +5,8 @@ ! * Ros4 ! ! * Rodas3 ! ! * Rodas4 ! +! * Rang3 ! +! * Rodas3.1 ! ! By default the code employs the KPP sparse linear algebra routines ! ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! ! ! @@ -205,6 +207,8 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ! = 3 : Ros4 ! = 4 : Rodas3 ! = 5 : Rodas4 +! = 6 : Rang3 +! = 7 : Rodas3.1 ! ! ICNTRL(4) -> maximum number of integration steps ! For ICNTRL(4)=0) the default value of 200000 is used @@ -347,6 +351,8 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & CALL Rodas4 CASE (6) CALL Rang3 + CASE (7) + CALL Rodas3_1 CASE DEFAULT PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) @@ -1367,6 +1373,73 @@ SUBROUTINE Rodas3 END SUBROUTINE Rodas3 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + SUBROUTINE Rodas3_1 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! --- A STIFFLY-STABLE METHOD, 4 stages, order 3 +! --- Updated coefficients by Mike Long (17 Jul 2025) +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + IMPLICIT NONE + + rosMethod = RD3 +!~~~> Name of the method + ros_Name = 'RODAS-3.1' +!~~~> Number of stages + ros_S = 4 + +!~~~> The coefficient matrices A and C are strictly lower triangular. +! The lower triangular (subdiagonal) elements are stored +! in row-wise order: +! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. +! The general mapping formula is: +! A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) +! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) + ros_A(1) = 0.000000000000000_dp + ros_A(2) = 0.646601929740551_dp + ros_A(3) = 0.409567801987914_dp + ros_A(4) = 0.646601929740551_dp + ros_A(5) = 0.409567801987914_dp + ros_A(6) = 1.000000000000000_dp + ros_C(1) = 4.198495621784201_dp + ros_C(2) = 3.711590161613010_dp + ros_C(3) = -1.787771994729384_dp + ros_C(4) = 4.458898153216104_dp + ros_C(5) = -2.024095448516552_dp + ros_C(6) = -2.626700600119396_dp +!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) +! or does it re-use the function evaluation from stage i-1 +! (ros_NewF(i)=FALSE) + ros_NewF(1) = .TRUE. + ros_NewF(2) = .FALSE. + ros_NewF(3) = .TRUE. + ros_NewF(4) = .TRUE. +!~~~> M_i = Coefficients for new step solution + ros_M(1) = 0.646601929740551_dp + ros_M(2) = 0.409567801987914_dp + ros_M(3) = 1.000000000000000_dp + ros_M(4) = 1.000000000000000_dp +!~~~> E_i = Coefficients for error estimator + ros_E(1) = 0.000000000000000_dp + ros_E(2) = 0.000000000000000_dp + ros_E(3) = 0.000000000000000_dp + ros_E(4) = 1.000000000000000_dp +!~~~> ros_ELO = estimator of local order - the minimum between the +! main and the embedded scheme orders plus 1 + ros_ELO = 3.000000000000000_dp +! ~~~> Y_stage_i ~ Y( T + H*Alpha_i ) + ros_Alpha(1) = 0.000000000000000_dp + ros_Alpha(2) = 0.000000000000000_dp + ros_Alpha(3) = 1.000000000000000_dp + ros_Alpha(4) = 1.000000000000000_dp +!~~~> Gamma_i = \sum_j gamma_{i,j} + ros_Gamma(1) = 0.515000000000000_dp + ros_Gamma(2) = 1.628546001287715_dp + ros_Gamma(3) = 0.000000000000000_dp + ros_Gamma(4) = 0.000000000000000_dp + + END SUBROUTINE Rodas3_1 + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Rodas4 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/int/rosenbrock_tlm.f90 b/int/rosenbrock_tlm.f90 index de5f3d4..6cebbde 100644 --- a/int/rosenbrock_tlm.f90 +++ b/int/rosenbrock_tlm.f90 @@ -6,6 +6,7 @@ ! * Ros4 ! ! * Rodas3 ! ! * Rodas4 ! +! * Rodas3.1 ! ! By default the code employs the KPP sparse linear algebra routines ! ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! ! ! @@ -233,6 +234,7 @@ SUBROUTINE RosenbrockTLM(N,Y,NTLM,Y_tlm, & ! = 3 : method is Ros4 ! = 4 : method is Rodas3 ! = 5 : method is Rodas4 +! = 7 : method is Rodas3.1 ! ! ICNTRL(4) -> maximum number of integration steps ! For ICNTRL(4)=0) the default value of 200000 is used @@ -348,6 +350,8 @@ SUBROUTINE RosenbrockTLM(N,Y,NTLM,Y_tlm, & CALL Rodas3 CASE (5) CALL Rodas4 + CASE (7) + CALL Rodas3_1 CASE DEFAULT PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) @@ -1225,6 +1229,74 @@ SUBROUTINE Rodas3 END SUBROUTINE Rodas3 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + SUBROUTINE Rodas3_1 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! --- A STIFFLY-STABLE METHOD, 4 stages, order 3 +! --- Updated coefficients by Mike Long (17 Jul 2025) +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + IMPLICIT NONE + + rosMethod = RD3 +!~~~> Name of the method + ros_Name = 'RODAS-3.1' +!~~~> Number of stages + ros_S = 4 + +!~~~> The coefficient matrices A and C are strictly lower triangular. +! The lower triangular (subdiagonal) elements are stored +! in row-wise order: +! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. +! The general mapping formula is: +! A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) +! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) + ros_A(1) = 0.000000000000000_dp + ros_A(2) = 0.646601929740551_dp + ros_A(3) = 0.409567801987914_dp + ros_A(4) = 0.646601929740551_dp + ros_A(5) = 0.409567801987914_dp + ros_A(6) = 1.000000000000000_dp + ros_C(1) = 4.198495621784201_dp + ros_C(2) = 3.711590161613010_dp + ros_C(3) = -1.787771994729384_dp + ros_C(4) = 4.458898153216104_dp + ros_C(5) = -2.024095448516552_dp + ros_C(6) = -2.626700600119396_dp +!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) +! or does it re-use the function evaluation from stage i-1 +! (ros_NewF(i)=FALSE) + ros_NewF(1) = .TRUE. + ros_NewF(2) = .FALSE. + ros_NewF(3) = .TRUE. + ros_NewF(4) = .TRUE. +!~~~> M_i = Coefficients for new step solution + ros_M(1) = 0.646601929740551_dp + ros_M(2) = 0.409567801987914_dp + ros_M(3) = 1.000000000000000_dp + ros_M(4) = 1.000000000000000_dp +!~~~> E_i = Coefficients for error estimator + ros_E(1) = 0.000000000000000_dp + ros_E(2) = 0.000000000000000_dp + ros_E(3) = 0.000000000000000_dp + ros_E(4) = 1.000000000000000_dp +!~~~> ros_ELO = estimator of local order - the minimum between the +! main and the embedded scheme orders plus 1 + ros_ELO = 3.000000000000000_dp +! ~~~> Y_stage_i ~ Y( T + H*Alpha_i ) + ros_Alpha(1) = 0.000000000000000_dp + ros_Alpha(2) = 0.000000000000000_dp + ros_Alpha(3) = 1.000000000000000_dp + ros_Alpha(4) = 1.000000000000000_dp +!~~~> Gamma_i = \sum_j gamma_{i,j} + ros_Gamma(1) = 0.515000000000000_dp + ros_Gamma(2) = 1.628546001287715_dp + ros_Gamma(3) = 0.000000000000000_dp + ros_Gamma(4) = 0.000000000000000_dp + + END SUBROUTINE Rodas3_1 + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Rodas4 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~