* Bug #13308 fixed - Xcos solvers: Crank-Nicolson implem and interface 77/14377/14
Paul Bignier [Tue, 8 Apr 2014 10:36:08 +0000 (12:36 +0200)]
Change-Id: I80a09f831a07d57618abfa99efa9fea7e57023b3

35 files changed:
scilab/CHANGES
scilab/modules/helptools/etc/images_md5.txt
scilab/modules/helptools/images/5-CrankNicolson_1.png [new file with mode: 0644]
scilab/modules/helptools/images/_LaTeX_5-CrankNicolson.xml_1.png [new file with mode: 0644]
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/Solvers/ODE/CRANI.dia.ref [new file with mode: 0644]
scilab/modules/scicos/tests/unit_tests/Solvers/ODE/CRANI.tst [new file with mode: 0644]
scilab/modules/xcos/etc/XConfiguration-xcos.xml
scilab/modules/xcos/examples/solvers/benchBasic.sce
scilab/modules/xcos/examples/solvers/benchKalman.sce
scilab/modules/xcos/examples/solvers/benchSine.sce
scilab/modules/xcos/examples/solvers/integCRANI.sce [new file with mode: 0644]
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-Prince.xml
scilab/modules/xcos/help/en_US/solvers/4-ImplicitRK.xml
scilab/modules/xcos/help/en_US/solvers/5-CrankNicolson.xml [new file with mode: 0644]
scilab/modules/xcos/help/en_US/solvers/6-IDA.xml
scilab/modules/xcos/help/en_US/solvers/7-DDaskr.xml
scilab/modules/xcos/help/en_US/solvers/8-Rootfinding.xml
scilab/modules/xcos/help/en_US/solvers/9-Comparisons.xml
scilab/modules/xcos/help/fr_FR/solvers/0-LSodar.xml
scilab/modules/xcos/help/fr_FR/solvers/1-CVode.xml
scilab/modules/xcos/help/fr_FR/solvers/2-Runge-Kutta.xml
scilab/modules/xcos/help/fr_FR/solvers/3-Dormand-Prince.xml
scilab/modules/xcos/help/fr_FR/solvers/4-ImplicitRK.xml
scilab/modules/xcos/help/fr_FR/solvers/5-CrankNicolson.xml [new file with mode: 0644]
scilab/modules/xcos/help/fr_FR/solvers/6-IDA.xml
scilab/modules/xcos/help/fr_FR/solvers/7-DDaskr.xml
scilab/modules/xcos/help/fr_FR/solvers/8-Racines.xml
scilab/modules/xcos/help/fr_FR/solvers/9-Comparaisons.xml
scilab/modules/xcos/src/java/org/scilab/modules/xcos/actions/dialog/SetupDialog.java

index 080d742..173c6f2 100644 (file)
@@ -266,6 +266,13 @@ Modified Functions
     "-" + mgetl("TMPDIR/test") + "-"
 
 
+Xcos
+====
+
+* Implicit fixed-size step ODE solver added: Crank-Nicolson 2(3).
+  Added to the CVode package, it also benefits from the CVode rootfinding feature.
+
+
 Obsolete functions or features
 ==============================
 * maxfiles function is now obsolete.
@@ -357,6 +364,8 @@ In 6.0.0:
 
 * Bug #13154       - In shellmode, completion now separates Files from Directories.
 
+* Bug #13308 fixed - Added the Crank-Nicolson solver to Xcos.
+
 * Bug #13465 fixed - The display of polyline .display_function and .display_function properties was not conventional
 
 * Bug #13468 fixed - Scilab hanged when incorrect format was used for file reading using mfscanf.
index 197fd5c..bcb0eee 100644 (file)
@@ -4,6 +4,7 @@
 3-Dormand-Prince_1.png=f03f6156c4442805827ba389af1f314e
 4-ImplicitRK_1.png=3baa7637cb0be94076c6e06b7ae34991
 4-RKImplicite_1.png=3baa7637cb0be94076c6e06b7ae34991
+5-CrankNicolson_1.png=e460c3c5dcdf57152d190386d91eac1d
 6-IDA_en_US_1.png=56248160f1b98de82bccc2bbd1850f5
 6-IDA_fr_FR_1.png=c0c4c44ea933b98ca3a7111ff32c5066
 6-IDA_ja_JP_1.png=56248160f1b98de82bccc2bbd1850f5
@@ -187,6 +188,7 @@ _LaTeX_2-Runge-Kutta.xml_1.png=e130ce3936e3d90b95bf43660f09cb4e
 _LaTeX_3-Dormand-Prince.xml_1.png=e130ce3936e3d90b95bf43660f09cb4e
 _LaTeX_4-ImplicitRK.xml_1.png=e130ce3936e3d90b95bf43660f09cb4e
 _LaTeX_4-RKImplicite.xml_1.png=e130ce3936e3d90b95bf43660f09cb4e
+_LaTeX_5-CrankNicolson.xml_1.png=e130ce3936e3d90b95bf43660f09cb4e
 _LaTeX_6-IDA.xml_1.png=dff7132b713410ac0eb11a45b1d627dc
 _LaTeX_6-IDA.xml_2.png=1aa0a943d7130da822dd44a06bd29841
 _LaTeX_6-IDA.xml_3.png=dc4144bc7f206b73cb21d076e514318
diff --git a/scilab/modules/helptools/images/5-CrankNicolson_1.png b/scilab/modules/helptools/images/5-CrankNicolson_1.png
new file mode 100644 (file)
index 0000000..de86bf2
Binary files /dev/null and b/scilab/modules/helptools/images/5-CrankNicolson_1.png differ
diff --git a/scilab/modules/helptools/images/_LaTeX_5-CrankNicolson.xml_1.png b/scilab/modules/helptools/images/_LaTeX_5-CrankNicolson.xml_1.png
new file mode 100644 (file)
index 0000000..e660e5b
Binary files /dev/null and b/scilab/modules/helptools/images/_LaTeX_5-CrankNicolson.xml_1.png differ
index 01cdcc7..27b679d 100644 (file)
@@ -113,6 +113,7 @@ enum Solver
     Dormand_Prince,
     Runge_Kutta,
     Implicit_Runge_Kutta,
+    Crank_Nicolson,
     IDA_BDF_Newton = 100,
     DDaskr_BDF_Newton = 101,
     DDaskr_BDF_GMRes = 102
@@ -880,6 +881,7 @@ int C2F(scicos)(double *x_in, int *xptr_in, double *z__,
             case Dormand_Prince:
             case Runge_Kutta:
             case Implicit_Runge_Kutta:
+            case Crank_Nicolson:
                 cossim(t0);
                 break;
             case IDA_BDF_Newton:
@@ -1472,6 +1474,7 @@ static void cossim(double *told)
         case Dormand_Prince:
         case Runge_Kutta:
         case Implicit_Runge_Kutta:
+        case Crank_Nicolson:
             ODEFree = &CVodeFree;
             ODE = &CVode;
             ODEReInit = &CVodeReInit;
@@ -1567,6 +1570,9 @@ static void cossim(double *told)
             case Implicit_Runge_Kutta:
                 ode_mem = CVodeCreate(CV_ImpRK, CV_FUNCTIONAL);
                 break;
+            case Crank_Nicolson:
+                ode_mem = CVodeCreate(CV_CRANI, CV_FUNCTIONAL);
+                break;
         }
 
         /*    ode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);*/
