Merge branch 'master' of vo20.nl:/git/uva

modsim/ass3/methods.c

@@ -20,28 +20,112 @@ occurs (overflow, NaN or such)
 
 int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
         f_ptr f, void *params) {
-    int i;
-    double slope[N];
+    int i, j, intervals = (t1 - t0) / dt;
+    double dy[N], t = t0;
 
-    if( f(t0, y0, slope, params) ) {
-        puts("Error calculating slope.");
-        return -1;
-    }
+    // Initialize with first values (specified)
+    for( j = 0; j < N; j++ )
+        y1[j] = y0[j];
+
+    // Estimate for each interval
+    for( i = 0; i < intervals; i++ ) {
+        if( f(t, y0, dy, params) ) {
+            puts("Error calculating slope.");
+            return -1;
+        }
 
-    for( i = 0; i < N; i++ )
-        y1[i] = y0[i] + slope[i] * dt;
+        for( j = 0; j < N; j++ ) {
+            y1[j] = y1[j] + dy[j] * dt;
+            // Value can be printed/logged here
+        }
+
+        t = t0 + i * dt;
+    }
 
-    return (t1 - t0) > dt ? Euler(t0 + dt, t1, dt, y1, y1, N, f, params) : 0;
+    return 0;
 }
-/*
+
 int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N,
         f_ptr f, void *params) {
+    int i, j, intervals = (t1 - t0) / dt;
+    double dy1[N], dy2[N], t = t0;
+
+    // Initialize with first values (specified)
+    for( j = 0; j < N; j++ )
+        y1[j] = y0[j];
+
+    // Estimate for each interval
+    for( i = 0; i < intervals; i++ ) {
+        if( f(t, y0, dy1, params) || f(t + dt, y1, dy2, params) ) {
+            puts("Error calculating slope.");
+            return -1;
+        }
+
+        for( j = 0; j < N; j++ ) {
+            y1[j] = y1[j] + .5 * (dy1[j] + dy2[j]) * dt;
+            // Value can be printed/logged here
+        }
+
+        t = t0 + i * dt;
+    }
 
     return 0;
 }
 
 int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
-        f_ptr f, void *params)
+        f_ptr f, void *params) {
+    int i, j, intervals = (t1 - t0) / dt;
+    double y[N], dy[N], e1[N], e2[N], e3[N], e4, t = t0, half_dt = .5 * dt;
+
+    // Initialize with first values (specified)
+    for( j = 0; j < N; j++ )
+        y[j] = y1[j] = y0[j];
+
+    // Estimate for each interval
+    for( i = 0; i < intervals; i++ ) {
+        if( f(t, y1, dy, params) ) {
+            puts("Error calculating first estimate.");
+            return -1;
+        }
+
+        for( j = 0; j < N; j++ ) {
+            e1[j] = dy[j] * dt;
+            y[j] = y1[j] + .5 * e1[j];
+        }
+
+        if( f(t + half_dt, y, dy, params) ) {
+            puts("Error calculating second estimate.");
+            return -1;
+        }
+
+        for( j = 0; j < N; j++ ) {
+            e2[j] = dy[j] * dt;
+            y[j] = y1[j] + .5 * e2[j];
+        }
+
+        if( f(t + half_dt, y, dy, params) ) {
+            puts("Error calculating third estimate.");
+            return -1;
+        }
+
+        for( j = 0; j < N; j++ ) {
+            e3[j] = dy[j] * dt;
+            y[j] = y0[j] + e3[j];
+        }
+
+        if( f(t + dt, y, dy, params) ) {
+            puts("Error calculating fourth estimate.");
+            return -1;
+        }
+
+        for( j = 0; j < N; j++ ) {
+            e4 = dy[j] * dt;
+            y1[j] = y1[j] + (e1[j] + 2 * (e2[j] + e3[j]) + e4) / 6;
+            // Value can be printed/logged here
+        }
+
+        t = t0 + i * dt;
+    }
 
     return 0;
-}*/
+}

modsim/ass3/methods.h

@@ -2,10 +2,8 @@ typedef int (*f_ptr)(double t, double *y, double *dy, void *params);
 
 int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
     f_ptr f, void *params);
-/*
 int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N,
     f_ptr f, void *params);
 
 int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
     f_ptr f, void *params);
-    */

modsim/ass3/test.c

@@ -36,9 +36,9 @@ int f_cos(double t, double *y, double *dy, void *params) {
 }
 
 int main(void) {
-    double y0[1] = {.506}, y1[1];
+    double y0[1] = {.0}, y1[1];
 
-    if( Euler(-100.0, 10.0, 0.001, y0, y1, 1, &f_cos, NULL) ) {
+    if( RungeKutta4(.0, 10.0, .000001, y0, y1, 1, &f_cos, NULL) ) {
         puts("FAIL");
     } else {
         printf("y1: %f\n", *y1);