1 #include <stdio.h>
2 #include <time.h>
3 #include <math.h>
4
lorenz(const double * x,double * restrict y)5 void lorenz(const double *x, double *restrict y) {
6 y[0] = 10.0 * (x[1] - x[0]);
7 y[1] = 28.0 * x[0] - x[1] - x[0] * x[2];
8 y[2] = x[0] * x[1] - (8.0 / 3.0) * x[2];
9 }
10
main(int argc,const char * argv[])11 int main(int argc, const char *argv[])
12 {
13 const int nb_steps = 20000000;
14 const double h = 1.0e-10;
15 const double h2 = 0.5 * h;
16 const double nb_loops = 21;
17 double x[3];
18 double y[3];
19 double f1[3];
20 double f2[3];
21 double f3[3];
22 double f4[3];
23 double min_time = 1E6;
24 clock_t begin, end;
25 double time_spent;
26
27 for (int j = 0; j < nb_loops; j++) {
28 x[0] = 8.5;
29 x[1] = 3.1;
30 x[2] = 1.2;
31 begin = clock();
32 for (int k = 0; k < nb_steps; k++) {
33 lorenz(x, f1);
34 for (int i = 0; i < 3; i++) {
35 y[i] = x[i] + h2 * f1[i];
36 }
37 lorenz(y, f2);
38 for (int i = 0; i < 3; i++) {
39 y[i] = x[i] + h2 * f2[i];
40 }
41 lorenz(y, f3);
42 for (int i = 0; i < 3; i++) {
43 y[i] = x[i] + h * f3[i];
44 }
45 lorenz(y, f4);
46 for (int i = 0; i < 3; i++) {
47 x[i] = x[i] + h * (f1[i] + 2 * (f2[i] + f3[i]) + f4[i]) / 6.0;
48 }
49 }
50 end = clock();
51 min_time = fmin(min_time, (double)(end-begin)/CLOCKS_PER_SEC);
52 printf("Result: %f\t runtime: %f\n", x[0], (double)(end-begin)/CLOCKS_PER_SEC);
53 }
54 printf("Minimal Runtime: %f\n", min_time);
55
56 return 0;
57 }
58