@@ -3643,9 +3649,9 @@ void callf(double *t, scicos_block *block, scicos_flag *flag)
     //sciprint("callf type=%d flag=%d\n",block->type,flagi);
     switch (block->type)
     {
-        /*******************/
-        /* function type 0 */
-        /*******************/
+            /*******************/
+            /* function type 0 */
+            /*******************/
         case 0 :
         {
             /* This is for compatibility */
index 3de4805..fe30759 100644 (file)
@@ -80,6 +80,7 @@ extern "C" {
 #define CV_DOPRI 3
 #define CV_ExpRK 4
 #define CV_ImpRK 5
+#define CV_CRANI 6
 
     /* iter */
 #define CV_FUNCTIONAL 1
index f33e0f5..366a13b 100644 (file)
@@ -265,6 +265,7 @@ 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 CVStepCRANI(CVodeMem cv_mem);
 
 static int CVsldet(CVodeMem cv_mem);
 
@@ -354,7 +355,7 @@ void *CVodeCreate(int lmm, int iter)
 
     /* Test inputs */
 
-    if ((lmm != CV_ADAMS) && (lmm != CV_BDF) && (lmm != CV_DOPRI) && (lmm != CV_ExpRK) && (lmm != CV_ImpRK))   /* Integration mode: ADAMS, BDF or RK-based */
+    if ((lmm != CV_ADAMS) && (lmm != CV_BDF) && (lmm != CV_DOPRI) && (lmm != CV_ExpRK) && (lmm != CV_ImpRK) && (lmm != CV_CRANI))   /* Integration mode: ADAMS, BDF or RK-based */
     {
         CVProcessError(NULL, 0, "CVODE", "CVodeCreate", MSGCV_BAD_LMM);
         return(NULL);
@@ -388,6 +389,9 @@ void *CVodeCreate(int lmm, int iter)
     /* 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;
 
+    /* If Crank-Nicolson is selected, then maxord = 3 to use the 2 extra vectors allocated (zn[2, 3]) */
+    maxord = (lmm == CV_CRANI) ? 3 : maxord;
+
     /* copy input parameters into cv_mem */
     cv_mem->cv_lmm  = lmm;
     cv_mem->cv_iter = iter;
@@ -1493,7 +1497,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
             }
 
         }
-        if ((lmm == CV_DOPRI) || (lmm == CV_ExpRK) || (lmm == CV_ImpRK))
+        if ((lmm == CV_DOPRI) || (lmm == CV_ExpRK) || (lmm == CV_ImpRK) || (lmm == CV_CRANI))
         {
             mxstep = CVHinFixed(cv_mem, tout, tret);
         }
@@ -1600,6 +1604,9 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
             case CV_ImpRK:
                 kflag = CVStepImpRK(cv_mem);
                 break;
+            case CV_CRANI:
+                kflag = CVStepCRANI(cv_mem);
+                break;
             default:
                 kflag = CVStep(cv_mem);
         }
@@ -2753,6 +2760,78 @@ static int CVStepImpRK(CVodeMem cv_mem)
 }
 
 /*
+ * CVStepCRANI
+ *
+ * This routine performs one internal cvode step using the Crank-Nicolson method, from tn to tn + h.
+ * In order to temporarily store the results, we use zn[2, 3], tempv and ftemp, which will represent the Ki in turn.
+ */
+
+static int CVStepCRANI(CVodeMem cv_mem)
+{
+    int retval, nb_iter;
+    realtype difference;
+
+    /* Coefficients */
+    realtype a11, a12, b1, b2, c1;
+    a11 = 0.5;
+    a12 = 0.5;
+    b1  = 0.5;
+    b2  = 0.5;
+    c1  = 1;
+
+    difference = 0;  /* Difference between two computed solutions */
+    nb_iter    = 1;  /* Iterations counter                        */
+    maxcor     = 30; /* Set maximum number of iterations          */
+
+    /* Here, we use zn[2, 3] to represent the Runge-Kutta coefficients K1, K2.
+     * 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_VLinearSum_Serial(ONE, zn[0], h * a11, zn[2], ftemp); /* ftemp = a11hK1 + Yn          */
+    N_VLinearSum_Serial(ONE, ftemp, h * a12, zn[3], ftemp); /* ftemp = a11hK1 + a12hK2 + Yn */
+
+    retval = f(tn + c1 * h, ftemp, zn[2], user_data); /* K1 = f(tn+c1h, Yn + a11hK1 + a12hK2) */
+    retval = f(tn, zn[0], zn[3], user_data);          /* K2 = f(tn, Yn)                       */
+
+    N_VLinearSum_Serial(b1, zn[2], b2, zn[3], ftemp); /* ftemp = b1K1 + b2K2            */
+    N_VLinearSum_Serial(ONE, zn[0], h, ftemp, tempv); /* y = Yn+1 = Yn + h(b1K1 + b2K2) */
+
+    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 = a11hK1 + Yn          */
+        N_VLinearSum_Serial(ONE, ftemp, h * a12, zn[3], ftemp); /* ftemp = a11hK1 + a12hK2 + Yn */
+
+        retval = f(tn + c1 * h, ftemp, zn[2], user_data); /* K1 = f(tn+c1h, Yn + a11hK1 + a12hK2) */
+        retval = f(tn, zn[0], zn[3], user_data);          /* K2 = f(tn, Yn)                       */
+
+        N_VLinearSum_Serial(b1, zn[2], b2, zn[3], ftemp); /* ftemp = b1K1 + b2K2            */
+        N_VLinearSum_Serial(ONE, zn[0], h, ftemp, y);     /* y = Yn+1 = Yn + h(b1K1 + b2K2) */
+
+        /* 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                             */
+            nst++;                                   /* Increment solver calls (complete step)   */
+            N_VScale (ONE, y, zn[0]);                /* Update Nordsziek array: - zn[0] = Yn+1   */
+            retval = f(tn, zn[0], zn[1], user_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
  *
  * This routine is called when a change in step size was decided upon,
diff --git a/scilab/modules/scicos/tests/unit_tests/Solvers/ODE/CRANI.dia.ref b/scilab/modules/scicos/tests/unit_tests/Solvers/ODE/CRANI.dia.ref
new file mode 100644 (file)
index 0000000..5fb0718
--- /dev/null
@@ -0,0 +1,38 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2016 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- XCOS TEST -->
+//
+// Testing the Crank-Nicolson method.
+// Import diagram
+assert_checktrue(importXcosDiagram("SCI/modules/xcos/tests/unit_tests/Solvers/ODE/Kalman.zcos"));
+Info = scicos_simulate(scs_m, list());
+for i=2:4  // 'max step size' = 10^-i, precision
+    // Start by updating the clock block period (sampling)
+    scs_m.objs(9).model.rpar(1) = 5*10^-i;
+    scs_m.objs(10).model.rpar(1) = 5*10^-i;
+    // Modify solver and 'max step size' + run CRANI + save results
+    scs_m.props.tol(7) = 5*10^-i; // 'max step size'
+    scs_m.props.tol(6) = 8;       // Solver
+    scs_m.props.tol(2) = 1d-12;   // reltol
+    try scicos_simulate(scs_m, Info); catch disp(lasterror()); end  // CRANI
+    rkval = res.values;   // Results
+    // Modify solver and reltol + run CVode + save results
+    scs_m.props.tol(6) = 4; scs_m.props.tol(2) = 1d-15;
+    try scicos_simulate(scs_m, Info); 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 = stdev(compa);
+    // Verifying closeness of the results
+    assert_checktrue(maxi <= 5*10^-i);
+    assert_checktrue(mea <= 10^-i);
+    assert_checktrue(stdeviation <= 10^-i);
+end
diff --git a/scilab/modules/scicos/tests/unit_tests/Solvers/ODE/CRANI.tst b/scilab/modules/scicos/tests/unit_tests/Solvers/ODE/CRANI.tst
new file mode 100644 (file)
index 0000000..08bac9e
--- /dev/null
@@ -0,0 +1,47 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2016 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- XCOS TEST -->
+//
+// Testing the Crank-Nicolson method.
+
+// Import diagram
+assert_checktrue(importXcosDiagram("SCI/modules/xcos/tests/unit_tests/Solvers/ODE/Kalman.zcos"));
+Info = scicos_simulate(scs_m, list());
+
+for i=2:4  // 'max step size' = 10^-i, precision
+
+    // Start by updating the clock block period (sampling)
+    scs_m.objs(9).model.rpar(1) = 5*10^-i;
+    scs_m.objs(10).model.rpar(1) = 5*10^-i;
+
+    // Modify solver and 'max step size' + run CRANI + save results
+    scs_m.props.tol(7) = 5*10^-i; // 'max step size'
+    scs_m.props.tol(6) = 8;       // Solver
+    scs_m.props.tol(2) = 1d-12;   // reltol
+    try scicos_simulate(scs_m, Info); catch disp(lasterror()); end  // CRANI
+    rkval = res.values;   // Results
+
+    // Modify solver and reltol + run CVode + save results
+    scs_m.props.tol(6) = 4; scs_m.props.tol(2) = 1d-15;
+    try scicos_simulate(scs_m, Info); 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 = stdev(compa);
+
+    // Verifying closeness of the results
+    assert_checktrue(maxi <= 5*10^-i);
+    assert_checktrue(mea <= 10^-i);
+    assert_checktrue(stdeviation <= 10^-i);
+
+end
index e9949d0..52205f7 100644 (file)
                 <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="7" description="Implicit RK45 - Implicit Runge-Kutta 4(5) - FIXED-POINT"/>
+                <solver code="8" description="CRANI - Crank-Nicolson 2(3) - FIXED-POINT"/>
                 <solver code="100" description="Sundials/IDA"/>
                 <solver code="101" description="DDaskr - BDF"/>
-                <!-- Available trace values lists from scicos -->        
+                <!-- Available trace values lists from scicos -->
                 <trace code="0" description="_(No trace nor debug printing)"/>
                 <trace code="1" description="_(Light Simulation trace (Discrete and Continuous part switches))"/>
                 <trace code="2" description="_(Per block execution trace and Debug block calls)"/>
index f26cbae..bdf7ccf 100644 (file)
@@ -14,10 +14,10 @@ tolerances(1) = 1d-13;
 tolerances(2) = 1d-13;
 [%tcur, %cpr, alreadyran, needstart, needcompile, %state0] = Info(:);
 
-solverName = ["LSodar", "CVode BDF/Newton", "CVode BDF/Functional", "CVode Adams/Newton", "CVode Adams/Functional", "Dormand-Prince", "Runge-Kutta", "implicit Runge-Kutta"];
+solverName = ["LSodar", "CVode BDF/Newton", "CVode BDF/Functional", "CVode Adams/Newton", "CVode Adams/Functional", "Dormand-Prince", "Runge-Kutta", "implicit Runge-Kutta", "Crank-Nicolson"];
 
 disp("--------------------------------");
-for solver = 0:7
+for solver = 0:8
 
     disp("Time for " + solverName(solver + 1) + ":");
     tolerances(6) = solver;
index 043c832..1a2c88f 100644 (file)
@@ -14,10 +14,10 @@ tolerances(1) = 1d-13;
 tolerances(2) = 1d-13;
 [%tcur, %cpr, alreadyran, needstart, needcompile, %state0] = Info(:);
 
-solverName = ["LSodar", "CVode BDF/Newton", "CVode BDF/Functional", "CVode Adams/Newton", "CVode Adams/Functional", "Dormand-Prince", "Runge-Kutta", "implicit Runge-Kutta"];
+solverName = ["LSodar", "CVode BDF/Newton", "CVode BDF/Functional", "CVode Adams/Newton", "CVode Adams/Functional", "Dormand-Prince", "Runge-Kutta", "implicit Runge-Kutta", "Crank-Nicolson"];
 
 disp("--------------------------------");
-for solver = 0:7
+for solver = 0:8
 
     disp("Time for " + solverName(solver + 1) + ":");
     tolerances(6) = solver;
index 4c16415..e9603fb 100644 (file)
@@ -12,10 +12,10 @@ tf = 50000;
 tolerances = scs_m.props.tol;
 [%tcur, %cpr, alreadyran, needstart, needcompile, %state0] = Info(:);
 
-solverName = ["LSodar", "CVode BDF/Newton", "CVode BDF/Functional", "CVode Adams/Newton", "CVode Adams/Functional", "Dormand-Prince", "Runge-Kutta", "implicit Runge-Kutta"];
+solverName = ["LSodar", "CVode BDF/Newton", "CVode BDF/Functional", "CVode Adams/Newton", "CVode Adams/Functional", "Dormand-Prince", "Runge-Kutta", "implicit Runge-Kutta", "Crank-Nicolson"];
 
 disp("--------------------------------");
-for solver = 0:7
+for solver = 0:8
 
     disp("Time for " + solverName(solver + 1) + ":");
     tolerances(6) = solver;
diff --git a/scilab/modules/xcos/examples/solvers/integCRANI.sce b/scilab/modules/xcos/examples/solvers/integCRANI.sce
new file mode 100644 (file)
index 0000000..87b3ae7
--- /dev/null
@@ -0,0 +1,29 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2016 - Scilab Enterprises - Paul Bignier
+//
+// This file is released under the 3-clause BSD license. See COPYING-BSD.
+
+// Import the diagram and augment the ending time
+loadScicos();
+loadXcosLibs();
+importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.zcos");
+scs_m.props.tf = 30000;
+
+solverName=["BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functional", "Crank-Nicolson"];
+
+for solver=1:5
+
+    // Select the solver (Crank-Nicolson is solver number 8)
+    scs_m.props.tol(6) = solver;
+    if (solver == 5) then scs_m.props.tol(6) = 8; end
+
+    // Set max step size and reltol if Crank-Nicolson
+    if (solver == 5) then scs_m.props.tol(7) = 0.01; scs_m.props.tol(2) = 1d-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 045c711..9f937ae 100644 (file)
@@ -177,6 +177,9 @@ Time for Adams / Functional:
                 <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index 45503d8..abbada9 100644 (file)
@@ -268,6 +268,9 @@ Time for Adams / Functional:
                 <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="IDA">IDA</link>
             </member>
             <member>
index 78bae81..80ce007 100644 (file)
@@ -41,7 +41,7 @@
             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 
+            It is an enhancement of the Euler method, which approximates
             <emphasis>
                 y<subscript>n+1</subscript>
             </emphasis>
             An important difference of <emphasis>Runge-Kutta</emphasis> with the previous methods is that it computes up to the fourth 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 
+            Here, the next value is determined by the present value
             <emphasis>
                 y<subscript>n</subscript>
             </emphasis>
             plus the weighted average of four 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>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</emphasis> is the increment based on the slope at the midpoint of the interval, using 
+                    <emphasis>k2</emphasis> is the increment based on the slope at the midpoint of the interval, using
                     <emphasis>
                         y<subscript>n</subscript> + h*k1/2
                     </emphasis>
                     ,
                 </listitem>
                 <listitem>
-                    <emphasis>k3</emphasis> is again the increment based on the slope at the midpoint, but now using 
+                    <emphasis>k3</emphasis> is again the increment based on the slope at the midpoint, but now using
                     <emphasis>
                         y<subscript>n</subscript> + h*k2/2
                     </emphasis>
                 </listitem>
                 <listitem>
-                    <emphasis>k4</emphasis> is the increment based on the slope at the end of the interval, using 
+                    <emphasis>k4</emphasis> is the increment based on the slope at the end of the interval, using
                     <emphasis>
                         y<subscript>n</subscript> + h*k3
                     </emphasis>
             </itemizedlist>
         </para>
         <para>
-            We can see that with the <emphasis>ki</emphasis>, we progress in the derivatives of 
+            We can see that with the <emphasis>ki</emphasis>, we progress in the derivatives of
             <emphasis>
                 y<subscript>n</subscript>
             </emphasis>
-            . So in <emphasis>k4</emphasis>, we are approximating 
+            . So in <emphasis>k4</emphasis>, we are approximating
             <emphasis>
                 y<superscript>(4)</superscript><subscript>n</subscript>
             </emphasis>
-            , thus making an error in 
+            , thus making an error in
             <emphasis>
                 O(h<superscript>5</superscript>)
             </emphasis>
             .
         </para>
         <para>
-            So the total error is 
+            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 
+            . 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>Runge-Kutta 4(5)</emphasis>, 
+            That error analysis baptized the method <emphasis>Runge-Kutta 4(5)</emphasis>,
             <emphasis>
                 O(h<superscript>5</superscript>)
             </emphasis>
-            per step, 
+            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 
+            Although the solver works fine for <link linkend="Simulatemenu_Menu_entries">max step size</link> up to
             <emphasis>
                 10<superscript>-3</superscript>
             </emphasis>
@@ -234,6 +234,9 @@ Time for Runge-Kutta:
                 <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index edee838..addcc40 100644 (file)
@@ -218,6 +218,9 @@ Time for Dormand-Prince:
                 <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index ec12b40..a386ed4 100644 (file)
                 <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,
+                        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,
+                        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.
+                        y<subscript>n</subscript> + a31*h*k1 + a32*h*k2 + a33*h*k3
                     </emphasis>
+                    .
                 </listitem>
             </itemizedlist>
         </para>
@@ -277,6 +278,9 @@ Time for implicit Runge-Kutta:
                 <link linkend="DoPri">Dormand-Prince 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
diff --git a/scilab/modules/xcos/help/en_US/solvers/5-CrankNicolson.xml b/scilab/modules/xcos/help/en_US/solvers/5-CrankNicolson.xml
new file mode 100644 (file)
index 0000000..3eaaae9
--- /dev/null
@@ -0,0 +1,346 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) Scilab Enterprises - 2016 - 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.1-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="CRANI">
+    
+    <refnamediv>
+        <refname>Crank-Nicolson 2(3)</refname>
+        <refpurpose>
+            <emphasis>Crank-Nicolson</emphasis> is a numerical solver based on the <emphasis>Runge-Kutta</emphasis>
+            scheme 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 role="description">
+        <title>Description</title>
+        <para>
+            <emphasis>Crank-Nicolson</emphasis> is a numerical solver based on the
+            <emphasis>Runge-Kutta</emphasis> scheme 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>
+            This makes the computation times unpredictable.
+            <emphasis>Runge-Kutta</emphasis>-based solvers do not adapt
+            to the complexity of the problem, but guarantee a stable computation time.
+        </para>
+        <para>
+            This method being implicit, 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>
+            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 <emphasis>Crank-Nicolson</emphasis>
+            with the previous methods is that it computes up to the second 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 two 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 midpoint of the interval, using
+                    <emphasis>
+                        y<subscript>n</subscript> + a11*h*k1/2 + a12*h*k2/2
+                    </emphasis>
+                    ,
+                </listitem>
+                <listitem>
+                    <emphasis>k2</emphasis> is the increment based on the slope at the midpoint of the interval, but now using
+                    <emphasis>
+                        y<subscript>n</subscript>
+                    </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 both <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>
+            While computing a new <emphasis>k2</emphasis> only requires one call to the derivative of
+            <emphasis>
+                y<subscript>n</subscript>
+            </emphasis>
+            ,thus making an error in
+            <emphasis>
+                O(h<superscript>2</superscript>)
+            </emphasis>
+            ,
+            <emphasis>k1</emphasis> requires two calls (one for its initial computation and one for the new one).
+            So in <emphasis>k1</emphasis>, we are approximating
+            <emphasis>
+                y<superscript>(2)</superscript><subscript>n</subscript>
+            </emphasis>
+            ,
+            thus making an error in
+            <emphasis>
+                O(h<superscript>3</superscript>)
+            </emphasis>
+            .
+        </para>
+        <para>
+            So the total error is
+            <emphasis>
+                number of steps * O(h<superscript>3</superscript>)
+            </emphasis>
+            .
+            And since <emphasis>number of steps = interval size / h</emphasis> by definition, the total error is in
+            <emphasis>
+                O(h<superscript>2</superscript>)
+            </emphasis>
+            .
+        </para>
+        <para>
+            That error analysis baptized the method <emphasis>Crank-Nicolson 2(3)</emphasis>:
+            <emphasis>
+                O(h<superscript>3</superscript>)
+            </emphasis>
+            per step,
+            <emphasis>
+                O(h<superscript>2</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 role="examples">
+        <title>Examples</title>
+        <para>
+            <link type="scilab" linkend="scilab.xcos/xcos/examples/solvers/ODE_Example.zcos">
+                <inlinemediaobject>
+                    <imageobject>
+                        <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.zcos" valign="middle"/>
+                    </imageobject>
+                </inlinemediaobject>
+            </link>
+            <scilab:image><![CDATA[
+loadScicos();
+loadXcosLibs();
+importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.zcos");
+scs_m.props.tol(2) = 10^-5;
+scs_m.props.tol(6) = 8;
+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 <emphasis>Crank-Nicolson</emphasis>
+            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.zcos");
+      scs_m.props.tf = 5000;
+
+      // Select the solver Crank-Nicolson and set the precision
+      scs_m.props.tol(2) = 10^-10;
+      scs_m.props.tol(6) = 8;
+      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 Crank-Nicolson:");
+      ]]></programlisting>
+        </para>
+        <para>
+            The Scilab console displays:
+            <screen><![CDATA[
+Time for Crank-Nicolson:
+ 8.911
+            ]]></screen>
+        </para>
+        <para>
+            Now, in the following script, we compare the time difference between <emphasis>Crank-Nicolson</emphasis>
+            and <emphasis>CVode</emphasis> by running the example with the five solvers in turn:
+            <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/integCRANI.sce">
+                Open the script
+            </link>
+        </para>
+        <para>
+            <screen><![CDATA[
+Time for BDF / Newton:
+ 18.894
+
+Time for BDF / Functional:
+ 18.382
+
+Time for Adams / Newton:
+ 10.368
+
+Time for Adams / Functional:
+ 9.815
+
+Time for Crank-Nicolson:
+ 6.652
+            ]]></screen>
+        </para>
+        <para>
+            These results show that on a nonstiff problem, for relatively same precision required and forcing the
+            same step size, <emphasis>Crank-Nicolson</emphasis> competes with <emphasis>Adams / Functional</emphasis>.
+        </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 role="see also">
+        <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="DoPri">Dormand-Prince 4(5)</link>
+            </member>
+            <member>
+                <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
+            </member>
+            <member>
+                <link linkend="DDaskr">DDaskr</link>
+            </member>
+            <member>
+                <link linkend="Comparisons">Comparisons</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 role="bibliography">
+        <title>Bibliography</title>
+        <para>
+            Advances in Computational Mathematics, Volume 6, Issue 1, 1996, Pages 207-226 
+            <ulink url="http://link.springer.com/article/10.1007%2FBF02127704">
+                A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type
+            </ulink>
+        </para>
+        <para>
+            Pages 81-82 <ulink url="http://repository.lib.ncsu.edu/ir/bitstream/1840.16/6653/1/etd.pdf">
+                A family of higher-order implicit time integration methods for unsteady compressible flows
+            </ulink>
+        </para>
+        <para>
+            <ulink url="https://computation.llnl.gov/casc/sundials/documentation/documentation.html">Sundials Documentation</ulink>
+        </para>
+    </refsection>
+    
+    <refsection role="history">
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.0.0</revnumber>
+                <revdescription>Crank-Nicolson 2(3) solver added</revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+    
+</refentry>
index f62564a..9a65f07 100644 (file)
@@ -179,6 +179,9 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
                 <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index 7f79c79..22a85ec 100644 (file)
@@ -204,6 +204,9 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
                 <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="Comparisons">Comparisons</link>
             </member>
             <member>
index 7355adc..3d25e3e 100644 (file)
@@ -158,6 +158,9 @@ Time without rootfinding:
                 <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index edd4e97..df4ede7 100644 (file)
@@ -71,6 +71,7 @@
                 <member>Runge-Kutta 4(5)</member>
                 <member>Dormand-Prince 4(5)</member>
                 <member>Implicit Runge-Kutta 4(5)</member>
+                <member>Crank-Nicolson 2(3)</member>
             </simplelist>
         </para>
         <para>
@@ -93,6 +94,7 @@
                 <member>CVode</member>
                 <member>IDA</member>
                 <member>Implicit Runge-Kutta 4(5)</member>
+                <member>Crank-Nicolson 2(3)</member>
             </simplelist>
             <emphasis role="bold">Explicit</emphasis> solvers:
             <simplelist>
@@ -145,7 +147,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
-            In the following script, we compare the time difference between the solvers by running the example with the eight solvers in turn (<emphasis>IDA</emphasis> doesn't handle this kind of problem):
+            In the following script, we compare the time difference between the solvers by running the example with the nine solvers in turn (<emphasis>IDA</emphasis> doesn't handle this kind of problem):
             <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/benchSine.sce">
                 Open the script
             </link>
@@ -176,6 +178,9 @@ Time for Runge-Kutta:
 
 Time for implicit Runge-Kutta:
  10.881
+
+Time for Crank-Nicolson:
+ 7.856
             ]]></screen>
         </para>
         <para>
@@ -212,6 +217,9 @@ Time for implicit Runge-Kutta:
                     <td align="center">
                         <emphasis>Implicit Runge-Kutta</emphasis>
                     </td>
+                    <td align="center">
+                        <emphasis>Crank-Nicolson</emphasis>
+                    </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -224,6 +232,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 1.3x </td>
                     <td align="center"> 0.75x </td>
                     <td align="center"> 1.08x </td>
+                    <td align="center"> .07x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -236,6 +245,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.4x </td>
                     <td align="center"> 0.25x </td>
                     <td align="center"> 0.35x </td>
+                    <td align="center"> 0.23x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -248,6 +258,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.4x </td>
                     <td align="center"> 0.25x </td>
                     <td align="center"> 0.4x </td>
+                    <td align="center"> 0.24x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -260,6 +271,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.75x </td>
                     <td align="center"> 0.45x </td>
                     <td align="center"> 0.6x </td>
+                    <td align="center"> 0.42x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -272,6 +284,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.8x </td>
                     <td align="center"> 0.5x </td>
                     <td align="center"> 0.7x </td>
+                    <td align="center"> 0.45x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -284,6 +297,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> </td>
                     <td align="center"> 0.6x </td>
                     <td align="center"> 0.8x </td>
+                    <td align="center"> 0.56x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -296,6 +310,20 @@ Time for implicit Runge-Kutta:
                     <td align="center"> </td>
                     <td align="center"> </td>
                     <td align="center"> 1.4x </td>
+                    <td align="center"> 0.95x </td>
+                </tr>
+                <tr>
+                    <td align="center">
+                        <emphasis>Implicit Runge-Kutta</emphasis>
+                    </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> 0.67x </td>
                 </tr>
             </informaltable>
         </para>
@@ -317,7 +345,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
-            In the following script, we compare the time difference between the solvers by running the example with the eight solvers in turn (<emphasis>IDA</emphasis> doesn't handle this kind of problem):
+            In the following script, we compare the time difference between the solvers by running the example with the nine solvers in turn (<emphasis>IDA</emphasis> doesn't handle this kind of problem):
             <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/benchBasic.sce">
                 Open the script
             </link>
@@ -348,6 +376,9 @@ Time for Runge-Kutta:
 
 Time for implicit Runge-Kutta:
  5.612
+
+Time for Crank-Nicolson:
+ 3.345
             ]]></screen>
         </para>
         <para>
@@ -381,6 +412,9 @@ Time for implicit Runge-Kutta:
                     <td align="center">
                         <emphasis>Implicit Runge-Kutta</emphasis>
                     </td>
+                    <td align="center">
+                        <emphasis>Crank-Nicolson</emphasis>
+                    </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -393,6 +427,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.2x </td>
                     <td align="center"> 0.17x </td>
                     <td align="center"> 0.5x </td>
+                    <td align="center"> 0.33x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -405,6 +440,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.05x </td>
                     <td align="center"> 0.2x </td>
+                    <td align="center"> 0.12x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -417,6 +453,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.07x </td>
                     <td align="center"> 0.2x </td>
+                    <td align="center"> 0.13x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -429,6 +466,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.15x </td>
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.4x </td>
+                    <td align="center"> 0.22x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -441,6 +479,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.2x </td>
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.5x </td>
+                    <td align="center"> 0.28x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -453,6 +492,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> </td>
                     <td align="center"> 0.7x </td>
                     <td align="center"> 2.4x </td>
+                    <td align="center"> 1.42x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -465,6 +505,20 @@ Time for implicit Runge-Kutta:
                     <td align="center"> </td>
                     <td align="center"> </td>
                     <td align="center"> 3.4x </td>
+                    <td align="center"> 2x </td>
+                </tr>
+                <tr>
+                    <td align="center">
+                        <emphasis>Implicit Runge-Kutta</emphasis>
+                    </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> 0.6x </td>
                 </tr>
             </informaltable>
         </para>
@@ -486,7 +540,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
-            In the following script, we compare the time difference between the solvers by running the example with the eight solvers in turn (<emphasis>IDA</emphasis> doesn't handle this kind of problem):
+            In the following script, we compare the time difference between the solvers by running the example with the nine solvers in turn (<emphasis>IDA</emphasis> doesn't handle this kind of problem):
             <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/benchKalman.sce">
                 Open the script
             </link>
@@ -517,6 +571,9 @@ Time for Runge-Kutta:
 
 Time for implicit Runge-Kutta:
  4
+
+Time for Crank-Nicolson:
+ 2.657
             ]]></screen>
         </para>
         <para>
@@ -550,6 +607,9 @@ Time for implicit Runge-Kutta:
                     <td align="center">
                         <emphasis>Implicit Runge-Kutta</emphasis>
                     </td>
+                    <td align="center">
+                        <emphasis>Crank-Nicolson</emphasis>
+                    </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -562,6 +622,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.4x </td>
+                    <td align="center"> 0.26x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -574,6 +635,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.06x </td>
                     <td align="center"> 0.05x </td>
                     <td align="center"> 0.2x </td>
+                    <td align="center"> 0.12x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -586,6 +648,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.08x </td>
                     <td align="center"> 0.06x </td>
                     <td align="center"> 0.25x </td>
+                    <td align="center"> 0.17x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -598,6 +661,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.07x </td>
                     <td align="center"> 0.3x </td>
+                    <td align="center"> 0.21x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -610,6 +674,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> 0.15x </td>
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.5x </td>
+                    <td align="center"> 0.3x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -622,6 +687,7 @@ Time for implicit Runge-Kutta:
                     <td align="center"> </td>
                     <td align="center"> 0.7x </td>
                     <td align="center"> 3.2x </td>
+                    <td align="center"> 2.1x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -634,6 +700,20 @@ Time for implicit Runge-Kutta:
                     <td align="center"> </td>
                     <td align="center"> </td>
                     <td align="center"> 4.6x </td>
+                    <td align="center"> 3.1x </td>
+                </tr>
+                <tr>
+                    <td align="center">
+                        <emphasis>Implicit Runge-Kutta</emphasis>
+                    </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> 0.66x </td>
                 </tr>
             </informaltable>
         </para>
@@ -807,6 +887,9 @@ Time for DDaskr - GMRes:
                 <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index 542f274..52f9f89 100644 (file)
@@ -177,6 +177,9 @@ Temps pour Adams / Functional :
                 <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index 8d2bd33..5818b58 100644 (file)
@@ -268,6 +268,9 @@ Temps pour Adams / Functional :
                 <link linkend="ImpRK">Runge-Kutta Implicit 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="IDA">IDA</link>
             </member>
             <member>
index e8f843a..1045ab4 100644 (file)
@@ -235,6 +235,9 @@ Temps pour Runge-Kutta :
                 <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index 15271cf..0bea1d3 100644 (file)
@@ -219,6 +219,9 @@ Temps pour Dormand-Prince :
                 <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index f740346..51aa6de 100644 (file)
             <emphasis>
                 10<superscript>-3</superscript>
             </emphasis>
-            , es erreurs d'arrondi surviennent parfois quand l'on approche
+            , les erreurs d'arrondi surviennent parfois quand l'on approche
             <emphasis>
                 4*10<superscript>-4</superscript>
             </emphasis>
@@ -278,6 +278,9 @@ Temps pour Runge-Kutta implicite :
                 <link linkend="DoPri">Dormand-Prince 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
diff --git a/scilab/modules/xcos/help/fr_FR/solvers/5-CrankNicolson.xml b/scilab/modules/xcos/help/fr_FR/solvers/5-CrankNicolson.xml
new file mode 100644 (file)
index 0000000..1d5a2e4
--- /dev/null
@@ -0,0 +1,344 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) Scilab Enterprises - 2016 - 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.1-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="fr_FR" xml:id="CRANI">
+    <refnamediv>
+        <refname>Crank-Nicolson 2(3)</refname>
+        <refpurpose>
+            <emphasis>Crank-Nicolson</emphasis> est un solveur numérique basé sur le schéma <emphasis>Runge-Kutta</emphasis> fournissant une méthode implicite efficace et stable pour résoudre des Problèmes à Valeur Initiale d'Equations Différentielles Ordinaires (EDOs). Appelé par <link linkend="xcos">xcos</link>.
+        </refpurpose>
+    </refnamediv>
+    
+    <refsection role="description">
+        <title>Description</title>
+        <para>
+            <emphasis>Crank-Nicolson</emphasis> est un solveur numérique basé sur le schéma
+            <emphasis>Runge-Kutta</emphasis> fournissant une méthode explicite efficace
+            pour résoudre des Problèmes à Valeur Initiale de la forme :
+        </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> et <emphasis>IDA</emphasis>
+            utilisent un pas variable dans l'intégration.
+        </para>
+        <para>
+            Cela a pour défaut de rendre les temps de calcul imprévisibles.
+            Les solveurs basés sur <emphasis>Runge-Kutta</emphasis> ne s'adaptent
+            pas à la complexité du problème, mais garantissent un temps de calcul stable.
+        </para>
+        <para>
+            Cette méthode étant implicite, elle peut être utilisée sur des problèmes raides.
+        </para>
+        <para>
+            C'est une amélioration de la méthode d'Euler implicite, qui approxime
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+            en calculant
+            <emphasis>
+                f(t<subscript>n</subscript>+h, y<subscript>n+1</subscript>)
+            </emphasis>
+            et en tronquant le développement de Taylor.
+        </para>
+        <para>
+            Par convention, pour utiliser des pas fixes, le programme commence par calculer un pas
+            <emphasis>h</emphasis> qui approche le paramètre de simulation
+            <link linkend="Simulatemenu_Menu_entries">max step size</link>.
+        </para>
+        <para>
+            Une différence notable de <emphasis>Crank-Nicolson</emphasis>
+            par rapport à <emphasis>Sundials</emphasis> est qu'il calcule
+            jusqu'à la dérivée seconde de <emphasis>y</emphasis>,
+            alors que les autres n'utilisent que des combinaisons linéaires de
+            <emphasis>y</emphasis> et <emphasis>y'</emphasis>.
+        </para>
+        <para>
+            Ici, la valeur suivante est déterminée par la valeur actuelle
+            <emphasis>
+                y<subscript>n</subscript>
+            </emphasis>
+            plus la moyenne pondérée de deux increments, où chaque incrément est le produit du pas,
+            <emphasis>h</emphasis>, et une estimation de la pente spécifiée par la fonction
+            <emphasis>f(t,y)</emphasis>. Ils sont distribués de manière à peu près égale sur l'intervalle.
+            <itemizedlist>
+                <listitem>
+                    <emphasis>k1</emphasis> est l'incrément basé sur la pente au milieu de l'intervalle, utilisant
+                    <emphasis>
+                        y<subscript>n</subscript>+ a11*h*k1/2 + a12*h*k2/2
+                    </emphasis>
+                    ,
+                </listitem>
+                <listitem>
+                    <emphasis>k2</emphasis> est l'incrément basé sur la pente au milieu de l'intervalle, mais utilisant
+                    <emphasis>
+                        y<subscript>n</subscript>
+                    </emphasis>
+                    ,
+                </listitem>
+            </itemizedlist>
+        </para>
+        <para>
+            On peut voir que le calcul d'un <emphasis>ki</emphasis>
+            requiert <emphasis>ki</emphasis>, nécéssitant ainsi l'utilisation
+            d'un solveur non linéraire (ici, itérations point-fixes).
+        </para>
+        <para>
+            D'abord, on initialise
+            <emphasis>
+                k0 = h * f(t<subscript>n</subscript>, y<subscript>n</subscript>)
+            </emphasis>
+            comme première estimation pour les deux <emphasis>ki</emphasis>,
+            pour obtenir de nouveaux <emphasis>ki</emphasis> et une première valeur pour
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+            .
+        </para>
+        <para>
+            Ensuite, on sauve les valeurs et on recalcule
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+            avec ces nouveaux <emphasis>ki</emphasis>.
+        </para>
+        <para>
+            Puis on compare les deux
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+            et on le recalcule jusqu'à ce que son écart avec le dernier
+            soit inférieur au paramètre de simulation <emphasis>reltol</emphasis>.
+        </para>
+        <para>
+            Ce processus ajoute un temps de calcul significatif à la méthode,
+            mais améliore grandement la stabilité.
+        </para>
+        <para>
+            Alors que le calcul d'un nouveau <emphasis>k2</emphasis> ne requiert qu'un appel à la dérivée de
+            <emphasis>
+                y<subscript>n</subscript>
+            </emphasis>
+            ,faisant donc une erreur en
+            <emphasis>
+                O(h<superscript>2</superscript>)
+            </emphasis>
+            ,
+            <emphasis>k1</emphasis> requiert deux appels (un pour sa valeur initiale et un pour sa nouvelle valeur).
+            Donc avec <emphasis>k1</emphasis>, on approxime
+            <emphasis>
+                y<superscript>(2)</superscript><subscript>n</subscript>
+            </emphasis>
+            ,
+            faisant donc une erreur en
+            <emphasis>
+                O(h<superscript>3</superscript>)
+            </emphasis>
+            .
+        </para>
+        <para>
+            L'erreur totale est donc
+            <emphasis>
+                nombre de pas * O(h<superscript>3</superscript>)
+            </emphasis>
+            .
+            Et puisque par définition <emphasis>nombre de pas = taille de l'intervalle / h</emphasis>, l'erreur totale est en
+            <emphasis>
+                O(h<superscript>2</superscript>)
+            </emphasis>
+            .
+        </para>
+        <para>
+            Cette analyse d'erreur a baptisé la méthode <emphasis>Crank-Nicolson 2(3)</emphasis>:
+            <emphasis>
+                O(h<superscript>3</superscript>)
+            </emphasis>
+            par pas,
+            <emphasis>
+                O(h<superscript>2</superscript>)
+            </emphasis>
+            au total.
+        </para>
+        <para>
+            Bien que le solveur fonctionne bien pour
+            <link linkend="Simulatemenu_Menu_entries">max step size</link> jusqu'à
+            <emphasis>
+                10<superscript>-3</superscript>
+            </emphasis>
+            ,
+            les erreurs d'arrondi surviennent parfois quand l'on approche
+            <emphasis>
+                4*10<superscript>-4</superscript>
+            </emphasis>
+            .
+            En effet, le scindage de l'intervalle ne peut pas être effectué
+            correctement et l'on obtient des résultats imprévisibles.
+        </para>
+    </refsection>
+    
+    <refsection role="examples">
+        <title>Examples</title>
+        <para>
+            <link type="scilab" linkend="scilab.xcos/xcos/examples/solvers/ODE_Example.zcos">
+                <inlinemediaobject>
+                    <imageobject>
+                        <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.zcos" valign="middle"/>
+                    </imageobject>
+                </inlinemediaobject>
+            </link>
+            <scilab:image><![CDATA[
+loadScicos();
+loadXcosLibs();
+importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.zcos");
+scs_m.props.tol(2) = 10^-5;
+scs_m.props.tol(6) = 8;
+scs_m.props.tol(7) = 10^-2;
+try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
+]]></scilab:image>
+        </para>
+        <para>
+            Le bloc intégrale retourne son état continu, on peut l'évaluer avec <emphasis>Crank-Nicolson</emphasis> en lançant l'exemple :
+        </para>
+        <para>
+            <programlisting language="example"><![CDATA[
+      // Import du diagramme et réglage du temps final
+      loadScicos();
+      loadXcosLibs();
+      importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.zcos");
+      scs_m.props.tf = 5000;
+
+      // Sélection de Crank-Nicolson et réglage de la précision
+      scs_m.props.tol(2) = 10^-10;
+      scs_m.props.tol(6) = 8;
+      scs_m.props.tol(7) = 10^-2;
+
+      // Lancement du chrono, de la simulation et affichage du temps
+      tic();
+      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
+      t = toc();
+      disp(t, "Temps pour Crank-Nicolson :");
+      ]]></programlisting>
+        </para>
+        <para>
+            La console Scilab affiche :
+            <screen><![CDATA[
+Temps pour Crank-Nicolson :
+ 4.615
+            ]]></screen>
+        </para>
+        <para>
+            Maintenant, dans le script suivant, on compare la différence de temps entre <emphasis>Crank-Nicolson</emphasis> et <emphasis>CVode</emphasis> en lançant l'exemple avec différents solveurs tour à tour :
+            <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/integCRANI.sce">
+                Ouverture du script
+            </link>
+        </para>
+        <para>
+            <screen><![CDATA[
+Temps pour BDF / Newton :
+ 18.894
+
+Temps pour BDF / Functional :
+ 18.382
+
+Temps pour Adams / Newton :
+ 10.368
+
+Temps pour Adams / Functional :
+ 9.815
+
+Temps pour Crank-Nicolson :
+ 4.652
+            ]]></screen>
+        </para>
+        <para>
+            Ces résultats montrent que pour un problème non-raide, pour à peu près la même précision demandée et en forçant la même taille de pas, <emphasis>Crank-Nicolson</emphasis> est plus rapide.
+        </para>
+    </refsection>
+    
+    <refsection role="see also">
+        <title>Voir Aussi</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="DoPri">Dormand-Prince 4(5)</link>
+            </member>
+            <member>
+                <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
+            </member>
+            <member>
+                <link linkend="DDaskr">DDaskr</link>
+            </member>
+            <member>
+                <link linkend="Comparaisons">Comparaisons</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 role="bibliography">
+        <title>Bibliographie</title>
+        <para>
+            Advances in Computational Mathematics, Volume 6, Issue 1, 1996,
+            Pages 207-226
+            <ulink url="http://link.springer.com/article/10.1007%2FBF02127704">
+                A practical method for numerical evaluation of
+                solutions of partial differential equations of the heat-conduction type
+            </ulink>
+        </para>
+        <para>
+            Pages 81-82
+            <ulink url="http://repository.lib.ncsu.edu/ir/bitstream/1840.16/6653/1/etd.pdf">
+                A family of higher-order implicit time integration methods for unsteady compressible flows
+            </ulink>
+        </para>
+        <para>
+            <ulink url="https://computation.llnl.gov/casc/sundials/documentation/documentation.html">Documentation Sundials</ulink>
+        </para>
+    </refsection>
+    
+    <refsection role="history">
+        <title>Historique</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.0.0</revnumber>
+                <revdescription>Crank-Nicolson 2(3) ajouté</revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+    
+</refentry>
index f9ed7dd..ae63269 100644 (file)
@@ -179,6 +179,9 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
                 <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index 808f349..0fbcb58 100644 (file)
@@ -204,6 +204,9 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
                 <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="Comparaisons">Comparaisons</link>
             </member>
             <member>
index 060db1f..3b92701 100644 (file)
@@ -159,6 +159,9 @@ Temps sans recherche de racine :
                 <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index 231fc94..1d88085 100644 (file)
@@ -71,6 +71,7 @@
                 <member>Runge-Kutta 4(5)</member>
                 <member>Dormand-Prince 4(5)</member>
                 <member>Runge-Kutta Implicite 4(5)</member>
+                <member>Crank-Nicolson 2(3)</member>
             </simplelist>
         </para>
         <para>
@@ -93,6 +94,7 @@
                 <member>CVode</member>
                 <member>IDA</member>
                 <member>Runge-Kutta Implicite 4(5)</member>
+                <member>Crank-Nicolson 2(3)</member>
             </simplelist>
             <emphasis role="bold">Explicites</emphasis>:
             <simplelist>
@@ -145,7 +147,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
-            Dans le script suivant, on compare la différence de temps entre les solveurs en lançant l'exemple avec les huit solveurs tout à tour (<emphasis>IDA</emphasis> n'est pas adapté à ce genre de problème) :
+            Dans le script suivant, on compare la différence de temps entre les solveurs en lançant l'exemple avec différents solveurs tour à tour (<emphasis>IDA</emphasis> n'est pas adapté à ce genre de problème) :
             <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/benchSine.sce">
                 Ouverture du script
             </link>
