Xcos solvers : Sundials Dormand-Price 4(5) implementation 58/9658/59
Clément DAVID [Tue, 2 Oct 2012 15:25:33 +0000 (17:25 +0200)]
Integrated an explicit Dormand-Price 4(5) solver in Sundials
Included tests, doc and corrections to the existing doc

Change-Id: I934a3453f6c22ace789f4ba2fc6a0faf4c86fc09

20 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/include/sundials/sundials_math.h
scilab/modules/scicos/src/scicos_sundials/src/cvode/cvode.c
scilab/modules/scicos/src/scicos_sundials/src/sundials/sundials_math.c
scilab/modules/scicos/tests/unit_tests/DoPri.dia.ref [new file with mode: 0644]
scilab/modules/scicos/tests/unit_tests/DoPri.tst [new file with mode: 0644]
scilab/modules/scicos/tests/unit_tests/LSodar.dia.ref
scilab/modules/scicos/tests/unit_tests/LSodar.tst
scilab/modules/xcos/examples/solvers/integDoPri.sce [new file with mode: 0644]
scilab/modules/xcos/examples/solvers/integLSodar.sce
scilab/modules/xcos/examples/solvers/integRK.sce
scilab/modules/xcos/examples/solvers/integSundials.sce
scilab/modules/xcos/help/en_US/solvers/CVode.xml
scilab/modules/xcos/help/en_US/solvers/Dormand-Price.xml [new file with mode: 0644]
scilab/modules/xcos/help/en_US/solvers/IDA.xml
scilab/modules/xcos/help/en_US/solvers/LSodar.xml
scilab/modules/xcos/help/en_US/solvers/Runge-Kutta.xml
scilab/modules/xcos/tests/unit_tests/DoPri_test.zcos [new file with mode: 0644]

index a05a44b..0f1b08e 100644 (file)
@@ -55,6 +55,8 @@ 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.
+
 * Bug #10040 fixed - VARIABLE_DELAY documentation does not fully describe the
                      initial value behavioral.
 
index f84ef7f..955fa37 100644 (file)
@@ -823,6 +823,10 @@ 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 */
+        {
+            cossim(t0);
+        }
         else if (C2F(cmsolver).solver == 100)  /* IDA  : Method:       , Nonlinear solver=  */
         {
             cossimdaskr(t0);
@@ -1374,6 +1378,9 @@ static void cossim(double *told)
             case 4:
                 cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
                 break;
+            case 5:
+                cvode_mem = CVodeCreate(CV_DOPRI, CV_FUNCTIONAL);
+                break;
         }
 
         /*    cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);*/
index 528fc17..fc044dd 100644 (file)
@@ -92,6 +92,7 @@ extern "C" {
   /* lmm */
 #define CV_ADAMS 1
 #define CV_BDF   2
+#define CV_DOPRI 3
 
   /* iter */
 #define CV_FUNCTIONAL 1
index a887954..c0fdb3f 100644 (file)
@@ -240,10 +240,12 @@ static int CVEwtSetSS(CVodeMem cv_mem, N_Vector ycur, N_Vector weight);
 static int CVEwtSetSV(CVodeMem cv_mem, N_Vector ycur, N_Vector weight);
 
 static int CVHin(CVodeMem cv_mem, realtype tout);
+static int CVHinFixed(CVodeMem cv_mem, realtype tout, realtype *tret);
 static realtype CVUpperBoundH0(CVodeMem cv_mem, realtype tdist);
 static int CVYddNorm(CVodeMem cv_mem, realtype hg, realtype *yddnrm);
 
 static int CVStep(CVodeMem cv_mem);
+static int CVStepDoPri(CVodeMem cv_mem);
 
 static int CVsldet(CVodeMem cv_mem);
 
@@ -333,7 +335,7 @@ void *CVodeCreate(int lmm, int iter)
 
   /* Test inputs */
 
-  if ((lmm != CV_ADAMS) && (lmm != CV_BDF)) {
+  if ((lmm != CV_ADAMS) && (lmm != CV_BDF) && (lmm != CV_DOPRI)) { /* Integration mode : ADAMS, BDF or DOPRI */
     CVProcessError(NULL, 0, "CVODE", "CVodeCreate", MSGCV_BAD_LMM);
     return(NULL);
   }
@@ -352,6 +354,9 @@ 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]) */
+  maxord = (lmm == CV_DOPRI) ? 8 : maxord;
+
   /* copy input parameters into cv_mem */
   cv_mem->cv_lmm  = lmm;
   cv_mem->cv_iter = iter;
