Scicos Sundials: minor fixes 86/14386/11
Paul Bignier [Thu, 10 Apr 2014 08:21:03 +0000 (10:21 +0200)]
 * Reordered and indented the coefficients, it is much clearer this way,
and used float divions for more accuracy.
Coefficients source: http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method

 * Used magic numbers instead of hard-coded constants, mainly to add consistency in the code

Aligned the comments.

Change-Id: I2a27b18c7e948f862f885fa19203bfa7a0249b36

12 files changed:
scilab/modules/scicos/src/scicos_sundials/src/cvode/cvode.c
scilab/modules/xcos/help/en_US/solvers/4-ImplicitRK.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-RKImplicite.xml
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

index e56e025..f33e0f5 100644 (file)
@@ -2492,94 +2492,94 @@ static int CVStepDoPri(CVodeMem cv_mem)
     /* Constant coefficients */
     realtype c2, c3, c4, c5, b1, b3, b4, b5, b6, d1, d3, d4, d5, d6, d7;
     realtype a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a73, a74, a75, a76;
-    a21 =  0.2;
-    c2 =  0.2;
-    a31 =  0.075;
-    c3 =  0.3;
-    a32 =  0.225;
-    c4 =  0.8;
-    a41 =  0.9777777777777778;
-    c5 =  0.8888888888888888;
-    a42 = -3.7333333333333333;
-    b1 =  0.0911458333333333;
-    a43 =  3.5555555555555556;
-    b3 =  0.4492362982929021;
-    a51 =  2.9525986892242036;
-    b4 =  0.6510416666666667;
-    a52 = -11.595793324188386;
-    b5 = -0.3223761792452830;
-    a53 =  9.8228928516994361;
-    b6 =  0.1309523809523810;
-    a54 = -0.2908093278463649;
-    d1 =  0.0899131944444444;
-    a61 =  2.8462752525252525;
-    d3 =  0.4534890685834082;
-    a62 = -10.757575757575758;
-    d4 =  0.6140625;
-    a63 =  8.9064227177434725;
-    d5 = -0.2715123820754717;
-    a64 =  0.2784090909090909;
-    d6 =  0.0890476190476190;
-    a65 = -0.2735313036020583;
-    d7 =  0.025;
-    a71 =  0.0911458333333333;
-    a73 =  0.4492362982929021;
-    a74 =  0.6510416666666667;
-    a75 = -0.3223761792452830;
-    a76 =  0.1309523809523810;
+    a21 =  1. / 5.;
+    a31 =  3. / 40.;
+    a32 =  9. / 40.;
+    a41 =  44. / 45.;
+    a42 = -56. / 15.;
+    a43 =  32. / 9.;
+    a51 =  19372. / 6561.;
+    a52 = -25360. / 2187.;
+    a53 =  64448. / 6561.;
+    a54 = -212. / 729.;
+    a61 =  9017. / 3168.;
+    a62 = -355. / 33.;
+    a63 =  46732. / 5247.;
+    a64 =  49. / 176.;
+    a65 = -5103. / 18656.;
+    a71 =  35. / 384.;
+    a73 =  500. / 1113.;
+    a74 =  125. / 192.;
+    a75 = -2187. / 6784.;
+    a76 =  11. / 84.;
+    b1  =  35. / 384.;
+    b3  =  500. / 1113.;
+    b4  =  125. / 192.;
+    b5  = -2187. / 6784.;
+    b6  =  11. / 84.;
+    c2  =  1. / 5.;
+    c3  =  3. / 10.;
+    c4  =  4. / 5.;
+    c5  =  8. / 9.;
+    d1  =  5179. / 57600.;
+    d3  =  7571. / 16695.;
+    d4  =  393. / 640.;
+    d5  = -92097. / 339200.;
+    d6  =  187. / 2100.;
+    d7  =  1. / 40.;
 
     /* K1, K2 */
