Finished integration implementations and test file.

modsim/ass3/methods.c

@@ -29,7 +29,7 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
 
     // Estimate for each interval
     for( i = 0; i < intervals; i++ ) {
-        if( f(t, y0, dy, params) ) {
+        if( f(t, y1, dy, params) ) {
             puts("Error calculating slope.");
             return -1;
         }
@@ -56,7 +56,7 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N,
 
     // Estimate for each interval
     for( i = 0; i < intervals; i++ ) {
-        if( f(t, y0, dy1, params) || f(t + dt, y1, dy2, params) ) {
+        if( f(t, y1, dy1, params) || f(t + dt, y1, dy2, params) ) {
             puts("Error calculating slope.");
             return -1;
         }
@@ -93,7 +93,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
             y[j] = y1[j] + .5 * e1[j];
         }
 
-        if( f(t + half_dt, y, dy, params) ) {
+        if( f(t + half_dt, y1, dy, params) ) {
             puts("Error calculating second estimate.");
             return -1;
         }
@@ -103,7 +103,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
             y[j] = y1[j] + .5 * e2[j];
         }
 
-        if( f(t + half_dt, y, dy, params) ) {
+        if( f(t + half_dt, y1, dy, params) ) {
             puts("Error calculating third estimate.");
             return -1;
         }
@@ -113,7 +113,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
             y[j] = y0[j] + e3[j];
         }
 
-        if( f(t + dt, y, dy, params) ) {
+        if( f(t + dt, y1, dy, params) ) {
             puts("Error calculating fourth estimate.");
             return -1;
         }

modsim/ass3/test.c

@@ -3,12 +3,20 @@
 #include <math.h>
 #include "methods.h"
 
+#define DT .001
+
 int f_test_1(double t, double *y, double *dy, void *params) {
     *dy = 1.0;
 
     return 0;
 }
 
+int f_test_2(double t, double *y, double *dy, void *params) {
+    *dy = t;
+
+    return 0;
+}
+
 int f_test_3(double t, double *y, double *dy, void *params) {
     if( isnan(*y) ) {
         return -1;
@@ -29,20 +37,31 @@ int f_test_4(double t, double *y, double *dy, void *params) {
     }
 }
 
-int f_cos(double t, double *y, double *dy, void *params) {
-    *dy = cos(t);
-
-    return 0;
+#define INTEGRATE(method, y_start, t0, t1, dy) { \
+    *y0 = y_start; \
+    if( method(t0, t1, DT, y0, y1, 1, &dy, NULL) ) \
+        puts("FAIL"); \
+    else \
+        printf("y1 for dy = " #dy " using %-12s %.10f\n", #method ":", *y1); \
+}
+#define COMPARE(y_start, t0, t1, dy) { \
+    INTEGRATE(Euler, y_start, t0, t1, dy); \
+    INTEGRATE(RungeKutta2, y_start, t0, t1, dy); \
+    INTEGRATE(RungeKutta4, y_start, t0, t1, dy); \
 }
 
 int main(void) {
-    double y0[1] = {.0}, y1[1];
+    double y0[1] = {.0}, y1[1] = {.0};
 
-    if( RungeKutta4(.0, 10.0, .000001, y0, y1, 1, &f_cos, NULL) ) {
-        puts("FAIL");
-    } else {
-        printf("y1: %f\n", *y1);
-    }
+    COMPARE(.0, .0, 10.0, f_test_1);
+    puts("");
+    COMPARE(.0, .0, 10.0, f_test_2);
+    puts("");
+    COMPARE(1.0, .0, 5.0, f_test_3);
+    puts("");
+    COMPARE(1.0, -1.0, 1.0, f_test_4);
+    puts("");
+    COMPARE(-1.0, 1.0, 10.0, f_test_4);
 
     return 0;
 }