Implemented both Runke-Kutta methods.

modsim/ass3/methods.c

@@ -21,34 +21,111 @@ 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, j, intervals = (t1 - t0) / dt;
-    double slope[N], t = t0;
+    double dy[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, slope, params) ) {
+        if( f(t, y0, dy, params) ) {
             puts("Error calculating slope.");
             return -1;
         }
 
-        for( j = 0; j < N; j++ )
-            y1[j] = y1[j] + slope[j] * 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 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

@@ -38,7 +38,7 @@ int f_cos(double t, double *y, double *dy, void *params) {
 int main(void) {
     double y0[1] = {.0}, y1[1];
 
-    if( Euler(.0, 10.0, .0001, 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);