-    retval = f(tn, zn[0], zn[2], user_data);                 /* zn[2] = K1 = f(Tn, Yn)                 */
-    N_VLinearSum_Serial (a21 * h, zn[2], 1., zn[0], y);      /* y = K1*a21*h + Yn                      */
-    retval = f(tn + c2 * h, y, zn[3], user_data);            /* zn[3] = K2 = f(Tn+c2*h, Yn + a21*h*K1) */
-    N_VScale(b1, zn[2], y);                                  /* y = b1*K1 + 0*K2                       */
-    N_VScale(d1, zn[2], tempv);                              /* tempv = d1*K1 + 0*K2                   */
+    retval = f(tn, zn[0], zn[2], user_data);             /* zn[2] = K1 = f(Tn, Yn)                 */
+    N_VLinearSum_Serial (a21 * h, zn[2], ONE, zn[0], y); /* y = K1*a21*h + Yn                      */
+    retval = f(tn + c2 * h, y, zn[3], user_data);        /* zn[3] = K2 = f(Tn+c2*h, Yn + a21*h*K1) */
+    N_VScale(b1, zn[2], y);                              /* y = b1*K1 + 0*K2                       */
+    N_VScale(d1, zn[2], tempv);                          /* tempv = d1*K1 + 0*K2                   */
 
     /* K3 */
-    N_VLinearSum_Serial (a31 * h, zn[2], 1., zn[0], ftemp);  /* ftemp = Yn + K1*a31*h                             */
-    N_VLinearSum_Serial (a32 * h, zn[3], 1., ftemp, ftemp);  /* ftemp += K2*a32*h                                 */
+    N_VLinearSum_Serial (a31 * h, zn[2], ONE, zn[0], ftemp); /* ftemp = Yn + K1*a31*h                             */
+    N_VLinearSum_Serial (a32 * h, zn[3], ONE, ftemp, ftemp); /* ftemp += K2*a32*h                                 */
     retval = f(tn + c3 * h, ftemp, zn[4], user_data);        /* zn[4] = K3 = f(Tn+c3*h, Yn + a31*h*K1 + a32*h*K2) */
