Xcos solvers : Sundials implicit Runge-Kutta 4(5) implem and interface 38/10038/13
Paul BIGNIER [Mon, 17 Dec 2012 15:51:36 +0000 (16:51 +0100)]
Integrated an implicit Runge-Kutta 4(5) solver in Sundials
to add a fixed-size step method
Set it available in the Xcos interface and added the doc

Change-Id: I934a3452f6c22ace789f4ba2fc6a0faf4c86fc19

18 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/ImpRK.dia.ref [new file with mode: 0644]
scilab/modules/scicos/tests/unit_tests/ImpRK.tst [new file with mode: 0644]
scilab/modules/xcos/etc/XConfiguration-xcos.xml
scilab/modules/xcos/examples/solvers/ODE_Example.xcos
scilab/modules/xcos/examples/solvers/integImpRK.sce [new file with mode: 0644]
scilab/modules/xcos/examples/solvers/integRK.sce
scilab/modules/xcos/help/en_US/solvers/0-LSodar.xml
scilab/modules/xcos/help/en_US/solvers/1-CVode.xml
scilab/modules/xcos/help/en_US/solvers/2-Runge-Kutta.xml
scilab/modules/xcos/help/en_US/solvers/3-Dormand-Price.xml
scilab/modules/xcos/help/en_US/solvers/4-ImplicitRK.xml [new file with mode: 0644]
scilab/modules/xcos/help/en_US/solvers/7-IDA.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

index cfa007b..a660e1d 100644 (file)
@@ -59,6 +59,8 @@ Xcos
 
 * 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.
 
+* Add an implicit Fixed-size stop ODE solver : Implicit Runge-Kutta 4(5). Also benefits from the CVode rootfinding feature.
+
 * Bug #10040 fixed - VARIABLE_DELAY documentation does not fully describe the
                      initial value behavioral.
 