@@ -176,6 +178,9 @@ Temps pour Runge-Kutta :
 
 Temps pour Runge-Kutta implicite :
  10.881
+
+Temps pour Crank-Nicolson :
+ 7.856
             ]]></screen>
         </para>
         <para>
@@ -212,6 +217,9 @@ Temps pour Runge-Kutta implicite :
                     <td align="center">
                         <emphasis>Runge-Kutta Implicite</emphasis>
                     </td>
+                    <td align="center">
+                        <emphasis>Crank-Nicolson</emphasis>
+                    </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -224,6 +232,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 1.3x </td>
                     <td align="center"> 0.75x </td>
                     <td align="center"> 1.08x </td>
+                    <td align="center"> .07x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -236,6 +245,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.4x </td>
                     <td align="center"> 0.25x </td>
                     <td align="center"> 0.35x </td>
+                    <td align="center"> 0.23x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -248,6 +258,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.4x </td>
                     <td align="center"> 0.25x </td>
                     <td align="center"> 0.4x </td>
+                    <td align="center"> 0.24x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -260,6 +271,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.75x </td>
                     <td align="center"> 0.45x </td>
                     <td align="center"> 0.6x </td>
+                    <td align="center"> 0.42x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -272,6 +284,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.8x </td>
                     <td align="center"> 0.5x </td>
                     <td align="center"> 0.7x </td>