-    N_VLinearSum_Serial (1., y, b3, zn[4], y);               /* y = b1*K1 + 0*K2 + b3*K3                          */
-    N_VLinearSum_Serial (1., tempv, d3, zn[4], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3                      */
+    N_VLinearSum_Serial (ONE, y, b3, zn[4], y);              /* y = b1*K1 + 0*K2 + b3*K3                          */
+    N_VLinearSum_Serial (ONE, tempv, d3, zn[4], tempv);      /* tempv = d1*K1 + 0*K2 + d3*K3                      */
 
     /* K4 */
-    N_VLinearSum_Serial (a41 * h, zn[2], 1., zn[0], ftemp);  /* ftemp = Yn + K1*a41*h                                        */
-    N_VLinearSum_Serial (a42 * h, zn[3], 1., ftemp, ftemp);  /* ftemp += K2*a42*h                                            */
-    N_VLinearSum_Serial (a43 * h, zn[4], 1., ftemp, ftemp);  /* ftemp += K3*a43*h                                            */
+    N_VLinearSum_Serial (a41 * h, zn[2], ONE, zn[0], ftemp); /* ftemp = Yn + K1*a41*h                                        */
+    N_VLinearSum_Serial (a42 * h, zn[3], ONE, ftemp, ftemp); /* ftemp += K2*a42*h                                            */
+    N_VLinearSum_Serial (a43 * h, zn[4], ONE, ftemp, ftemp); /* ftemp += K3*a43*h                                            */
     retval = f(tn + c4 * h, ftemp, zn[5], user_data);        /* zn[5] = K4 = f(Tn+c4*h, Yn + a41*h*K1 + a42*h*K2 + a43*h*K3) */
-    N_VLinearSum_Serial (1., y, b4, zn[5], y);               /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4                             */
-    N_VLinearSum_Serial (1., tempv, d4, zn[5], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4                         */
+    N_VLinearSum_Serial (ONE, y, b4, zn[5], y);              /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4                             */
+    N_VLinearSum_Serial (ONE, tempv, d4, zn[5], tempv);      /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4                         */
 
     /* K5 */
-    N_VLinearSum_Serial (a51 * h, zn[2], 1., zn[0], ftemp);  /* ftemp = Yn + K1*a51*h                                                   */
-    N_VLinearSum_Serial (a52 * h, zn[3], 1., ftemp, ftemp);  /* ftemp += K2*a52*h                                                       */
-    N_VLinearSum_Serial (a53 * h, zn[4], 1., ftemp, ftemp);  /* ftemp += K3*a53*h                                                       */
-    N_VLinearSum_Serial (a54 * h, zn[5], 1., ftemp, ftemp);  /* ftemp += K4*a54*h                                                       */
+    N_VLinearSum_Serial (a51 * h, zn[2], ONE, zn[0], ftemp); /* ftemp = Yn + K1*a51*h                                                   */
+    N_VLinearSum_Serial (a52 * h, zn[3], ONE, ftemp, ftemp); /* ftemp += K2*a52*h                                                       */
+    N_VLinearSum_Serial (a53 * h, zn[4], ONE, ftemp, ftemp); /* ftemp += K3*a53*h                                                       */
+    N_VLinearSum_Serial (a54 * h, zn[5], ONE, ftemp, ftemp); /* ftemp += K4*a54*h                                                       */
     retval = f(tn + c5 * h, ftemp, zn[6], user_data);        /* zn[6] = K5 = f(Tn+c5*h, Yn + a51*h*K1 + a52*h*K2 + a53*h*K3 + a54*h*K4) */
-    N_VLinearSum_Serial (1., y, b5, zn[6], y);               /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4 + b5*K5                                */
-    N_VLinearSum_Serial (1., tempv, d5, zn[6], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5                            */
+    N_VLinearSum_Serial (ONE, y, b5, zn[6], y);              /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4 + b5*K5                                */
+    N_VLinearSum_Serial (ONE, tempv, d5, zn[6], tempv);      /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5                            */
 
     /* K6 */
-    N_VLinearSum_Serial (a61 * h, zn[2], 1., zn[0], ftemp);  /* ftemp = Yn + K1*a61*h                                                           */
-    N_VLinearSum_Serial (a62 * h, zn[3], 1., ftemp, ftemp);  /* ftemp += K2*a62*h                                                               */
-    N_VLinearSum_Serial (a63 * h, zn[4], 1., ftemp, ftemp);  /* ftemp += K3*a63*h                                                               */
-    N_VLinearSum_Serial (a64 * h, zn[5], 1., ftemp, ftemp);  /* ftemp += K3*a64*h                                                               */
-    N_VLinearSum_Serial (a65 * h, zn[6], 1., ftemp, ftemp);  /* ftemp += K3*a65*h                                                               */
+    N_VLinearSum_Serial (a61 * h, zn[2], ONE, zn[0], ftemp); /* ftemp = Yn + K1*a61*h                                                           */
+    N_VLinearSum_Serial (a62 * h, zn[3], ONE, ftemp, ftemp); /* ftemp += K2*a62*h                                                               */
+    N_VLinearSum_Serial (a63 * h, zn[4], ONE, ftemp, ftemp); /* ftemp += K3*a63*h                                                               */
+    N_VLinearSum_Serial (a64 * h, zn[5], ONE, ftemp, ftemp); /* ftemp += K3*a64*h                                                               */
+    N_VLinearSum_Serial (a65 * h, zn[6], ONE, ftemp, ftemp); /* ftemp += K3*a65*h                                                               */
     retval = f(tn + h, ftemp, zn[7], user_data);             /* zn[7] = K6 = f(Tn+h, Yn + a61*h*K1 + a62*h*K2 + a63*h*K3 + a64*h*K4 + a65*h*K5) */
-    N_VLinearSum_Serial (1., y, b6, zn[7], y);               /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4 + b5*K5 + b6*K6                                */
-    N_VLinearSum_Serial (1., tempv, d6, zn[7], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5 + d6*K6                            */
+    N_VLinearSum_Serial (ONE, y, b6, zn[7], y);              /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4 + b5*K5 + b6*K6                                */
+    N_VLinearSum_Serial (ONE, tempv, d6, zn[7], tempv);      /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5 + d6*K6                            */
 
     /* K7 */
-    N_VLinearSum_Serial (a71 * h, zn[2], 1., zn[0], ftemp);  /* ftemp = Yn + K1*a71*h                                                                    */
-    N_VLinearSum_Serial (a73 * h, zn[4], 1., ftemp, ftemp);  /* ftemp += K3*a73*h                                                                        */
-    N_VLinearSum_Serial (a74 * h, zn[5], 1., ftemp, ftemp);  /* ftemp += K3*a74*h                                                                        */
-    N_VLinearSum_Serial (a75 * h, zn[6], 1., ftemp, ftemp);  /* ftemp += K3*a75*h                                                                        */
-    N_VLinearSum_Serial (a76 * h, zn[7], 1., ftemp, ftemp);  /* ftemp += K3*a76*h                                                                        */
+    N_VLinearSum_Serial (a71 * h, zn[2], ONE, zn[0], ftemp); /* ftemp = Yn + K1*a71*h                                                                    */
+    N_VLinearSum_Serial (a73 * h, zn[4], ONE, ftemp, ftemp); /* ftemp += K3*a73*h                                                                        */
+    N_VLinearSum_Serial (a74 * h, zn[5], ONE, ftemp, ftemp); /* ftemp += K3*a74*h                                                                        */
+    N_VLinearSum_Serial (a75 * h, zn[6], ONE, ftemp, ftemp); /* ftemp += K3*a75*h                                                                        */
+    N_VLinearSum_Serial (a76 * h, zn[7], ONE, ftemp, ftemp); /* ftemp += K3*a76*h                                                                        */
     retval = f(tn + h, ftemp, zn[8], user_data);             /* zn[8] = K7 = f(Tn+h, Yn + a71*h*K1 + 0*h*K2 + a73*h*K3 + a74*h*K4 + a75*h*K5 + a76*h*K6) */
-    N_VLinearSum_Serial (1., tempv, d7, zn[8], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5 + d6*K6 + d7*K7                             */
+    N_VLinearSum_Serial (ONE, tempv, d7, zn[8], tempv);      /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5 + d6*K6 + d7*K7                             */
 
     /* Yn+1 */
-    N_VLinearSum_Serial(1., zn[0], h, y, zn[0]);             /* zn[0] = Yn+1 = Yn + h*y */
+    N_VLinearSum_Serial(ONE, zn[0], h, y, zn[0]); /* zn[0] = Yn+1 = Yn + h*y */
 
     /* Check for errors in the evaluations of f thanks to retval */
     if (retval < 0)
@@ -2623,19 +2623,19 @@ static int CVStepExpRK(CVodeMem cv_mem)
     int retval;
 
     retval = f(tn, zn[0], ftemp, user_data);             /* ftemp = K1                            */
-    N_VLinearSum_Serial (h / 2, ftemp, 1., zn[0], y);    /* y = K1*h/2 + Yn                       */
-    retval = f(tn + h / 2, y, tempv, user_data);         /* tempv = K2 = f(Tn+h/2, Yn + (h/2)*K1) */
-    N_VLinearSum_Serial (1., ftemp, 2., tempv, y);       /* y = K1 + 2K2                          */
+    N_VLinearSum_Serial (h / TWO, ftemp, ONE, zn[0], y); /* y = K1*h/2 + Yn                       */
+    retval = f(tn + h / TWO, y, tempv, user_data);       /* tempv = K2 = f(Tn+h/2, Yn + (h/2)*K1) */
+    N_VLinearSum_Serial (ONE, ftemp, TWO, tempv, y);     /* y = K1 + 2K2                          */
 
-    N_VLinearSum_Serial (h / 2, tempv, 1., zn[0], ftemp); /* ftemp = Yn + K2*h/2                   */
-    retval = f(tn + h / 2, ftemp, tempv, user_data);     /* tempv = K3 = f(Tn+h/2, Yn + (h/2)*K2) */
-    N_VLinearSum_Serial (1., y, 2., tempv, y);           /* y = K1 + 2K2 + 2K3                    */
+    N_VLinearSum_Serial (h / TWO, tempv, ONE, zn[0], ftemp); /* ftemp = Yn + K2*h/2                   */
+    retval = f(tn + h / TWO, ftemp, tempv, user_data);       /* tempv = K3 = f(Tn+h/2, Yn + (h/2)*K2) */
+    N_VLinearSum_Serial (ONE, y, TWO, tempv, y);             /* y = K1 + 2K2 + 2K3                    */
 
-    N_VLinearSum_Serial (h, tempv, 1., zn[0], ftemp);    /* ftemp = Yn + K3*h                 */
-    retval = f(tn + h, ftemp, tempv, user_data);         /* tempv = K4 = f(Tn+h/2, Yn + h*K3) */
-    N_VLinearSum_Serial (1., y, 1., tempv, y);           /* y = K1 + 2K2 + 2K3 + K4           */
+    N_VLinearSum_Serial (h, tempv, ONE, zn[0], ftemp); /* ftemp = Yn + K3*h                 */
+    retval = f(tn + h, ftemp, tempv, user_data);       /* tempv = K4 = f(Tn+h/2, Yn + h*K3) */
+    N_VLinearSum_Serial (ONE, y, ONE, tempv, y);       /* y = K1 + 2K2 + 2K3 + K4           */
 
-    N_VLinearSum_Serial(1., zn[0], h / 6, y, zn[0]);     /* zn[0] = Yn+1 = Yn + y*h/6 */
+    N_VLinearSum_Serial(ONE, zn[0], h / 6., y, zn[0]); /* zn[0] = Yn+1 = Yn + y*h/6 */
 
     /* Check for errors in the evaluations of f thanks to retval */
     if (retval < 0)
@@ -2654,7 +2654,7 @@ static int CVStepExpRK(CVodeMem cv_mem)
     nst++;
 
     /* Update the Nordsieck history array */
-    retval = f(tn, zn[0], zn[1], user_data);  /* zn[1] = y'(tn) */
+    retval = f(tn, zn[0], zn[1], user_data); /* zn[1] = y'(tn) */
 
     return(CV_SUCCESS);
 }
@@ -2686,59 +2686,60 @@ static int CVStepImpRK(CVodeMem cv_mem)
     c2  =  0.434296446908075;
     c3  =  0.758519768667167;
 
-    difference = 0;  /* Difference between two computed solutions */
-    nb_iter    = 1;  /* Iterations counter                        */
-    maxcor     = 30; /* Set maximum number of iterations          */
+    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, 4] to represent the Runge-Kutta coefficients K1, K2, K3.
-     * Set zn[1] = h*y'(tn) as the first guess for the K[i]. */
+     * 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], user_data);        /* K1 = f(tn+c1h, Yn + a11K1) */
+    N_VLinearSum_Serial(ONE, zn[0], h * a11, zn[2], ftemp); /* ftemp = a11hK1 + Yn         */
+    retval = f(tn + c1 * h, ftemp, zn[2], user_data);       /* K1 = f(tn+c1h, Yn + a11hK1) */
 
-    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], user_data);        /* K2 = f(tn+c2h, K2)  */
+    N_VLinearSum_Serial(h * a21, zn[2], h * a22, zn[3], ftemp); /* K2 = a21hK1 + a22hK2 */
+    N_VLinearSum_Serial(ONE, zn[0], ONE, ftemp, ftemp);         /* K2 = Yn + K2         */
+    retval = f(tn + c2 * h, ftemp, zn[3], user_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], user_data);        /* K3 = f(tn+c3h, K3)  */
+    N_VLinearSum_Serial(h * a32, zn[3], h * a33, zn[4], ftemp); /* K3 = a32hK2 + a33hK3 */
+    N_VLinearSum_Serial(h * a31, zn[2], ONE, ftemp, ftemp);     /* K3 = a31hK1 + K3     */
+    N_VLinearSum_Serial(ONE, zn[0], ONE, ftemp, ftemp);         /* K3 = Yn + K3         */
+    retval = f(tn + c3 * h, ftemp, zn[4], user_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 */
+    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 + hK3 */
 
     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], user_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], user_data);         /* K2 = f(tn+c2h, K2)         */