@@ -1284,8 +1289,10 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
       return(CV_FIRST_RHSFUNC_ERR);
     }
 
-    /* Set initial h (from H0 or CVHin). */
+    /* Set initial h (from H0 or CVHin). 
+       If Dormand-Price is selected, then we choose to set h to hmax (= 1/hmax_inv). */
 
+    if ((lmm == CV_ADAMS) || (lmm == CV_BDF)) {
     h = hin;
     if ( (h != ZERO) && ((tout-tn)*h < ZERO) ) {
       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_BAD_H0);
@@ -1311,6 +1318,11 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     rh = ABS(h)*hmax_inv;
     if (rh > ONE) h /= rh;
     if (ABS(h) < hmin) h *= hmin/ABS(h);
+       }
+
+    else { /* Compute the fixed step size h, and set the max number of steps */
+      mxstep = CVHinFixed(cv_mem, tout, tret);
+    }
 
     /* Check for approach to tstop */
 
@@ -1525,7 +1537,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);
   } /* end stopping tests block */  
 
   /*
@@ -1557,9 +1569,9 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
 
       if (ewtsetOK != 0) {
 
-        if (itol == CV_WF) 
+        if (itol == CV_WF)
           CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_EWT_NOW_FAIL, tn);
-        else 
+        else
           CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_EWT_NOW_BAD, tn);
        
         istate = CV_ILL_INPUT;
@@ -1603,7 +1615,10 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     }
 
     /* Call CVStep to take a step */