+                    <td align="center"> 0.45x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -284,6 +297,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> </td>
                     <td align="center"> 0.6x </td>
                     <td align="center"> 0.8x </td>
+                    <td align="center"> 0.56x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -296,6 +310,20 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> </td>
                     <td align="center"> </td>
                     <td align="center"> 1.4x </td>
+                    <td align="center"> 0.95x </td>
+                </tr>
+                <tr>
+                    <td align="center">
+                        <emphasis>Runge-Kutta Implicite</emphasis>
+                    </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> 0.67x </td>
                 </tr>
             </informaltable>
         </para>
@@ -317,7 +345,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
-            Dans le script suivant, on compare la différence de temps entre les solveurs en lançant l'exemple avec les huit solveurs tout à tour (<emphasis>IDA</emphasis> n'est pas adapté à ce genre de problème) :
+            Dans le script suivant, on compare la différence de temps entre les solveurs en lançant l'exemple avec les neuf solveurs tout à tour (<emphasis>IDA</emphasis> n'est pas adapté à ce genre de problème) :
             <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/benchBasic.sce">
                 Ouverture du script
             </link>
@@ -348,6 +376,9 @@ Temps pour Runge-Kutta :
 
 Temps pour Runge-Kutta implicite :
  5.612
+
+Temps pour Crank-Nicolson :
+ 3.345
             ]]></screen>
         </para>
         <para>