+        N_VLinearSum_Serial(ONE, zn[0], h * a11, zn[2], ftemp);     /* ftemp = a11hK1 + Yn         */
+        retval = f(tn + c1 * h, ftemp, zn[2], user_data);           /* K1 = f(tn+c1h, Yn + a11hK1) */
+        N_VLinearSum_Serial(h * a21, zn[2], h * a22, zn[3], ftemp); /* K2 = a21hK1 + a22hK2        */
+        N_VLinearSum_Serial(ONE, zn[0], ONE, ftemp, ftemp);         /* K2 = Yn + K2                */
+        retval = f(tn + c2 * h, ftemp, zn[3], user_data);           /* K2 = f(tn+c2h, hK2)         */
 
-        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], user_data);         /* K3 = f(tn+c3h, K3)  */
+        N_VLinearSum_Serial(h * a32, zn[3], h * a33, zn[4], ftemp); /* K3 = a32hK2 + a33hK3 */
+        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], user_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 */
+        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 + hK3 */
 
         /* Convergence test */
-        N_VLinearSum_Serial(ONE, tempv, -ONE, y, ftemp);   /* ftemp = tempv-y       */
-        difference = N_VMaxNorm(ftemp);                     /* max = Max(ABS(ftemp)) */
+        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], user_data);  /*                         - zn[1] = Y'(tn) */
-            N_VScale (h, zn[1], zn[1]);               /* Scale zn[1] by h                         */
+            tn += h;                                 /* Increment tn                             */
+            nst++;                                   /* Increment the solver calls               */
+            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 */
index 75efc2e..ec12b40 100644 (file)
@@ -58,7 +58,7 @@
             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>implicit Runge-Kutta</emphasis> 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>.