-    kflag = CVStep(cv_mem);
+    if (lmm == CV_DOPRI)
+    kflag = CVStepDoPri(cv_mem);
+       else
+       kflag = CVStep(cv_mem);
 
     /* Process failed step cases, and exit loop */
     if (kflag != CV_SUCCESS) {
@@ -2176,6 +2191,46 @@ static int CVHin(CVodeMem cv_mem, realtype tout)
 }
 
 /*
+ * CVHinFixed
+ *
+ * This routine computes the fixed step size h.
+ * The objective is to approach hmax (= 1/hmax_inv) with h by trying to split the time interval (t-*told) into hmax-long parts.
+ * - if t-*told is smaller than hmax, then set h = t-*told (one iteration)
+ * - 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.
+ * This procedure returns that number (minus 2 because nstloc starts at 0).
+ */
+
+static int CVHinFixed(CVodeMem cv_mem, realtype tout, realtype *tret)
+{
+  long int n_points;
+  realtype interval_size, hmax, test_div, floor_test;
+
+  hmax = 1./hmax_inv;
+  interval_size = tout-*tret;
+
+  if (interval_size <= hmax) {  /* "Small" interval, h is the size of it */
+    n_points = 2;
+    h = interval_size;
+  }
+  else {
+    test_div = interval_size/hmax;
+    floor_test = FLOOR(test_div);
+    if (test_div-floor_test <= TINY) {  /* t-*told divisible by hmax, cutting the interval into hmax-long parts */
+      n_points = floor_test+1;
+      h = interval_size/(n_points-1);
+    }
+    else {  /* Adding a point and h < hmax to fit the interval */
+      n_points = floor_test+2;
+      h = interval_size/(n_points-1);
+    }
+  }
+  return(n_points-1);
+}
+
+/*
  * CVUpperBoundH0
  *
  * This routine sets an upper bound on abs(h0) based on
@@ -2328,6 +2383,117 @@ static int CVStep(CVodeMem cv_mem)
 }
 
 /*
+ * CVStepDoPri
+ *
+ * This routine performs one internal cvode step using the Dormand-Price 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.
+ * While ftempv is just going be an accumulator, tempv will represent the "correction solution".
+ */
+
+static int CVStepDoPri(CVodeMem cv_mem)
+{
+  int retval;
+
+  /* Constant coefficients */
+  realtype c2, c3, c4, c5, b1, b3, b4, b5, b6, d1, d3, d4, d5, d6, d7;
+  realtype a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a73, a74, a75, a76;
+  a21 =  0.2;                   c2 =  0.2;
+  a31 =  0.075;                 c3 =  0.3;
+  a32 =  0.225;                 c4 =  0.8;
+  a41 =  0.9777777777777778;    c5 =  0.8888888888888888;
+  a42 = -3.7333333333333333;    b1 =  0.0911458333333333;
+  a43 =  3.5555555555555556;    b3 =  0.4492362982929021;
+  a51 =  2.9525986892242036;    b4 =  0.6510416666666667;
+  a52 = -11.595793324188386;    b5 = -0.3223761792452830;
+  a53 =  9.8228928516994361;    b6 =  0.1309523809523810;
+  a54 = -0.2908093278463649;    d1 =  0.0899131944444444;
+  a61 =  2.8462752525252525;    d3 =  0.4534890685834082;
+  a62 = -10.757575757575758;    d4 =  0.6140625;
+  a63 =  8.9064227177434725;    d5 = -0.2715123820754717;
+  a64 =  0.2784090909090909;    d6 =  0.0890476190476190;
+  a65 = -0.2735313036020583;    d7 =  0.025;
+  a71 =  0.0911458333333333;
+  a73 =  0.4492362982929021;
+  a74 =  0.6510416666666667;
+  a75 = -0.3223761792452830;
+  a76 =  0.1309523809523810;
+
+  /* K1, K2 */
+  retval = f(tn, zn[0], zn[2], f_data);                    /* zn[2] = K1 = f(Tn, Yn), */
+  N_VLinearSum_Serial (a21*h, zn[2], 1., zn[0], y);        /* y = K1*a21*h + Yn */
+  retval = f(tn + c2*h, y, zn[3], f_data);                 /* zn[3] = K2 = f(Tn+c2*h, Yn + a21*h*K1), */
+  N_VScale(b1, zn[2], y);                                  /* y = b1*K1 + 0*K2, */
+  N_VScale(d1, zn[2], tempv);                              /* y = d1*K1 + 0*K2, */
+
+  /* K3 */
+  N_VLinearSum_Serial (a31*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a31*h, */
+  N_VLinearSum_Serial (a32*h, zn[3], 1., ftemp, ftemp);    /* ftemp += K2*a32*h, */
+  retval = f(tn + c3*h, ftemp, zn[4], f_data);             /* zn[4] = K3 = f(Tn+c3*h, Yn + a31*h*K1 + a32*h*K2), */
+  N_VLinearSum_Serial (1., y, b3, zn[4], y);               /* y = b1*K1 + 0*K2 + b3*K3, */
+  N_VLinearSum_Serial (1., tempv, d3, zn[4], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3, */
+
+  /* K4 */
+  N_VLinearSum_Serial (a41*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a41*h, */
+  N_VLinearSum_Serial (a42*h, zn[3], 1., ftemp, ftemp);    /* ftemp += K2*a42*h, */
+  N_VLinearSum_Serial (a43*h, zn[4], 1., ftemp, ftemp);    /* ftemp += K3*a43*h, */
+  retval = f(tn + c4*h, ftemp, zn[5], f_data);             /* zn[5] = K4 = f(Tn+c4*h, Yn + a41*h*K1 + a42*h*K2 + a43*h*K3), */
+  N_VLinearSum_Serial (1., y, b4, zn[5], y);               /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4, */
+  N_VLinearSum_Serial (1., tempv, d4, zn[5], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4, */
+
+  /* K5 */
+  N_VLinearSum_Serial (a51*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a51*h, */
+  N_VLinearSum_Serial (a52*h, zn[3], 1., ftemp, ftemp);    /* ftemp += K2*a52*h, */
+  N_VLinearSum_Serial (a53*h, zn[4], 1., ftemp, ftemp);    /* ftemp += K3*a53*h, */
+  N_VLinearSum_Serial (a54*h, zn[5], 1., ftemp, ftemp);    /* ftemp += K4*a54*h, */
+  retval = f(tn + c5*h, ftemp, zn[6], f_data);             /* zn[6] = K5 = f(Tn+c5*h, Yn + a51*h*K1 + a52*h*K2 + a53*h*K3 + a54*h*K4), */
+  N_VLinearSum_Serial (1., y, b5, zn[6], y);               /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4 + b5*K5, */
+  N_VLinearSum_Serial (1., tempv, d5, zn[6], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5, */
+
+  /* K6 */
+  N_VLinearSum_Serial (a61*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a61*h, */
+  N_VLinearSum_Serial (a62*h, zn[3], 1., ftemp, ftemp);    /* ftemp += K2*a62*h, */
+  N_VLinearSum_Serial (a63*h, zn[4], 1., ftemp, ftemp);    /* ftemp += K3*a63*h, */
+  N_VLinearSum_Serial (a64*h, zn[5], 1., ftemp, ftemp);    /* ftemp += K3*a64*h, */
+  N_VLinearSum_Serial (a65*h, zn[6], 1., ftemp, ftemp);    /* ftemp += K3*a65*h, */
+  retval = f(tn + h, ftemp, zn[7], f_data);                /* zn[7] = K6 = f(Tn+h, Yn + a61*h*K1 + a62*h*K2 + a63*h*K3 + a64*h*K4 + a65*h*K5), */
+  N_VLinearSum_Serial (1., y, b6, zn[7], y);               /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4 + b5*K5 + b6*K6, */
+  N_VLinearSum_Serial (1., tempv, d6, zn[7], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5 + d6*K6, */
+
+  /* K7 */
+  N_VLinearSum_Serial (a71*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a71*h, */
+  N_VLinearSum_Serial (a73*h, zn[4], 1., ftemp, ftemp);    /* ftemp += K3*a73*h, */
+  N_VLinearSum_Serial (a74*h, zn[5], 1., ftemp, ftemp);    /* ftemp += K3*a74*h, */
+  N_VLinearSum_Serial (a75*h, zn[6], 1., ftemp, ftemp);    /* ftemp += K3*a75*h, */
+  N_VLinearSum_Serial (a76*h, zn[7], 1., ftemp, ftemp);    /* ftemp += K3*a76*h, */
+  retval = f(tn + h, ftemp, zn[8], f_data);                /* zn[8] = K7 = f(Tn+h, Yn + a71*h*K1 + 0*h*K2 + a73*h*K3 + a74*h*K4 + a75*h*K5 + a76*h*K6), */
+  N_VLinearSum_Serial (1., tempv, d7, zn[8], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5 + d6*K6 + d7*K7, */
+
+  /* Yn+1 */
+  N_VLinearSum_Serial(1., zn[0], h, y, zn[0]);             /* zn[0] = Yn+1 = Yn + h*y */
+
+  /* Check for errors in the evaluations of f thanks to retval */
+  if (retval < 0) {
+    CVProcessError(cv_mem, CV_RHSFUNC_FAIL, "Dormand-Price", "CVStepDoPri", MSGCV_RHSFUNC_FAILED, tn);
+    return(CV_RHSFUNC_FAIL);
+  }
+  if (retval > 0) {
+    CVProcessError(cv_mem, CV_FIRST_RHSFUNC_ERR, "Dormand-Price", "CVStepDoPri", 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
  *
  * This routine is called when a change in step size was decided upon,
index 8bc9d59..fd018ba 100644 (file)
@@ -92,3 +92,17 @@ realtype RExp(realtype x)
   return(expl(x));
 #endif
 }
+
+realtype FLOOR(realtype x)
+{
+#if defined(SUNDIALS_USE_GENERIC_MATH)
+  return(floor((double) x));
+#elif defined(SUNDIALS_DOUBLE_PRECISION)
+  return((floor(x));
+#elif defined(SUNDIALS_SINGLE_PRECISION)
+  return(floor(x));
+#elif defined(SUNDIALS_EXTENDED_PRECISION)
+  return(floor(x));
+#endif
+}
+
diff --git a/scilab/modules/scicos/tests/unit_tests/DoPri.dia.ref b/scilab/modules/scicos/tests/unit_tests/DoPri.dia.ref
new file mode 100644 (file)
index 0000000..5c56fb6
--- /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/DoPri.tst");
+//  or test_run('scicos', 'DoPri', ['no_check_error_output']);
+// Import diagram
+loadScicos();
+loadXcosLibs();
+assert_checktrue(importXcosDiagram("SCI/modules/xcos/tests/unit_tests/DoPri_test.zcos"));
+for i=2:4  // 'max step size' = 5*10^-i, precision
+ // Start by updating the clock block period (sampling)
+ scs_m.objs(8).model.rpar(1) = 5*10^(-i);
+ scs_m.objs(9).model.rpar(1) = 5*10^(-i);
+ // Modify solver and 'max step size' + run DoPri + save results
+ scs_m.props.tol(7) = 5*10^(-i); scs_m.props.tol(6) = 5;           // 'max step size' + solver
+ try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;   // DoPri
+ doprival = 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(doprival-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+2));
+ assert_checktrue(mea <= 10^-(i+2));
+ assert_checktrue(stdeviation <= 10^-(i+2));
+end
diff --git a/scilab/modules/scicos/tests/unit_tests/DoPri.tst b/scilab/modules/scicos/tests/unit_tests/DoPri.tst
new file mode 100644 (file)
index 0000000..1464a19
--- /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/DoPri.tst");
+//  or test_run('scicos', 'DoPri', ['no_check_error_output']);
+
+// Import diagram
+loadScicos();
+loadXcosLibs();
+assert_checktrue(importXcosDiagram("SCI/modules/xcos/tests/unit_tests/DoPri_test.zcos"));
+
+for i=2:4  // 'max step size' = 5*10^-i, precision
+
+ // Start by updating the clock block period (sampling)
+ scs_m.objs(8).model.rpar(1) = 5*10^(-i);
+ scs_m.objs(9).model.rpar(1) = 5*10^(-i);
+
+ // Modify solver and 'max step size' + run DoPri + save results
+ scs_m.props.tol(7) = 5*10^(-i); scs_m.props.tol(6) = 5;           // 'max step size' + solver
+ try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;   // DoPri
+ doprival = 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(doprival-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+2));
+ assert_checktrue(mea <= 10^-(i+2));
+ assert_checktrue(stdeviation <= 10^-(i+2));
+
+end
index 83584fa..56408e9 100644 (file)
@@ -1,6 +1,6 @@
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2012 - Scilab Enterprises
+// Copyright (C) 2012 - Scilab Enterprises - Paul Bignier
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
@@ -12,10 +12,10 @@ loadScicos();
 loadXcosLibs();
 assert_checktrue(importXcosDiagram("SCI/modules/xcos/tests/unit_tests/LSodar_test.zcos"));
 // Set solver to LSodar + run LSodar + save results
-scs_m.props.tol(6) = 0;                                                                                        // Set solver to LSodar
-try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;        // Run LSodar
-lsodarval = res.values;                // Results
-time = res.time;                       // Time
+scs_m.props.tol(6) = 0;                                         // Set solver to LSodar
+try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end; // Run LSodar
+lsodarval = res.values;  // Results
+time = res.time;         // Time
 // Set solver to CVode BDF/Newton + run + save results
 scs_m.props.tol(6) = 1;
 try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
index 39d2539..359960f 100644 (file)
@@ -1,6 +1,6 @@
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2012 - Scilab Enterprises
+// Copyright (C) 2012 - Scilab Enterprises - Paul Bignier
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
diff --git a/scilab/modules/xcos/examples/solvers/integDoPri.sce b/scilab/modules/xcos/examples/solvers/integDoPri.sce
new file mode 100644 (file)
index 0000000..10abf40
--- /dev/null
@@ -0,0 +1,34 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - Scilab Enterprises - Paul Bignier
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING,
+// which you should have received as part of this distribution.
+// The terms are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+
+// Run with exec("SCI/modules/xcos/examples/solvers/integDoPri.sce");
+
+// Import the diagram and augment the ending time
+loadScicos();
+loadXcosLibs();
+importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
+scs_m.props.tf = 30000;
+
+solverName=["BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functional", "Dormand-Price"];
+
+for solver=1:5
+
+ // Select the solver
+ scs_m.props.tol(6) = solver;
+
+ // Set max step size if DoPri
+ if (solver == 5) scs_m.props.tol(7) = 0.01;
+
+ // Start the timer, launch the simulation and display time
+ tic();
+ try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
+ t = toc();
+ disp(t, "Time for " + solverName(solver) + " :");
+
+end
index b2dd54f..cb64e4b 100644 (file)
@@ -1,5 +1,5 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2009 - Paul Bignier
+// Copyright (C) 2012 - Scilab Enterprises - Paul Bignier
 //
 // This file must be used under the terms of the CeCILL.
 // This source file is licensed as described in the file COPYING,
index a8e93f3..52822aa 100644 (file)
@@ -1,10 +1,10 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2009 - Paul Bignier
+// Copyright (C) 2012 - Scilab Enterprises - Paul Bignier
 //
 // This file must be used under the terms of the CeCILL.
-// This source file is licensed as described in the file COPYING, which
-// you should have received as part of this distribution.  The terms
-// are also available at
+// This source file is licensed as described in the file COPYING,
+// which you should have received as part of this distribution.
+// The terms are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 // Run with exec("SCI/modules/xcos/examples/solvers/integRK.sce");
@@ -17,18 +17,18 @@ scs_m.props.tf = 30000;
 
 solverName=["BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functional", "Runge-Kutta"];
 
-for solver=0:4
+for solver=1:5
 
  // Select the solver
- scs_m.props.tol(6) = solver+1;
+ scs_m.props.tol(6) = solver;
 
  // Set max step size if Runge-Kutta
- if ((solver+1) == 5) scs_m.props.tol(7) = 0.01;
+ if (solver == 5) scs_m.props.tol(7) = 0.01;
 
  // Start the timer, launch the simulation and display time
  tic();
  try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
  t = toc();
- disp(t, "Time for " + solverName(solver+1) + " :");
+ disp(t, "Time for " + solverName(solver) + " :");
 
 end
index c238b6c..ea5f34f 100644 (file)
@@ -1,10 +1,10 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2009 - Paul Bignier
+// Copyright (C) 2012 - Scilab Enterprises - Paul Bignier
 //
 // This file must be used under the terms of the CeCILL.
-// This source file is licensed as described in the file COPYING, which
-// you should have received as part of this distribution.  The terms
-// are also available at
+// This source file is licensed as described in the file COPYING,
+// which you should have received as part of this distribution.
+// The terms are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 // Run with exec("SCI/modules/xcos/examples/solvers/integ.sce");
@@ -17,15 +17,15 @@ scs_m.props.tf = 30000;
 
 solverName=["BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functional"];
 
-for solver=0:3
+for solver=1:4
 
  // Select the solver
- scs_m.props.tol(6) = solver+1;
+ scs_m.props.tol(6) = solver;
 
  // Start the timer, launch the simulation and display time
  tic();
  try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
  t = toc();
- disp(t, "Time for " + solverName(solver+1) + " :");
+ disp(t, "Time for " + solverName(solver) + " :");
 
 end
index dac1a61..23561e1 100644 (file)
@@ -262,6 +262,9 @@ Time for Adams / Functional :
                 <link linkend="RK">Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="DoPri">Dormand-Price 4(5)</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
diff --git a/scilab/modules/xcos/help/en_US/solvers/Dormand-Price.xml b/scilab/modules/xcos/help/en_US/solvers/Dormand-Price.xml
new file mode 100644 (file)
index 0000000..99fa747
--- /dev/null
@@ -0,0 +1,249 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) Scilab Enterprises - 2012 - Paul Bignier
+ *
+ * This file must be used under the terms of the CeCILL.
+ * This source file is licensed as described in the file COPYING, which
+ * you should have received as part of this distribution.
+ * The terms are also available at
+ * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ -->
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg"  xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en_US" xml:id="DoPri">
+    <refnamediv>
+        <refname>Dormand-Price 4(5)</refname>
+        <refpurpose>
+            <emphasis>Dormand-Price</emphasis> is a numerical solver providing an efficient explicit method to solve Ordinary Differential Equations (ODEs) Initial Value Problems.
+        </refpurpose>
+    </refnamediv>
+    <refsection>
+        <title>Description</title>
+        <para>
+            Called by <link linkend="xcos">xcos</link>, <emphasis>Dormand-Price</emphasis> is a numerical solver providing an efficient fixed-size step method to solve Initial Value Problems of the form :
+        </para>
+        <para>
+            <latex>
+                \begin{eqnarray}
+                \dot{y} = f(t,y), \hspace{3 mm} y(t_0) = y_0, \hspace{3 mm} y \in R^N
+                \end{eqnarray}
+            </latex>
+        </para>
+        <para>
+            <emphasis>CVode</emphasis> and <emphasis>IDA</emphasis> use variable-size steps for the integration.
+        </para>
+        <para>
+            A drawback of that is the unpredictable computation time. With DoPri, we do not adapt to the complexity of the problem, but we guarantee a stable computation time.
+        </para>
+        <para>
+            As of now, this method is explicit, so it is not concerned with Newton or Functional iterations, and not advised for stiff problems.
+        </para>
+        <para>
+            It is an enhancement of the Euler method, which approximates 
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+            by truncating the Taylor expansion.
+        </para>
+        <para>
+            By convention, to use fixed-size steps, the program first computes a fitting <emphasis>h</emphasis> that approaches the simulation parameter <link linkend="Simulatemenu_Menu_entries">max step size</link>.
+        </para>
+        <para>
+            An important difference of DoPri with the previous methods is that it computes up to the seventh derivative of <emphasis>y</emphasis>, while the others only use linear combinations of <emphasis>y</emphasis> and <emphasis>y'</emphasis>.
+        </para>
+        <para>
+            Here, the next value is determined by the present value 
+            <emphasis>
+                y<subscript>n</subscript>
+            </emphasis>
+            plus the weighted average of six increments, where each increment is the product of the size of the interval, <emphasis>h</emphasis>, and an estimated slope specified by the function <emphasis>f(t,y)</emphasis> :
+            <itemizedlist>
+                <listitem>
+                    <emphasis>k1</emphasis> is the increment based on the slope at the beginning of the interval, using 
+                    <emphasis>
+                        y<subscript>n</subscript>
+                    </emphasis>
+                    (Euler's method),
+                </listitem>
+                <listitem>
+                    <emphasis>k2, k3, k4</emphasis> and <emphasis>k5</emphasis> are the increments based on the slope at respectively <emphasis>0.2, 0.3, 0.8</emphasis> and <emphasis>0.9</emphasis> of the interval, using combinations of each other,
+                </listitem>
+                <listitem>
+                    <emphasis>k6</emphasis> is the increment based on the slope at the end, also using combinations of the other <emphasis>ki</emphasis>.
+                </listitem>
+            </itemizedlist>
+        </para>
+        <para>
+            We can see that with the <emphasis>ki</emphasis>, we progress in the derivatives of 
+            <emphasis>
+                y<subscript>n</subscript>
+            </emphasis>
+            . In the computation of the <emphasis>ki</emphasis>, we deliberately use coefficients that yield an error in 
+            <emphasis>
+                O(h<superscript>5</superscript>)
+            </emphasis>
+            at every step.
+        </para>
+        <para>
+            So the total error is 
+            <emphasis>
+                number of steps * O(h<superscript>5</superscript>)
+            </emphasis>
+            . And since <emphasis>number of steps = interval size / h</emphasis> by definition, the total error is in 
+            <emphasis>
+                O(h<superscript>4</superscript>)
+            </emphasis>
+            .
+        </para>
+        <para>
+            That error analysis baptized the method <emphasis>Dopri 4(5)</emphasis> : 
+            <emphasis>
+                O(h<superscript>5</superscript>)
+            </emphasis>
+            per step, 
+            <emphasis>
+                O(h<superscript>4</superscript>)
+            </emphasis>
+            in total.
+        </para>
+        <para>
+            Althought the solver works fine for <link linkend="Simulatemenu_Menu_entries">max step size</link> up to 
+            <emphasis>
+                10<superscript>-3</superscript>
+            </emphasis>
+            , rounding errors sometimes come into play as it approaches 
+            <emphasis>
+                4*10<superscript>-4</superscript>
+            </emphasis>
+            . Indeed, the interval splitting cannot be done properly and we get capricious results.
+        </para>
+    </refsection>
+    <refsection>
+        <title>Examples</title>
+        <para>
+            <link type="scilab" linkend="scilab.xcos/xcos/examples/solvers/ODE_Example.xcos">
+                <inlinemediaobject>
+                    <imageobject>
+                        <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.xcos" valign="middle"/>
+                    </imageobject>
+                </inlinemediaobject>
+            </link>
+            <scilab:image><![CDATA[
+loadScicos();
+loadXcosLibs();
+importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
+scs_m.props.tol(6) = 5;
+scs_m.props.tol(7) = 10^-2;
+try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+]]></scilab:image>
+        </para>
+        <para>
+            The integral block returns its continuous state, we can evaluate it with DoPri by running the example :
+        </para>
+        <para>
+            <programlisting language="example"><![CDATA[
+      // Import the diagram and set the ending time
+      loadScicos();
+      loadXcosLibs();
+      importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
+      scs_m.props.tf = 10000;
+
+      // Select the DoPri solver and set the precision
+      scs_m.props.tol(6) = 5;
+      scs_m.props.tol(7) = 10^-2;
+
+      // Start the timer, launch the simulation and display time
+      tic();
+      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+      t = toc();
+      disp(t, "Time for DoPri :");
+      ]]></programlisting>
+        </para>
+        <para>
+            The Scilab console displays :
+            <screen><![CDATA[
+Time for DoPri :
+ 31.501
+            ]]></screen>
+        </para>
+        <para>
+            Now, in the following script, we compare the time difference between DoPri and Sundials by running the example with the five solvers in turn :
+            <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/integDoPri.sce">
+                Open the script
+            </link>
+        </para>
+        <para>
+            <screen><![CDATA[
+Time for BDF / Newton :
+ 37.882
+
+Time for BDF / Functional :
+ 32.815
+
+Time for Adams / Newton :
+ 30.108
+
+Time for Adams / Functional :
+ 29.242
+
+Time for DoPri :
+ 31.501
+            ]]></screen>
+        </para>
+        <para>
+            These results show that on a nonstiff problem, for relatively same precision required and forcing the same step size, DoPri's computational overhead is significant. Its error to the solution is althought much smaller than the regular Runge-Kutta 4(5), for a small overhead in time.
+        </para>
+        <para>
+            Variable step-size ODE solvers are not appropriate for deterministic real-time applications because the computational overhead of taking a time step varies over the course of an application.
+        </para>
+    </refsection>
+    <refsection>
+        <title>See Also</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="LSodar">LSodar</link>
+            </member>
+            <member>
+                <link linkend="CVode">CVode</link>
+            </member>
+            <member>
+                <link linkend="IDA">IDA</link>
+            </member>
+            <member>
+                <link linkend="RK">Runge-Kutta 4(5)</link>
+            </member>
+            <member>
+                <link linkend="ode">ode</link>
+            </member>
+            <member>
+                <link linkend="ode_discrete">ode_discrete</link>
+            </member>
+            <member>
+                <link linkend="ode_root">ode_root</link>
+            </member>
+            <member>
+                <link linkend="odedc">odedc</link>
+            </member>
+            <member>
+                <link linkend="impl">impl</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection>
+        <title>Bibliography</title>
+        <para>
+            Journal of Computational and Applied Mathematics, Volume 15, Issue 2, 2 June 1986, Pages 203-211 <ulink url="http://dl.acm.org/citation.cfm?id=9958.9963&amp;coll=DL&amp;dl=GUIDE">Dormand-Price Method</ulink>
+        </para>
+        <para>
+            <ulink url="https://computation.llnl.gov/casc/sundials/documentation/documentation.html">Sundials Documentation</ulink>
+        </para>
+    </refsection>
+    <refsection>
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.4.1</revnumber>
+                <revdescription>Dormand-Price 4(5) solver added</revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+</refentry>
index b814b56..7de1c10 100644 (file)
@@ -170,6 +170,9 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
                 <link linkend="RK">Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="DoPri">Dormand-Price 4(5)</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
index f223379..faa8d01 100644 (file)
@@ -168,6 +168,9 @@ Time for Adams / Newton :
                 <link linkend="RK">Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="DoPri">Dormand-Price 4(5)</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
index 39c2b22..1cfbd61 100644 (file)
@@ -225,6 +225,9 @@ Time for Runge-Kutta :
                 <link linkend="LSodar">LSodar</link>
             </member>
             <member>
+                <link linkend="DoPri">Dormand-Price 4(5)</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
diff --git a/scilab/modules/xcos/tests/unit_tests/DoPri_test.zcos b/scilab/modules/xcos/tests/unit_tests/DoPri_test.zcos
new file mode 100644 (file)
index 0000000..60a3be7
Binary files /dev/null and b/scilab/modules/xcos/tests/unit_tests/DoPri_test.zcos differ