Xcos solvers : Sundials explicit Runge-Kutta 4(5) implem and interface 31/10031/8
Paul BIGNIER [Thu, 13 Dec 2012 09:25:49 +0000 (10:25 +0100)]
Integrated an explicit Runge-Kutta 4(5) solver in Sundials
to add a fixed-size step method
Set it available in the Xcos interface and the doc

Change-Id: I934a3451f6c22ace789f4ba2fc6a0faf4c86fc29

12 files changed:
scilab/CHANGES_5.4.X
scilab/modules/scicos/src/c/scicos.c
scilab/modules/scicos/src/scicos_sundials/include/cvode/cvode.h
scilab/modules/scicos/src/scicos_sundials/src/cvode/cvode.c
scilab/modules/scicos/tests/unit_tests/ExpRK.dia.ref [new file with mode: 0644]
scilab/modules/scicos/tests/unit_tests/ExpRK.tst [new file with mode: 0644]
scilab/modules/xcos/etc/XConfiguration-xcos.xml
scilab/modules/xcos/examples/solvers/integRK.sce
scilab/modules/xcos/help/en_US/solvers/2-Runge-Kutta.xml
scilab/modules/xcos/src/java/org/scilab/modules/xcos/actions/dialog/SetupDialog.java
scilab/modules/xcos/src/java/org/scilab/modules/xcos/graph/ScicosParameters.java
scilab/modules/xcos/tests/unit_tests/RK_test.zcos [new file with mode: 0644]

index b1e0182..2e56ea9 100644 (file)
@@ -57,7 +57,7 @@ Xcos
 
 * 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.
index 955fa37..9bc5851 100644 (file)
@@ -823,7 +823,11 @@ int C2F(scicos)(double *x_in, int *xptr_in, double *z__,
         {
             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);
         }
@@ -1381,6 +1385,9 @@ static void cossim(double *told)
             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);*/
index fc044dd..0437717 100644 (file)
@@ -93,6 +93,7 @@ extern "C" {
 #define CV_ADAMS 1
 #define CV_BDF   2
 #define CV_DOPRI 3
+#define CV_ExpRK 4
 
   /* iter */
 #define CV_FUNCTIONAL 1
index 172e26c..847493a 100644 (file)
@@ -247,6 +247,7 @@ static int CVYddNorm(CVodeMem cv_mem, realtype hg, realtype *yddnrm);
 
 static int CVStep(CVodeMem cv_mem);
 static int CVStepDoPri(CVodeMem cv_mem);
+static int CVStepExpRK(CVodeMem cv_mem);
 
 static int CVsldet(CVodeMem cv_mem);
 
@@ -336,7 +337,7 @@ void *CVodeCreate(int lmm, int iter)
 
   /* 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);
   }
@@ -355,9 +356,12 @@ void *CVodeCreate(int lmm, int iter)
 
   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;
@@ -1291,7 +1295,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     }
 
     /* 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;
@@ -1538,7 +1542,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
       }
 
     } /* 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 */  
 
   /*
@@ -1616,10 +1620,16 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     }
 
     /* 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) {
@@ -2200,7 +2210,7 @@ static int CVHin(CVodeMem cv_mem, realtype tout)
  * - 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).
  */
 
@@ -2386,7 +2396,7 @@ static int CVStep(CVodeMem cv_mem)
 /*
  * 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.
@@ -2494,6 +2504,62 @@ static int CVStepDoPri(CVodeMem cv_mem)
   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
  *
diff --git a/scilab/modules/scicos/tests/unit_tests/ExpRK.dia.ref b/scilab/modules/scicos/tests/unit_tests/ExpRK.dia.ref
new file mode 100644 (file)
index 0000000..acdbd7f
--- /dev/null
@@ -0,0 +1,37 @@
+// =============================================================================
+// 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
diff --git a/scilab/modules/scicos/tests/unit_tests/ExpRK.tst b/scilab/modules/scicos/tests/unit_tests/ExpRK.tst
new file mode 100644 (file)
index 0000000..fb934d1
--- /dev/null
@@ -0,0 +1,49 @@
+
+// =============================================================================
+// 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
index e44d31a..1d6d58b 100644 (file)
@@ -24,7 +24,8 @@
                 <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)"/>
index 885ba60..1167a04 100644 (file)
@@ -19,10 +19,11 @@ solverName=["BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functional",
 
 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
index ca278f6..28891ff 100644 (file)
 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>
@@ -164,7 +164,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
       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
index 702d43b..a72ee25 100644 (file)
@@ -175,11 +175,11 @@ public class SetupDialog extends JDialog {
 
         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);
index 16e9740..df65868 100644 (file)
@@ -297,12 +297,13 @@ public class ScicosParameters implements Serializable, Cloneable {
 
     /**
      * <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
      *
      *
@@ -314,12 +315,13 @@ public class ScicosParameters implements Serializable, Cloneable {
 
     /**
      * <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
      *
      *
diff --git a/scilab/modules/xcos/tests/unit_tests/RK_test.zcos b/scilab/modules/xcos/tests/unit_tests/RK_test.zcos
new file mode 100644 (file)
index 0000000..36c3a80
Binary files /dev/null and b/scilab/modules/xcos/tests/unit_tests/RK_test.zcos differ