+            An important difference of <emphasis>implicit 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
index 9c60841..edd4e97 100644 (file)
         <para>
             To put it simply, the <emphasis>Explicit</emphasis> go straight to a computation of the solution, whereas the <emphasis>Implicit</emphasis> concentrate on stability, involving more operations, following the tolerances.
         </para>
-        <title>So how to choose ?</title>
+        <title>So how to choose?</title>
         <para>
             Because it is not possible to know for sure whether a solver will be efficient for a given system or not, the best way is to run the most probable one on it and to compare their results.
         </para>
                 <emphasis role="bold">Precision:</emphasis> <emphasis>CVode</emphasis>,
             </para>
             <para>
-                <emphasis role="bold">Predictable time:</emphasis> Fixed-size step.
+                <emphasis role="bold">Predictable time:</emphasis> Fixed-size step,
             </para>
             <para>
                 <emphasis role="bold">Simulation time:</emphasis> <emphasis>LSodar</emphasis>,
             </para>
             <para>
-                <emphasis role="bold">Automated:</emphasis> <emphasis>LSodar</emphasis>,
+                <emphasis role="bold">Automated:</emphasis> <emphasis>LSodar</emphasis>.
             </para>
         </para>
     </refsection>
index 939c9f4..542f274 100644 (file)
@@ -174,7 +174,7 @@ Temps pour Adams / Functional :
                 <link linkend="DoPri">Dormand-Prince 4(5)</link>
             </member>
             <member>