@@ -381,6 +412,9 @@ Temps pour Runge-Kutta implicite :
                     <td align="center">
                         <emphasis>Runge-Kutta Implicite</emphasis>
                     </td>
+                    <td align="center">
+                        <emphasis>Crank-Nicolson</emphasis>
+                    </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -393,6 +427,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.2x </td>
                     <td align="center"> 0.17x </td>
                     <td align="center"> 0.5x </td>
+                    <td align="center"> 0.33x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -405,6 +440,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.05x </td>
                     <td align="center"> 0.2x </td>
+                    <td align="center"> 0.12x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -417,6 +453,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.07x </td>
                     <td align="center"> 0.2x </td>
+                    <td align="center"> 0.13x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -429,6 +466,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.15x </td>
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.4x </td>
+                    <td align="center"> 0.22x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -441,6 +479,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.2x </td>
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.5x </td>
+                    <td align="center"> 0.28x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -453,6 +492,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> </td>
                     <td align="center"> 0.7x </td>
                     <td align="center"> 2.4x </td>
+                    <td align="center"> 1.42x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -465,6 +505,20 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> </td>
                     <td align="center"> </td>
                     <td align="center"> 3.4x </td>
+                    <td align="center"> 2x </td>
+                </tr>
+                <tr>
+                    <td align="center">
+                        <emphasis>Implicit Runge-Kutta</emphasis>
+                    </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> 0.6x </td>
                 </tr>
             </informaltable>
         </para>
