Bug #10574 fixed - Runge-Kutta-Method failed for vector [x, 1] with x > 9.
[scilab.git] / scilab / modules / differential_equations / src / c / rk4.c
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2007 - INRIA
4  * ...
5  * 
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.  The terms
9  * are also available at    
10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11  *
12  */
13 /*--------------------------------------------------------------------------*/
14 #include "rk4.h"
15 #include "stack-c.h"
16 /*--------------------------------------------------------------------------*/
17 /* Runge-Kutta (RK4) method */
18 /* http://media4.obspm.fr/public/DEA/cours/chapitre3/souschapitre2/section4/page5/section3_2_4_5.html */
19 /*     The original version has been modified to replace statically
20  *    allocated arrays yt, dym and dyt by rwork arguments parts
21  *    array + blas use. Serge Steer INRIA- feb 2012*/
22
23 /*--------------------------------------------------------------------------*/
24 int C2F(rk4)(double *y, double *dydx, int *n,double *x, double *h, double *yout,void (*derivs)(), double *rwork)
25 {
26     double d = 0.0;
27     int i;
28     double h6 = 0.0, hh = 0.0, xh = 0.0;
29
30     double *yt = rwork;
31     double *dym = rwork+*n;
32     double *dyt = rwork+2*(*n);
33  
34     C2F(ierode).iero = 0;
35     hh = *h * 0.5;
36     h6 = *h / 6.0;
37     xh = *x + hh;
38     for (i = 0; i < *n; ++i) yt[i] = y[i] + hh * dydx[i];
39     (*derivs)(n, &xh, yt, dyt);
40
41     if (C2F(ierode).iero > 0) return 0;
42
43     for (i = 0; i < *n; ++i) yt[i] = y[i] + hh * dyt[i];
44     (*derivs)(n, &xh, yt, dym);
45
46     if (C2F(ierode).iero > 0) return 0;
47
48     for (i = 0; i < *n; ++i) 
49     {
50                 yt[i] = y[i] + *h * dym[i];
51                 dym[i] = dyt[i] + dym[i];
52     }
53     d = *x + *h;
54     (*derivs)(n, &d, yt, dyt);
55
56     if (C2F(ierode).iero > 0) return 0;
57
58     for (i = 0; i < *n; ++i) yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + dym[i] * 2.0);
59     return 0;
60 }
61 /*--------------------------------------------------------------------------*/