-                <link linkend="RKImp">Runge-Kutta Implicite 4(5)</link>
+                <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
                 <link linkend="DDaskr">DDaskr</link>
index fadb1d9..8d2bd33 100644 (file)
@@ -69,7 +69,7 @@
             </para>
         </para>
         <para>
-            Ces méthodes implicites sont caractérisées par leur order <emphasis>q</emphasis>, qui indique le nombre de points intermédiaires requis pour le calcul de
+            Ces méthodes implicites sont caractérisées par leur ordre <emphasis>q</emphasis>, qui indique le nombre de points intermédiaires requis pour le calcul de
             <emphasis>
                 y<subscript>n+1</subscript>
             </emphasis>
             </emphasis>
         </para>
         <para>
-            La fonction est appelée entre les activations, parce-qu'une activation est susceptible de modifier le système.
+            Le solveur est appelé entre les activations, parce-qu'une activation est susceptible de modifier le système.
         </para>
         <para>
             Suivant la criticalité de l'événement (son effet sur le problème continu), soit le solveur poursuit avec temps initial et final différents comme si rien ne s'était passé, soit, si le système a été modifié, il faut "réinitialiser à froid" le solveur et le relancer.
@@ -265,7 +265,7 @@ Temps pour Adams / Functional :
                 <link linkend="DoPri">Dormand-Prince 4(5)</link>
             </member>
             <member>
-                <link linkend="RKImp">Runge-Kutta Implicit 4(5)</link>
+                <link linkend="ImpRK">Runge-Kutta Implicit 4(5)</link>
             </member>
             <member>
                 <link linkend="IDA">IDA</link>