@@ -486,7 +540,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
 ]]></scilab:image>
         </para>
         <para>
-            Dans le script suivant, on compare la différence de temps entre les solveurs en lançant l'exemple avec les huit solveurs tout à tour (<emphasis>IDA</emphasis> n'est pas adapté à ce genre de problème) :
+            Dans le script suivant, on compare la différence de temps entre les solveurs en lançant l'exemple avec les neuf solveurs tout à tour (<emphasis>IDA</emphasis> n'est pas adapté à ce genre de problème) :
             <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/benchKalman.sce">
                 Ouverture du script
             </link>
@@ -517,6 +571,9 @@ Temps pour Runge-Kutta :
 
 Temps pour Runge-Kutta implicite :
  4
+
+Temps pour Crank-Nicolson :
+ 2.657
             ]]></screen>
         </para>
         <para>
@@ -550,6 +607,9 @@ Temps pour Runge-Kutta implicite :
                     <td align="center">
                         <emphasis>Runge-Kutta Implicite</emphasis>
                     </td>
+                    <td align="center">
+                        <emphasis>Crank-Nicolson</emphasis>
+                    </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -562,6 +622,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.4x </td>
+                    <td align="center"> 0.26x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -574,6 +635,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.06x </td>
                     <td align="center"> 0.05x </td>
                     <td align="center"> 0.2x </td>
