* Add a new ODE solver : LSodar. Automatically switches methods to efficiently solve both stiff and nonstiff problems. Includes a rootfinding feature.
-* Add a new ODE solver : Dormand-Price 4(5). Included in the CVode package, so benefits from the rootfinding feature.
+* Add two new Fixed-size step ODE solvers : Dormand-Price 4(5) then Runge-Kutta 4(5). Included in the CVode package, so benefit from the rootfinding feature.
* Bug #10040 fixed - VARIABLE_DELAY documentation does not fully describe the
initial value behavioral.
{
cossim(t0);
}
- else if (C2F(cmsolver).solver == 5) /* DOPRI: Method: Dormand-Price, Nonlinear solver= FUNCTIONAL */
+ else if (C2F(cmsolver).solver == 5) /* DOPRI: Method: Dormand-Price, Nonlinear solver= */
+ {
+ cossim(t0);
+ }
+ else if (C2F(cmsolver).solver == 6) /* RK45: Method: Runge-Kutta, Nonlinear solver= */
{
cossim(t0);
}
case 5:
cvode_mem = CVodeCreate(CV_DOPRI, CV_FUNCTIONAL);
break;
+ case 6:
+ cvode_mem = CVodeCreate(CV_ExpRK, CV_FUNCTIONAL);
+ break;
}
/* cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);*/
#define CV_ADAMS 1
#define CV_BDF 2
#define CV_DOPRI 3
+#define CV_ExpRK 4
/* iter */
#define CV_FUNCTIONAL 1
static int CVStep(CVodeMem cv_mem);
static int CVStepDoPri(CVodeMem cv_mem);
+static int CVStepExpRK(CVodeMem cv_mem);
static int CVsldet(CVodeMem cv_mem);
/* Test inputs */
- if ((lmm != CV_ADAMS) && (lmm != CV_BDF) && (lmm != CV_DOPRI)) { /* Integration mode : ADAMS, BDF or DOPRI */
+ if ((lmm != CV_ADAMS) && (lmm != CV_BDF) && (lmm != CV_DOPRI) && (lmm != CV_ExpRK)) { /* Integration mode : ADAMS, BDF or RK-based */
CVProcessError(NULL, 0, "CVODE", "CVodeCreate", MSGCV_BAD_LMM);
return(NULL);
}
maxord = (lmm == CV_ADAMS) ? ADAMS_Q_MAX : BDF_Q_MAX;
- /* If Runge-Kutta is selected, then maxord = 9 to use the 7 extra vectors allocated (zn[2, 3, 4, 5, 6, 7, 8]) */
+ /* If Dormand-Price is selected, then maxord = 8 to use the 7 extra vectors allocated (zn[2, 3, 4, 5, 6, 7, 8]) */
maxord = (lmm == CV_DOPRI) ? 8 : maxord;
+ /* If Runge-Kutta is selected, then maxord = 1 */
+ maxord = (lmm == CV_ExpRK) ? 1 : maxord;
+
/* copy input parameters into cv_mem */
cv_mem->cv_lmm = lmm;
cv_mem->cv_iter = iter;
}
/* Set initial h (from H0 or CVHin).
- If Dormand-Price is selected, then we choose to set h to hmax (= 1/hmax_inv). */
+ If an RK-based method is selected, then we choose to set h to hmax (= 1/hmax_inv). */
if ((lmm == CV_ADAMS) || (lmm == CV_BDF)) {
h = hin;
}
} /* end of istop tests block */
- if (lmm==CV_DOPRI) mxstep = CVHinFixed(cv_mem, tout, tret);
+ if ((lmm==CV_DOPRI) || (lmm==CV_ExpRK)) mxstep = CVHinFixed(cv_mem, tout, tret);
} /* end stopping tests block */
/*
}
/* Call CVStep to take a step */
- if (lmm == CV_DOPRI)
- kflag = CVStepDoPri(cv_mem);
- else
+ switch (lmm) {
+ case CV_DOPRI:
+ kflag = CVStepDoPri(cv_mem);
+ break;
+ case CV_ExpRK:
+ kflag = CVStepExpRK(cv_mem);
+ break;
+ default:
kflag = CVStep(cv_mem);
+ }
/* Process failed step cases, and exit loop */
if (kflag != CV_SUCCESS) {
* - if it is divisible by hmax, then set h = hmax.
* - if it is not, then "add an integration point" by setting h < hmax, just enough to fit the interval
*
- * Dormand-Price being a fixed step size method, we know the maximum number of steps to take.
+ * RK-based methods using Fixed-size step, we know the maximum number of steps to take.
* This procedure returns that number (minus 2 because nstloc starts at 0).
*/
/*
* CVStepDoPri
*
- * This routine performs one internal cvode step using the Dormand-Price method, from tn to tn + h.
+ * This routine performs one internal cvode step using the Dormand-Price 4(5) method, from tn to tn + h.
* Proceed to computing the K[i] coefficients, build the final solution, increment tn and return CV_SUCCESS.
*
* In order to temporarily store the results, we use tempv and ftemp.
return(CV_SUCCESS);
}
+
+/*
+ * CVStepExpRK
+ *
+ * This routine performs one internal cvode step using the Runge-Kutta 4(5) method, from tn to tn + h.
+ * Proceed to :
+ * - K1 = F(Tn, Yn),
+ * - K2 = F(Tn + h/2, Yn + (h/2)*K1),
+ * - K3 = F(Tn + h/2, Yn + (h/2)*K2),
+ * - K4 = F(Tn + h, Yn + h*K3),
+ * - Yn+1 = Yn + (h/6)*(K1 + 2K2 + 2K3 + K4)
+ * - increment tn
+ *
+ * In order to temporarily store the results, we use tempv and ftemp, which will represent the Ki in turn.
+ */
+
+static int CVStepExpRK(CVodeMem cv_mem)
+{
+ int retval;
+
+ retval = f(tn, zn[0], ftemp, f_data); /* ftemp = K1, */
+ N_VLinearSum_Serial (h/2, ftemp, 1., zn[0], y); /* y = K1*h/2 + Yn */
+ retval = f(tn + h/2, y, tempv, f_data); /* tempv = K2 = f(Tn+h/2, Yn + (h/2)*K1), */
+ N_VLinearSum_Serial (1., ftemp, 2., tempv, y); /* y = K1 + 2K2 */
+
+ N_VLinearSum_Serial (h/2, tempv, 1., zn[0], ftemp); /* ftemp = Yn + K2*h/2, */
+ retval = f(tn + h/2, ftemp, tempv, f_data); /* tempv = K3 = f(Tn+h/2, Yn + (h/2)*K2), */
+ N_VLinearSum_Serial (1., y, 2., tempv, y); /* y = K1 + 2K2 + 2K3, */
+
+ N_VLinearSum_Serial (h, tempv, 1., zn[0], ftemp); /* ftemp = Yn + K3*h, */
+ retval = f(tn + h, ftemp, tempv, f_data); /* tempv = K4 = f(Tn+h/2, Yn + h*K3), */
+ N_VLinearSum_Serial (1., y, 1., tempv, y); /* y = K1 + 2K2 + 2K3 + K4, */
+
+ N_VLinearSum_Serial(1., zn[0], h/6, y, zn[0]); /* zn[0] = Yn+1 = Yn + y*h/6 */
+
+ /* Check for errors in the evaluations of f thanks to retval */
+ if (retval < 0) {
+ CVProcessError(cv_mem, CV_RHSFUNC_FAIL, "Runge-Kutta", "CVStepExpRK", MSGCV_RHSFUNC_FAILED, tn);
+ return(CV_RHSFUNC_FAIL);
+ }
+ if (retval > 0) {
+ CVProcessError(cv_mem, CV_FIRST_RHSFUNC_ERR, "Runge-Kutta", "CVStepExpRK", MSGCV_RHSFUNC_FIRST);
+ return(CV_FIRST_RHSFUNC_ERR);
+ }
+
+ /* Increment tn => take a step. Increment solver calls as well */
+ tn += h;
+ nst++;
+
+ /* Update the Nordsieck history array */
+ retval = f(tn, zn[0], zn[1], f_data); /* zn[1] = y'(tn) */
+
+ return(CV_SUCCESS);
+}
+
+
/*
* CVAdjustParams
*
--- /dev/null
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - Scilab Enterprises - Paul Bignier
+//
+// This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- ENGLISH IMPOSED -->
+// Execute with exec("SCI/modules/scicos/tests/unit_tests/ExpRK.tst");
+// or test_run('scicos', 'ExpRK', ['no_check_error_output']);
+// Import diagram
+loadScicos();
+loadXcosLibs();
+assert_checktrue(importXcosDiagram("SCI/modules/xcos/tests/unit_tests/RK_test.zcos"));
+for i=2:4 // 'max step size' = 10^-i, precision
+ // Start by updating the clock block period (sampling)
+ scs_m.objs(7).model.rpar(1) = 5*10^(-i);
+ scs_m.objs(8).model.rpar(1) = 5*10^(-i);
+ // Modify solver and 'max step size' + run ExpRK + save results
+ scs_m.props.tol(7) = 5*10^(-i); scs_m.props.tol(6) = 6; // 'max step size' + solver
+ try scicos_simulate(scs_m); catch disp(lasterror()); end; // ExpRK
+ rkval = res.values; // Results
+ time = res.time; // Time
+ // Modify solver and 'max step size' + run CVode + save results
+ scs_m.props.tol(7) = 0; scs_m.props.tol(6) = 1;
+ try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
+ cvval = res.values;
+ // Compare results
+ compa = abs(rkval-cvval);
+ // Extract mean, standard deviation, maximum
+ mea = mean(compa);
+ [maxi, indexMaxi] = max(compa);
+ stdeviation = st_deviation(compa);
+ // Verifying closeness of the results
+ assert_checktrue(maxi <= 10^-(i+1));
+ assert_checktrue(mea <= 10^-(i+1));
+ assert_checktrue(stdeviation <= 10^-(i+2));
+end
--- /dev/null
+
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - Scilab Enterprises - Paul Bignier
+//
+// This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- ENGLISH IMPOSED -->
+
+// Execute with exec("SCI/modules/scicos/tests/unit_tests/ExpRK.tst");
+// or test_run('scicos', 'ExpRK', ['no_check_error_output']);
+
+// Import diagram
+loadScicos();
+loadXcosLibs();
+assert_checktrue(importXcosDiagram("SCI/modules/xcos/tests/unit_tests/RK_test.zcos"));
+
+for i=2:4 // 'max step size' = 10^-i, precision
+
+ // Start by updating the clock block period (sampling)
+ scs_m.objs(7).model.rpar(1) = 5*10^(-i);
+ scs_m.objs(8).model.rpar(1) = 5*10^(-i);
+
+ // Modify solver and 'max step size' + run ExpRK + save results
+ scs_m.props.tol(7) = 5*10^(-i); scs_m.props.tol(6) = 6; // 'max step size' + solver
+ try scicos_simulate(scs_m); catch disp(lasterror()); end; // ExpRK
+ rkval = res.values; // Results
+ time = res.time; // Time
+
+ // Modify solver and 'max step size' + run CVode + save results
+ scs_m.props.tol(7) = 0; scs_m.props.tol(6) = 1;
+ try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
+ cvval = res.values;
+
+ // Compare results
+ compa = abs(rkval-cvval);
+
+ // Extract mean, standard deviation, maximum
+ mea = mean(compa);
+ [maxi, indexMaxi] = max(compa);
+ stdeviation = st_deviation(compa);
+
+ // Verifying closeness of the results
+ assert_checktrue(maxi <= 10^-(i+1));
+ assert_checktrue(mea <= 10^-(i+1));
+ assert_checktrue(stdeviation <= 10^-(i+2));
+
+end
<solver code="2" description="Sundials/CVODE - BDF - FUNCTIONAL"/>
<solver code="3" description="Sundials/CVODE - ADAMS - NEWTON"/>
<solver code="4" description="Sundials/CVODE - ADAMS - FUNCTIONAL"/>
- <solver code="5" description="DOPRI5 - Dormand-Price 4(5)"/>
+ <solver code="5" description="DOPRI5 - Dormand-Prince 4(5)"/>
+ <solver code="6" description="RK45 - Runge-Kutta 4(5)"/>
<solver code="100" description="Sundials/IDA"/>
<!-- Available trace values lists from scicos -->
<trace code="0" description="_(No trace nor debug printing)"/>
for solver=1:5
- // Select the solver
+ // Select the solver (Runge-Kutta is solver number 6)
scs_m.props.tol(6) = solver;
+ if (solver == 5) then scs_m.props.tol(6) = 6; end
- // Set max step size if Runge-Kutta
+ // Set max step size if Runge-Kutta is selected
if (solver == 5) then scs_m.props.tol(7) = 0.01; end
// Start the timer, launch the simulation and display time
loadScicos();
loadXcosLibs();
importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
-scs_m.props.tol(6) = 5;
+scs_m.props.tol(6) = 6;
scs_m.props.tol(7) = 10^-2;
try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
]]></scilab:image>
scs_m.props.tf = 10000;
// Select the solver Runge-Kutta and set the precision
- scs_m.props.tol(6) = 5;
+ scs_m.props.tol(6) = 6;
scs_m.props.tol(7) = 10^-2;
// Start the timer, launch the simulation and display time
JLabel solverLabel = new JLabel(XcosMessages.SOLVER_CHOICE);
final String[] solvers = new String[] { "LSodar", "Sundials/CVODE - BDF - NEWTON", "Sundials/CVODE - BDF - FUNCTIONAL",
- "Sundials/CVODE - ADAMS - NEWTON", "Sundials/CVODE - ADAMS - FUNCTIONAL", "DOPRI5 - Dormand-Price 4(5)", "Sundials/IDA"
+ "Sundials/CVODE - ADAMS - NEWTON", "Sundials/CVODE - ADAMS - FUNCTIONAL", "DOPRI5 - Dormand-Prince 4(5)", "RK45 - Runge-Kutta 4(5)", "Sundials/IDA"
};
final String[] solversTooltips = new String[] { "Method: dynamic, Nonlinear solver= dynamic", "Method: BDF, Nonlinear solver= NEWTON",
"Method: BDF, Nonlinear solver= FUNCTIONAL", "Method: ADAMS, Nonlinear solver= NEWTON", "Method: ADAMS, Nonlinear solver= FUNCTIONAL",
- "Method: Fixed step", "Sundials/IDA"
+ "Method: Fixed step", "Method: Fixed step", "Sundials/IDA"
};
solver = new JComboBox(solvers);
/**
* <ul>
- * <li>0 : LSodar
+ * <li>0 : LSodar : Method: DYNAMIC, Nonlinear solver= DYNAMIC
* <li>1 : Sundials/CVODE : Method: BDF, Nonlinear solver= FUNCTIONAL
* <li>2 : Sundials/CVODE : Method: BDF, Nonlinear solver= FUNCTIONAL
* <li>3 : Sundials/CVODE : Method: ADAMS, Nonlinear solver= NEWTON
* <li>4 : Sundials/CVODE : Method: ADAMS, Nonlinear solver= FUNCTIONAL
- * <li>5 : DOPRI5 : Method: Runge-Kutta 4(5)
+ * <li>5 : DOPRI5 : Method: Dormand-Prince 4(5)
+ * <li>6 : RK45 : Method: Runge-Kutta 4(5)
* <li>100 : Sundials/IDA
*
*
/**
* <ul>
- * <li>0 : Lsodar : Method: BDF, Nonlinear solver= NEWTON
+ * <li>0 : LSodar : Method: DYNAMIC, Nonlinear solver= DYNAMIC
* <li>1 : Sundials/CVODE : Method: BDF, Nonlinear solver= FUNCTIONAL
* <li>2 : Sundials/CVODE : Method: BDF, Nonlinear solver= FUNCTIONAL
* <li>3 : Sundials/CVODE : Method: ADAMS, Nonlinear solver= NEWTON
* <li>4 : Sundials/CVODE : Method: ADAMS, Nonlinear solver= FUNCTIONAL
- * <li>5 : DOPRI5 : Method: Runge-Kutta 4(5)
+ * <li>5 : DOPRI5 : Method: Dormand-Prince 4(5)
+ * <li>6 : RK45 : Method: Runge-Kutta 4(5)
* <li>100 : Sundials/IDA
*
*
projects / scilab.git / commitdiff