index eca044c..e8f843a 100644 (file)
@@ -41,7 +41,7 @@
             Cette méthode est explicite, donc non concernée par les itérations fonctionnelles ou de Newton, et est déconseillée pour les problèmes raides.
         </para>
         <para>
-            C'est une amélioration de la méthode d'Euler, qui approxime 
+            C'est une amélioration de la méthode d'Euler, qui approxime
             <emphasis>
                 y<subscript>n+1</subscript>
             </emphasis>
             Une différence notable de <emphasis>Runge-Kutta</emphasis> par rapport à <emphasis>Sundials</emphasis> est qu'il calcule jusqu'à la dérivée quatrième 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 
+            Ici, la valeur suivante
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+             est déterminée par la valeur actuelle
             <emphasis>
                 y<subscript>n</subscript>
             </emphasis>
             plus la moyenne pondérée de quatre incréments, 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> :
             <itemizedlist>
                 <listitem>
-                    <emphasis>k1</emphasis> est l'incrément basé sur la pente au début de l'intervalle, utilisant 
+                    <emphasis>k1</emphasis> est l'incrément basé sur la pente au début de l'intervalle, utilisant
                     <emphasis>
                         y<subscript>n</subscript>
                     </emphasis>
                     (méthode d'Euler),
                 </listitem>
                 <listitem>
-                    <emphasis>k2</emphasis> est l'incrément basé sur la pente au milieu de l'intervalle, utilisant 
+                    <emphasis>k2</emphasis> est l'incrément basé sur la pente au milieu de l'intervalle, utilisant
                     <emphasis>
                         y<subscript>n</subscript> + h*k1/2
                     </emphasis>
                     ,
                 </listitem>
                 <listitem>
-                    <emphasis>k3</emphasis> est à nouveau l'incrément basé sur la pente au milieu de l'intervalle, utilisant 
+                    <emphasis>k3</emphasis> est à nouveau l'incrément basé sur la pente au milieu de l'intervalle, utilisant
                     <emphasis>
                         y<subscript>n</subscript> + h*k2/2
                     </emphasis>
                 </listitem>
                 <listitem>
-                    <emphasis>k4</emphasis> est l'incrément basé sur la pente à la fin de l'intervalle, utilisant 
+                    <emphasis>k4</emphasis> est l'incrément basé sur la pente à la fin de l'intervalle, utilisant
                     <emphasis>
                         y<subscript>n</subscript> + h*k3
                     </emphasis>
             </itemizedlist>
         </para>
         <para>
-            On peut voir qu'avec les <emphasis>ki</emphasis>, on progresse dans les dérivées de 
+            On peut voir qu'avec les <emphasis>ki</emphasis>, on progresse dans les dérivées de
             <emphasis>
                 y<subscript>n</subscript>
             </emphasis>
-            . Donc à <emphasis>k4</emphasis>, on approxime 
+            . Donc à <emphasis>k4</emphasis>, on approxime
             <emphasis>
                 y<superscript>(4)</superscript><subscript>n</subscript>
             </emphasis>
-            , faisant donc une erreur en 
+            , faisant donc une erreur en
             <emphasis>
                 O(h<superscript>5</superscript>)
             </emphasis>
             .
         </para>
         <para>
-            L'erreur totale est donc 
+            L'erreur totale est donc
             <emphasis>
                 nombre de pas * O(h<superscript>5</superscript>)
             </emphasis>
-            . Et puisque par définition <emphasis>nombre de pas = taille de l'intervalle / h</emphasis>, l'erreur totale est en 
+            . Et puisque par définition <emphasis>nombre de pas = taille de l'intervalle / h</emphasis>, l'erreur totale est en
             <emphasis>
                 O(h<superscript>4</superscript>)
             </emphasis>
             .
         </para>
         <para>
-            Cette analyse d'erreur a baptisé la méthode <emphasis>Runge-Kutta 4(5)</emphasis>, 
+            Cette analyse d'erreur a baptisé la méthode "<emphasis>Runge-Kutta 4(5)</emphasis>",
             <emphasis>
                 O(h<superscript>5</superscript>)
             </emphasis>
-            par pas, 
+            par pas,
             <emphasis>
                 O(h<superscript>4</superscript>)
             </emphasis>
             au total.
         </para>
         <para>
-            Bien que le solveur fonctionne bien pour <link linkend="Simulatemenu_Menu_entries">max step size</link> allant juqu'à 
+            Bien que le solveur fonctionne bien pour <link linkend="Simulatemenu_Menu_entries">max step size</link> allant juqu'à
             <emphasis>
                 10<superscript>-3</superscript>
             </emphasis>
@@ -228,7 +232,7 @@ Temps pour Runge-Kutta :
                 <link linkend="DoPri">Dormand-Prince 4(5)</link>
             </member>
             <member>
-                <link linkend="RKImp">Runge-Kutta Implicite 4(5)</link>
+                <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
                 <link linkend="DDaskr">DDaskr</link>
index ca06654..15271cf 100644 (file)
             Une différence notable de <emphasis>Runge-Kutta</emphasis> par rapport à <emphasis>Sundials</emphasis> est qu'il calcule jusqu'à la septième dérivée 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
+            Ici, la valeur suivante
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+             est déterminée par la valeur actuelle
             <emphasis>
                 y<subscript>n</subscript>
             </emphasis>
@@ -212,7 +216,7 @@ Temps pour Dormand-Prince :
                 <link linkend="RK">Runge-Kutta 4(5)</link>
             </member>
             <member>
-                <link linkend="RKImp">Runge-Kutta Implicite 4(5)</link>
+                <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
                 <link linkend="DDaskr">DDaskr</link>
index 384e2a3..9d8b5fa 100644 (file)
@@ -12,7 +12,7 @@
  * For more information, see the COPYING file which you should have received
  * along with this program.
  -->
-<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="RKImp">
+<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="ImpRK">
     <refnamediv>
         <refname>Implicit Runge-Kutta 4(5)</refname>
         <refpurpose>
             Une différence notable de <emphasis>Runge-Kutta implicite</emphasis> par rapport à <emphasis>Sundials</emphasis> est qu'il calcule jusqu'à la dérivée quatrième 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
+            Ici, la valeur suivante
+            <emphasis>
+                y<subscript>n+1</subscript>
+            </emphasis>
+             est déterminée par la valeur actuelle
             <emphasis>
                 y<subscript>n</subscript>
             </emphasis>
             <emphasis>
                 a<subscript>ij</subscript>
             </emphasis>
-            ) nous fait gagner un order, produisant ainsi une erreur par pas de temps en
+            ) nous fait gagner un ordre, produisant ainsi une erreur par pas de temps en
             <emphasis>
                 O(h<superscript>5</superscript>)
             </emphasis>