+                    <td align="center"> 0.12x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -586,6 +648,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.08x </td>
                     <td align="center"> 0.06x </td>
                     <td align="center"> 0.25x </td>
+                    <td align="center"> 0.17x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -598,6 +661,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.07x </td>
                     <td align="center"> 0.3x </td>
+                    <td align="center"> 0.21x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -610,6 +674,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> 0.15x </td>
                     <td align="center"> 0.1x </td>
                     <td align="center"> 0.5x </td>
+                    <td align="center"> 0.3x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -622,6 +687,7 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> </td>
                     <td align="center"> 0.7x </td>
                     <td align="center"> 3.2x </td>
+                    <td align="center"> 2.1x </td>
                 </tr>
                 <tr>
                     <td align="center">
@@ -634,6 +700,20 @@ Temps pour Runge-Kutta implicite :
                     <td align="center"> </td>
                     <td align="center"> </td>
                     <td align="center"> 4.6x </td>
+                    <td align="center"> 3.1x </td>
+                </tr>
+                <tr>
+                    <td align="center">
+                        <emphasis>Implicit Runge-Kutta</emphasis>
+                    </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> </td>
+                    <td align="center"> 0.66x </td>
                 </tr>
             </informaltable>
         </para>
@@ -807,6 +887,9 @@ Temps pour DDaskr - GMRes :
                 <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