index 9bc5851..c65bbb8 100644 (file)
@@ -831,6 +831,10 @@ int C2F(scicos)(double *x_in, int *xptr_in, double *z__,
         {
             cossim(t0);
         }
+        else if (C2F(cmsolver).solver == 7)   /*  ImpRK45: Method: Runge-Kutta, Nonlinear solver= FIXED-POINT */
+        {
+            cossim(t0);
+        }
         else if (C2F(cmsolver).solver == 100)  /* IDA  : Method:       , Nonlinear solver=  */
         {
             cossimdaskr(t0);
@@ -1388,6 +1392,9 @@ static void cossim(double *told)
             case 6:
                 cvode_mem = CVodeCreate(CV_ExpRK, CV_FUNCTIONAL);
                 break;
+            case 7:
+                cvode_mem = CVodeCreate(CV_ImpRK, CV_FUNCTIONAL);
+                break;
         }
 
         /*    cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);*/
index 0437717..0d41535 100644 (file)
@@ -94,6 +94,7 @@ extern "C" {
 #define CV_BDF   2
 #define CV_DOPRI 3
 #define CV_ExpRK 4
+#define CV_ImpRK 5
 
   /* iter */
 #define CV_FUNCTIONAL 1
index 847493a..1ccc3dd 100644 (file)
@@ -248,6 +248,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 CVStepImpRK(CVodeMem cv_mem);
 
 static int CVsldet(CVodeMem cv_mem);
 
@@ -337,7 +338,7 @@ void *CVodeCreate(int lmm, int iter)
 
   /* Test inputs */
 
-  if ((lmm != CV_ADAMS) && (lmm != CV_BDF) && (lmm != CV_DOPRI) && (lmm != CV_ExpRK)) { /* Integration mode : ADAMS, BDF or RK-based */
+  if ((lmm != CV_ADAMS) && (lmm != CV_BDF) && (lmm != CV_DOPRI) && (lmm != CV_ExpRK) && (lmm != CV_ImpRK)) { /* Integration mode : ADAMS, BDF or RK-based */
     CVProcessError(NULL, 0, "CVODE", "CVodeCreate", MSGCV_BAD_LMM);
     return(NULL);
   }
@@ -362,6 +363,9 @@ void *CVodeCreate(int lmm, int iter)
   /* If Runge-Kutta is selected, then maxord = 1 */
   maxord = (lmm == CV_ExpRK) ? 1 : maxord;
 
+  /* If implicit Runge-Kutta is selected, then maxord = 4 to use the 3 extra vectors allocated (zn[2, 3, 4]) */
+  maxord = (lmm == CV_ImpRK) ? 4 : maxord;
+
   /* copy input parameters into cv_mem */
   cv_mem->cv_lmm  = lmm;
   cv_mem->cv_iter = iter;
@@ -1542,7 +1546,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
       }
 
     } /* end of istop tests block */
-    if ((lmm==CV_DOPRI) || (lmm==CV_ExpRK)) mxstep = CVHinFixed(cv_mem, tout, tret);
+    if ((lmm==CV_DOPRI) || (lmm==CV_ExpRK) || (lmm==CV_ImpRK)) mxstep = CVHinFixed(cv_mem, tout, tret);
   } /* end stopping tests block */  
 
   /*
@@ -1627,6 +1631,9 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     case CV_ExpRK:
         kflag = CVStepExpRK(cv_mem);
         break;
+    case CV_ImpRK:
+        kflag = CVStepImpRK(cv_mem);
+        break;
     default:
        kflag = CVStep(cv_mem);
     }
@@ -2559,6 +2566,94 @@ static int CVStepExpRK(CVodeMem cv_mem)
   return(CV_SUCCESS);
 }
 
+/*
+ * CVStepImpRK
+ *
+ * This routine performs one internal cvode step using the implicit Runge-Kutta method, from tn to tn + h.
+ * In order to temporarily store the results, we use zn[2, 3, 4], tempv and ftemp, which will represent the Ki in turn.
+ */
+
+static int CVStepImpRK(CVodeMem cv_mem)
+{
+  int retval, nb_iter;
+  realtype difference;
+
+  /* Coefficients */
+  realtype a11, a21, a22, a31, a32, a33, b1, b2, b3, c1, c2, c3;
+  a11 =  0.377847764031163;
+  a21 =  0.385232756462588;
+  a22 =  0.461548399939329;
+  a31 =  0.675724855841358;
+  a32 = -0.061710969841169;
+  a33 =  0.241480233100410;
+  b1  =  0.750869573741408;
+  b2  = -0.362218781852651;
+  b3  =  0.611349208111243;
+  c1  =  0.257820901066211;
+  c2  =  0.434296446908075;
+  c3  =  0.758519768667167;
+
+  difference = 0;
+  nb_iter    = 1;
+  maxcor     = 30; /* Set maximum number of iterations */
+
+  /* Here, we use zn[2, 3, 4] to represent the Runge-Kutta coefficients K1, K2, K3.
+   * Set zn[1] = h*y'(tn) as the first guess for the K[i]. */
+  N_VScale (ONE, zn[1], zn[2]);
+  N_VScale (ONE, zn[1], zn[3]);
+  N_VScale (ONE, zn[1], zn[4]);
+
+  N_VLinearSum_Serial(ONE, zn[0], h*a11, zn[2], ftemp);   /* ftemp = a11K1 + Yn, */
+  retval = f(tn + c1*h, ftemp, zn[2], f_data);            /* K1 = f(tn+c1h, Yn + a11K1), */
+
+  N_VLinearSum_Serial(h*a21, zn[2], h*a22, zn[3], ftemp);  /* K2 = a21K1 + a22K2, */
+  N_VLinearSum_Serial(ONE, zn[0], ONE, ftemp, ftemp);      /* K2 = Yn + K2, */
+  retval = f(tn + c2*h, ftemp, zn[3], f_data);             /* K2 = f(tn+c2h, K2), */
+
+  N_VLinearSum_Serial(h*a32, zn[3], h*a33, zn[4], ftemp); /* K3 = a32K2 + a33K3, */
+  N_VLinearSum_Serial(h*a31, zn[2], ONE, ftemp, ftemp);    /* K3 = a31K1 + K3, */
+  N_VLinearSum_Serial(ONE, zn[0], ONE, ftemp, ftemp);      /* K3 = Yn + K3, */
+  retval = f(tn + c3*h, ftemp, zn[4], f_data);             /* K3 = f(tn+c3h, K3), */
+
+  N_VLinearSum_Serial(b2, zn[3], b3, zn[4], ftemp);       /* K3 = b2K2 + b3K3, */
+  N_VLinearSum_Serial(b1, zn[2], ONE, ftemp, ftemp);       /* K3 = b1K1 + K3, */
+  N_VLinearSum_Serial(ONE, zn[0], h, ftemp, tempv);        /* y = Yn+1 = Yn + K3 */
+
+  while (nb_iter <= maxcor) {  /* Same operations as above, but with K[i] updated and store result in y to compare with tempv */
+
+    N_VLinearSum_Serial(ONE, zn[0], h*a11, zn[2], ftemp);     /* ftemp = a11K1 + Yn, */
+    retval = f(tn + c1*h, ftemp, zn[2], f_data);              /* K1 = f(tn+c1h, Yn + a11K1), */
+    N_VLinearSum_Serial(h*a21, zn[2], h*a22, zn[3], ftemp);   /* K2 = a21K1 + a22K2, */
+    N_VLinearSum_Serial(ONE, zn[0], ONE, ftemp, ftemp);       /* K2 = Yn + K2, */
+    retval = f(tn + c2*h, ftemp, zn[3], f_data);              /* K2 = f(tn+c2h, K2), */
+
+    N_VLinearSum_Serial(h*a32, zn[3], h*a33, zn[4], ftemp);   /* K3 = a32K2 + a33K3, */
+    N_VLinearSum_Serial(h*a31, zn[2], ONE, ftemp, ftemp);      /* K3 = a31K1 + K3, */
+    N_VLinearSum_Serial(ONE, zn[0], ONE, ftemp, ftemp);        /* K3 = Yn + K3, */
+    retval = f(tn + c3*h, ftemp, zn[4], f_data);               /* K3 = f(tn+c3h, K3), */
+
+    N_VLinearSum_Serial(b2, zn[3], b3, zn[4], ftemp);     /* K3 = b2K2 + b3K3, */
+    N_VLinearSum_Serial(b1, zn[2], ONE, ftemp, ftemp);     /* K3 = b1K1 + K3, */
+    N_VLinearSum_Serial(ONE, zn[0], h, ftemp, y);          /* y = Yn+1 = Yn + K3 */
+
+    /* Convergence test */
+    N_VLinearSum_Serial(ONE, tempv, -ONE, y, ftemp);   /* ftemp = tempv-y, */
+    difference = N_VMaxNorm(ftemp);                    /* max = Max(ABS(ftemp)), */
+    if (difference < reltol) {  /* Converged */
+      tn += h;                               /* Increment tn, */
+      N_VScale (ONE, y, zn[0]);              /* Update Nordsziek array : - zn[0] = Yn+1, */
+      retval = f(tn, zn[0], zn[1], f_data);  /*                                                        - zn[1] = Y'(tn), */
+      N_VScale (h, zn[1], zn[1]);            /* Scale zn[1] by h */
+      return (CV_SUCCESS);
+    }
+    else {  /* Not converged yet, put y in tempv and reiterate */
+      N_VScale(ONE, y, tempv);
+      nb_iter++;
+    }
+  }
+  /* End of while : maxiter attained, we consider that the algorithm has diverged */
+  return (CONV_FAIL);
+}
 
 /*
  * CVAdjustParams
diff --git a/scilab/modules/scicos/tests/unit_tests/ImpRK.dia.ref b/scilab/modules/scicos/tests/unit_tests/ImpRK.dia.ref
new file mode 100644 (file)
index 0000000..92d0d68
--- /dev/null
@@ -0,0 +1,38 @@
+// =============================================================================
+// 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/ImpRK.tst");
+//  or test_run('scicos', 'ImpRK', ['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 ImpRK + save results
+ scs_m.props.tol(7) = 5*(10^-i); scs_m.props.tol(6) = 7;         // 'max step size' + solver
+ scs_m.props.tol(2) = 1.0e-12;                                   // reltol
+ try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end  // ImpRK
+ rkval = res.values;   // Results
+ time = res.time;      // Time
+ // Modify solver and reltol + run CVode + save results
+ scs_m.props.tol(6) = 1; scs_m.props.tol(2) = 1.0e-15;
+ 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);
+ assert_checktrue(mea <= 10^-i);
+ assert_checktrue(stdeviation <= 10^-i);
+end
diff --git a/scilab/modules/scicos/tests/unit_tests/ImpRK.tst b/scilab/modules/scicos/tests/unit_tests/ImpRK.tst
new file mode 100644 (file)
index 0000000..3232c74
--- /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/ImpRK.tst");
+//  or test_run('scicos', 'ImpRK', ['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 ImpRK + save results
+ scs_m.props.tol(7) = 5*(10^-i); scs_m.props.tol(6) = 7;         // 'max step size' + solver
+ scs_m.props.tol(2) = 1.0e-12;                                   // reltol
+ try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end  // ImpRK
+ rkval = res.values;   // Results
+ time = res.time;      // Time
+
+ // Modify solver and reltol + run CVode + save results
+ scs_m.props.tol(6) = 1; scs_m.props.tol(2) = 1.0e-15;
+ 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);
+ assert_checktrue(mea <= 10^-i);
+ assert_checktrue(stdeviation <= 10^-i);
+
+end
index 1d6d58b..10a42bb 100644 (file)
@@ -26,6 +26,7 @@
                 <solver code="4" description="Sundials/CVODE - ADAMS - FUNCTIONAL"/>
                 <solver code="5" description="DOPRI5 - Dormand-Prince 4(5)"/>
                 <solver code="6" description="RK45 - Runge-Kutta 4(5)"/>
+                <solver code="7" description="Implicit RK45 - Runge-Kutta 4(5) - FIXED-POINT"/>
                 <solver code="100" description="Sundials/IDA"/>
                 <!-- Available trace values lists from scicos -->        
                 <trace code="0" description="_(No trace nor debug printing)"/>
index 8069e08..007e27d 100644 (file)
Binary files a/scilab/modules/xcos/examples/solvers/ODE_Example.xcos and b/scilab/modules/xcos/examples/solvers/ODE_Example.xcos differ
diff --git a/scilab/modules/xcos/examples/solvers/integImpRK.sce b/scilab/modules/xcos/examples/solvers/integImpRK.sce
new file mode 100644 (file)
index 0000000..ac051f4
--- /dev/null
@@ -0,0 +1,35 @@
+// 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/integRK.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", "implicit Runge-Kutta"];
+
+for solver=1:5
+
+ // Select the solver (implicit Runge-Kutta is solver number 7)
+ scs_m.props.tol(6) = solver;
+ if (solver == 5) then scs_m.props.tol(6) = 7; end
+
+ // Set max step size and reltol if implicit Runge-Kutta
+ if (solver == 5) then scs_m.props.tol(7) = 0.01; scs_m.props.tol(2) = 1.0e-10; end
+
+ // 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 1167a04..74cf7bf 100644 (file)
@@ -28,7 +28,7 @@ for solver=1:5
 
  // Start the timer, launch the simulation and display time
  tic();
- try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
+ try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end
  t = toc();
  disp(t, "Time for " + solverName(solver) + " :");
 
index 81f3c4c..8601f06 100644 (file)
@@ -89,7 +89,7 @@ loadScicos();
 loadXcosLibs();
 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
 scs_m.props.tol(6) = 0;
-try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
@@ -101,14 +101,14 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
       loadScicos();
       loadXcosLibs();
       importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
-      scs_m.props.tf = 10000;
+      scs_m.props.tf = 5000;
 
       // Select the LSodar solver
       scs_m.props.tol(6) = 0;
 
       // Start the timer, launch the simulation and display time
       tic();
-      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
       t = toc();
       disp(t, "Time for LSodar :");
       ]]></programlisting>
@@ -117,7 +117,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
             The Scilab console displays :
             <screen><![CDATA[
 Time for LSodar :
- 30.584
+ 16.836
             ]]></screen>
         </para>
         <para>
@@ -129,19 +129,19 @@ Time for LSodar :
         <para>
             <screen><![CDATA[
 Time for LSodar :
- 5.001
+ 6.115
 
 Time for BDF / Newton :
- 12.943
+ 19.894
 
 Time for BDF / Functional :
- 12.415
+ 19.382
 
 Time for Adams / Newton :
- 8.651
+ 11.255
 
 Time for Adams / Functional :
- 7.962
+ 10.233
             ]]></screen>
         </para>
         <para>
@@ -171,6 +171,9 @@ Time for Adams / Functional :
                 <link linkend="DoPri">Dormand-Price 4(5)</link>
             </member>
             <member>
+                <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
index 7be923f..f3d94ff 100644 (file)
@@ -192,7 +192,7 @@ loadScicos();
 loadXcosLibs();
 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
 scs_m.props.tol(6) = 1;
-try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
@@ -204,14 +204,14 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
       loadScicos();
       loadXcosLibs();
       importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
-      scs_m.props.tf = 10000;
+      scs_m.props.tf = 5000;
 
       // Select the solver BDF / Newton
       scs_m.props.tol(6) = 1;
 
       // Start the timer, launch the simulation and display time
       tic();
-      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
       t = toc();
       disp(t, "Time for BDF / Newton :");
       ]]></programlisting>
@@ -220,7 +220,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
             The Scilab console displays :
             <screen><![CDATA[
 Time for BDF / Newton :
- 33.997
+ 19.454
             ]]></screen>
         </para>
         <para>
@@ -233,16 +233,16 @@ Time for BDF / Newton :
             Results :
             <screen language="example"><![CDATA[
 Time for BDF / Newton :
- 12.998
+ 19.894
 
 Time for BDF / Functional :
- 12.256
+ 19.382
 
 Time for Adams / Newton :
- 8.577
+ 11.255
 
 Time for Adams / Functional :
- 8.002
+ 10.233
             ]]></screen>
         </para>
         <para>
@@ -265,6 +265,9 @@ Time for Adams / Functional :
                 <link linkend="DoPri">Dormand-Price 4(5)</link>
             </member>
             <member>
+                <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
index 28891ff..034fd7d 100644 (file)
@@ -149,7 +149,7 @@ loadXcosLibs();
 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
 scs_m.props.tol(6) = 6;
 scs_m.props.tol(7) = 10^-2;
-try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
@@ -161,7 +161,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
       loadScicos();
       loadXcosLibs();
       importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
-      scs_m.props.tf = 10000;
+      scs_m.props.tf = 5000;
 
       // Select the solver Runge-Kutta and set the precision
       scs_m.props.tol(6) = 6;
@@ -169,7 +169,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
 
       // Start the timer, launch the simulation and display time
       tic();
-      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
       t = toc();
       disp(t, "Time for Runge-Kutta :");
       ]]></programlisting>
@@ -178,7 +178,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
             The Scilab console displays :
             <screen><![CDATA[
 Time for Runge-Kutta :
- 18.254
+ 15.624
             ]]></screen>
         </para>
         <para>
@@ -190,26 +190,26 @@ Time for Runge-Kutta :
         <para>
             <screen><![CDATA[
 Time for BDF / Newton :
- 25.457
+ 19.894
 
 Time for BDF / Functional :
- 24.893
+ 19.382
 
 Time for Adams / Newton :
- 20.469
+ 11.255
 
 Time for Adams / Functional :
- 20.049
+ 10.233
 
 Time for Runge-Kutta :
- 18.254
+ 5.502
             ]]></screen>
         </para>
         <para>
             These results show that on a nonstiff problem, for relatively same precision required and forcing the same step size, Runge-Kutta is faster.
         </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.
+            Variable-size step 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>
@@ -228,6 +228,9 @@ Time for Runge-Kutta :
                 <link linkend="DoPri">Dormand-Price 4(5)</link>
             </member>
             <member>
+                <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
index a23bd06..51ed626 100644 (file)
@@ -133,7 +133,7 @@ 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;
+try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
@@ -153,7 +153,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
 
       // Start the timer, launch the simulation and display time
       tic();
-      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
       t = toc();
       disp(t, "Time for Dormand-Price :");
       ]]></programlisting>
@@ -162,7 +162,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
             The Scilab console displays :
             <screen><![CDATA[
 Time for Dormand-Price :
- 17.126
+ 17.226
             ]]></screen>
         </para>
         <para>
@@ -174,26 +174,26 @@ Time for Dormand-Price :
         <para>
             <screen><![CDATA[
 Time for BDF / Newton :
- 12.852
+ 19.894
 
 Time for BDF / Functional :
- 12.343
+ 19.382
 
 Time for Adams / Newton :
- 8.75
+ 11.255
 
 Time for Adams / Functional :
- 7.789
+ 10.233
 
 Time for DoPri :
- 7.606
+ 8.654
             ]]></screen>
         </para>
         <para>
             These results show that on a nonstiff problem, for relatively same precision required and forcing the same step size, Dormand-Price's computational overhead (compared to Runge-Kutta) is significant and is close to <emphasis>Adams/Functional</emphasis>. 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.
+            Variable-size step 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>
@@ -212,6 +212,9 @@ Time for DoPri :
                 <link linkend="RK">Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
diff --git a/scilab/modules/xcos/help/en_US/solvers/4-ImplicitRK.xml b/scilab/modules/xcos/help/en_US/solvers/4-ImplicitRK.xml
new file mode 100644 (file)
index 0000000..8523aca
--- /dev/null
@@ -0,0 +1,305 @@
+<?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="ImpRK">
+    <refnamediv>
+        <refname>Implicit Runge-Kutta 4(5)</refname>
+        <refpurpose>
+            Runge-Kutta is a numerical solver providing an efficient and stable implicit method to solve Ordinary Differential Equations (ODEs) Initial Value Problems. Called by <link linkend="xcos">xcos</link>.
+        </refpurpose>
+    </refnamediv>
+    <refsection>
+        <title>Description</title>
+        <para>
+            Runge-Kutta is a numerical solver providing an efficient and stable 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 Runge-Kutta, we do not adapt to the complexity of the problem, but we guarantee a stable computation time.
+        </para>
+        <para>
+            This method being explicit, it can be used on stiff problems.
+        </para>
+        <para>
+            It is an enhancement of the backward Euler method, which approximates 
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+            by computing 
+            <emphasis>
+                f(t<subscript>n</subscript>+h, y<subscript>n+1</subscript>)
+            </emphasis>
+            and truncating the Taylor expansion.
+        </para>
+        <para>
+            The implemented scheme is inspired from the "Low-Dispersion Low-Dissipation Implicit Runge-Kutta Scheme" (see bottom for link).
+        </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 implicit Runge-Kutta with the previous methods is that it computes up to the fourth derivative of <emphasis>y</emphasis>, while the others mainly 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 three 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>. They are distributed approximately equally on the interval.
+            <itemizedlist>
+                <listitem>
+                    <emphasis>k1</emphasis> is the increment based on the slope near the quarter of the interval, using 
+                    <emphasis>
+                        y<subscript>n</subscript>+ a11*h*k1,
+                    </emphasis>
+                    ,
+                </listitem>
+                <listitem>
+                    <emphasis>k2</emphasis> is the increment based on the slope near the midpoint of the interval, using 
+                    <emphasis>
+                        y<subscript>n</subscript> + a21*h*k1 + a22*h*k2,
+                    </emphasis>
+                    ,
+                </listitem>
+                <listitem>
+                    <emphasis>k3</emphasis> is the increment based on the slope near the third quarter of the interval, using 
+                    <emphasis>
+                        y<subscript>n</subscript> + a31*h*k1 + a32*h*k2 + a33*h*k3.
+                    </emphasis>
+                </listitem>
+            </itemizedlist>
+        </para>
+        <para>
+            We see that the computation of <emphasis>ki</emphasis> requires <emphasis>ki</emphasis>, thus necessitating the use of a nonlinear solver (here, fixed-point iterations).
+        </para>
+        <para>
+            First, we set 
+            <emphasis>
+                k0 = h * f(t<subscript>n</subscript>, y<subscript>n</subscript>)
+            </emphasis>
+            as first guess for all the <emphasis>ki</emphasis>, to get updated <emphasis>ki</emphasis> and a first value for 
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+            .
+        </para>
+        <para>
+            Next, we save and recompute 
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+            with those new <emphasis>ki</emphasis>.
+        </para>
+        <para>
+            Then, we compare the two 
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+            and recompute it until its difference with the last computed one is inferior to the simulation parameter <emphasis>reltol</emphasis>.
+        </para>
+        <para>
+            This process adds a significant computation time to the method, but greatly improves stability.
+        </para>
+        <para>
+            We can see that with the <emphasis>ki</emphasis>, we progress in the derivatives of 
+            <emphasis>
+                y<subscript>n</subscript>
+            </emphasis>
+            . So in <emphasis>k3</emphasis>, we are approximating 
+            <emphasis>
+                y<superscript>(3)</superscript><subscript>n</subscript>
+            </emphasis>
+            , thus making an error in 
+            <emphasis>
+                O(h<superscript>4</superscript>)
+            </emphasis>
+            . But choosing the right coefficients in the computation of the <emphasis>ki</emphasis> (notably the 
+            <emphasis>
+                a<subscript>ij</subscript>
+            </emphasis>
+            ) makes us gain an order, thus making a per-step total error in 
+            <emphasis>
+                O(h<superscript>5</superscript>)
+            </emphasis>
+            .
+        </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>implicit Runge-Kutta 4(5)</emphasis> : 
+            <emphasis>
+                O(h<superscript>5</superscript>)
+            </emphasis>
+            per step, 
+            <emphasis>
+                O(h<superscript>4</superscript>)
+            </emphasis>
+            in total.
+        </para>
+        <para>
+            Although 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 we approach 
+            <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(2) = 10^-5;
+scs_m.props.tol(6) = 7;
+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 Runge-Kutta 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 = 5000;
+
+      // Select the solver Runge-Kutta and set the precision
+      scs_m.props.tol(2) = 10^-10;
+      scs_m.props.tol(6) = 7;
+      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 implicit Runge-Kutta :");
+      ]]></programlisting>
+        </para>
+        <para>
+            The Scilab console displays :
+            <screen><![CDATA[
+Time for implicit Runge-Kutta :
+ 17.003
+            ]]></screen>
+        </para>
+        <para>
+            Now, in the following script, we compare the time difference between implicit RK and CVode by running the example with the five solvers in turn :
+            <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/integImpRK.sce">
+                Open the script
+            </link>
+        </para>
+        <para>
+            <screen><![CDATA[
+Time for BDF / Newton :
+ 19.894
+
+Time for BDF / Functional :
+ 19.382
+
+Time for Adams / Newton :
+ 11.255
+
+Time for Adams / Functional :
+ 10.233
+
+Time for implicit Runge-Kutta :
+ 8.112
+            ]]></screen>
+        </para>
+        <para>
+            These results show that on a nonstiff problem, for relatively same precision required and forcing the same step size, implicit Runge-Kutta competes with Adams / Functional.
+        </para>
+        <para>
+            Variable-size step 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="RK">Runge-Kutta 4(5)</link>
+            </member>
+            <member>
+                <link linkend="DoPri">Dormand-Prince 4(5)</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>
+            ACM SIGNUM Newsletter, Volume 15, Issue 4, December 1980, Pages 10-11 <ulink url="http://dl.acm.org/citation.cfm?id=1218052.1218054&amp;coll=DL&amp;dl=GUIDE&amp;CFID=213153447&amp;CFTOKEN=26634355">LSode - LSodi</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>Implicit Runge-Kutta 4(5) solver added</revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+</refentry>
index 1b129b0..6a1da87 100644 (file)
@@ -153,7 +153,7 @@ function messagebox(msg, title)
  disp(title);
  disp(msg);
 endfunction
-try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
     </refsection>
@@ -173,6 +173,9 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
                 <link linkend="DoPri">Dormand-Price 4(5)</link>
             </member>
             <member>
+                <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
index a72ee25..fa12a24 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-Prince 4(5)", "RK45 - Runge-Kutta 4(5)", "Sundials/IDA"
+                                                "Sundials/CVODE - ADAMS - NEWTON", "Sundials/CVODE - ADAMS - FUNCTIONAL", "DOPRI5 - Dormand-Prince 4(5)", "RK45 - Runge-Kutta 4(5)", "Implicit 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", "Method: Fixed step", "Sundials/IDA"
+                "Method: Fixed step", "Method: Fixed step", "Method: Fixed step, Nonlinear solver= FIXED-POINT", "Sundials/IDA"
                                                       };
 
         solver = new JComboBox(solvers);
index df65868..29c342d 100644 (file)
@@ -304,6 +304,7 @@ public class ScicosParameters implements Serializable, Cloneable {
      * <li>4 : Sundials/CVODE : Method: ADAMS, Nonlinear solver= FUNCTIONAL
      * <li>5 : DOPRI5 : Method: Dormand-Prince 4(5)
      * <li>6 : RK45 : Method: Runge-Kutta 4(5)
+     * <li>7 : Implicit RK45 : Method: Runge-Kutta 4(5), Nonlinear solver= Fixed-point
      * <li>100 : Sundials/IDA
      *
      *
@@ -322,6 +323,7 @@ public class ScicosParameters implements Serializable, Cloneable {
      * <li>4 : Sundials/CVODE : Method: ADAMS, Nonlinear solver= FUNCTIONAL
      * <li>5 : DOPRI5 : Method: Dormand-Prince 4(5)
      * <li>6 : RK45 : Method: Runge-Kutta 4(5)
+     * <li>7 : Implicit RK45 : Method: Runge-Kutta 4(5), Nonlinear solver= FIXED-POINT
      * <li>100 : Sundials/IDA
      *
      *