index dd794d5..f9ed7dd 100644 (file)
@@ -176,7 +176,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
                 <link linkend="DoPri">Dormand-Prince 4(5)</link>
             </member>
             <member>
-                <link linkend="RKImp">Runge-Kutta Implicite 4(5)</link>
+                <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
                 <link linkend="DDaskr">DDaskr</link>
index e108d0b..808f349 100644 (file)
@@ -201,7 +201,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
                 <link linkend="DoPri">Dormand-Prince 4(5)</link>
             </member>
             <member>
-                <link linkend="RKImp">Runge-Kutta Implicite 4(5)</link>
+                <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
                 <link linkend="Comparaisons">Comparaisons</link>
index e9f533f..060db1f 100644 (file)
@@ -156,7 +156,7 @@ Temps sans recherche de racine :
                 <link linkend="DoPri">Dormand-Prince 4(5)</link>
             </member>
             <member>
-                <link linkend="RKImp">Runge-Kutta Implicite 4(5)</link>
+                <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
                 <link linkend="DDaskr">DDaskr</link>
index 231240b..231fc94 100644 (file)
@@ -804,7 +804,7 @@ Temps pour DDaskr - GMRes :
                 <link linkend="DoPri">Dormand-Prince 4(5)</link>
             </member>
             <member>
-                <link linkend="RKImp">Runge-Kutta Implicite 4(5)</link>
+                <link linkend="ImpRK">Runge-Kutta Implicite 4(5)</link>
             </member>
             <member>
                 <link linkend="DDaskr">DDaskr</link>