+                <link linkend="CRANI">Crank-Nicolson 2(3)</link>
+            </member>
+            <member>
                 <link linkend="DDaskr">DDaskr</link>
             </member>
             <member>
index b5afea2..a472c62 100644 (file)
@@ -169,7 +169,12 @@ public class SetupDialog extends JDialog {
                              "Method: Fixed step",
                              EnumSet.of(SetupDialog.SolverModifiers.MAX_STEP_SIZE)),
 
-        new SolverDescriptor(7, "Implicit RK45 - Runge-Kutta 4(5)",
+        new SolverDescriptor(7, "Implicit RK45 - Implicit Runge-Kutta 4(5)",
+                             "Method: Fixed step, Nonlinear solver= FIXED-POINT",
+                             EnumSet.of(SetupDialog.SolverModifiers.MAX_STEP_SIZE,
+                                        SetupDialog.SolverModifiers.RELATIVE_TOLERANCE)),
+
+        new SolverDescriptor(8, "CRANI - Crank-Nicolson 2(3)",
                              "Method: Fixed step, Nonlinear solver= FIXED-POINT",
                              EnumSet.of(SetupDialog.SolverModifiers.MAX_STEP_SIZE,
                                         SetupDialog.SolverModifiers.RELATIVE_TOLERANCE)),