1 #include "RungeKutta.hpp"
2
3 #include <iomanip>
4
5 using namespace std;
6
7 template <class E>
RungeKutta(int dim)8 RungeKutta<E>::RungeKutta(int dim)
9 {
10 m = dim;
11 k1 = new double[m];
12 k2 = new double[m];
13 k3 = new double[m];
14 k4 = new double[m];
15 y = new double[m];
16 y_tmp = new double[m];
17 }
18
19 template <class E>
~RungeKutta()20 RungeKutta<E>::~RungeKutta()
21 {
22 delete [] k1;
23 delete [] k2;
24 delete [] k3;
25 delete [] k4;
26 delete [] y;
27 delete [] y_tmp;
28 }
29
30 template <class E>
set(E * tmp,double * y0,double beg,double end)31 void RungeKutta<E>::set ( E * tmp , double * y0 , double beg , double end )
32 {
33 unit=tmp;
34 x0=beg; xn=end;
35 x=x0;
36 h=double(xn-x0)/double(N_INTER);
37 for (i=0;i<m;i++) {y[i]=y0[i];}
38 success=true;
39 }
40
41 template <class E>
run()42 bool RungeKutta<E>::run() {
43 for(j=0;j<MAX_ITER_RK;j++) {
44 //Avoid going out of x interval
45 if (x+h >xn) {
46 h = xn-x;
47 j = MAX_ITER_RK;
48 }
49
50 //Compute k1, k2, k3, k4
51 for(i=0;i<m;i++)
52 k1[i] = h*unit->f(i, x, y);
53 for(i=0;i<m;i++)
54 y_tmp[i] = y[i]+k1[i]/2.0;
55 for(i=0;i<m;i++)
56 k2[i] = h*unit->f(i, x+h/2.0, y_tmp);
57 for(i=0;i<m;i++)
58 y_tmp[i] = y[i]+k2[i]/2.0;
59 for(i=0;i<m;i++)
60 k3[i] = h*unit->f(i, x+h/2.0, y_tmp);
61 for(i=0;i<m;i++)
62 y_tmp[i] = y[i]+k3[i];
63 for ( i = 0 ; i < m ; i++ )
64 k4[i] = h*unit->f ( i , x+h , y_tmp );
65 //Compute the new y
66 for(i=0;i<m;i++)
67 y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
68 x += h;
69 }
70
71 if ( x < xn-EPS ) {// MODIF SEB (le EPS)
72 success=false;
73
74 // cout.setf(ios::fixed);
75 // cout << setprecision(12);
76 // cout << "x=" << x << " < xn=" << xn << " diff=" << xn-x << endl;
77 }
78
79 return success;
80 }
81