ModSim ass3: Changed integration methods.

modsim/ass3/compare_osc.sh

+./main 0
+cat report/osc_euler | ./plot.py report/osc_euler.pdf
+cat report/osc_rk2 | ./plot.py report/osc_rk2.pdf
+cat report/osc_rk4 | ./plot.py report/osc_rk4.pdf
+rm report/osc_euler report/osc_rk2 report/osc_rk4

modsim/ass3/integration.c

 #include <stdlib.h>
 #include <stdio.h>
+#include <math.h>
 #include "integration.h"
 #include "logger.h"
 
@@ -21,29 +22,27 @@ 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 dy[N], t = t0;
+    int i;
+    double dy[N], t;
 
     // Initialize with first values (specified)
-    for( j = 0; j < N; j++ )
-        y1[j] = y0[j];
+    for( i = 0; i < N; i++ )
+        y1[i] = y0[i];
 
     logger(t0, N, y0);
 
     // Estimate for each interval
-    for( i = 0; i < intervals; i++ ) {
+    for( t = t0; t < t1; t = t + dt > t1 ? t1 : t + dt ) {
         if( f(t, y1, dy, params) ) {
             puts("Error calculating slope.");
             return -1;
         }
 
-        for( j = 0; j < N; j++ )
-            y1[j] = y1[j] + dy[j] * dt;
+        for( i = 0; i < N; i++ )
+            y1[i] = y1[i] + dy[i] * dt;
 
         // Value can be printed/logged here
         logger(t, N, y1);
-
-        t = t0 + i * dt;
     }
 
     return 0;
@@ -51,29 +50,27 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
 
 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;
+    int i;
+    double dy1[N], dy2[N], t;
 
     // Initialize with first values (specified)
-    for( j = 0; j < N; j++ )
-        y1[j] = y0[j];
+    for( i = 0; i < N; i++ )
+        y1[i] = y0[i];
 
     logger(t0, N, y0);
 
     // Estimate for each interval
-    for( i = 0; i < intervals; i++ ) {
+    for( t = t0; t < t1; t = t + dt > t1 ? t1 : t + dt ) {
         if( f(t, y1, 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;
+        for( i = 0; i < N; i++ )
+            y1[i] = y1[i] + .5 * (dy1[i] + dy2[i]) * dt;
 
         // Value can be printed/logged here
         logger(t, N, y1);
-
-        t = t0 + i * dt;
     }
 
     return 0;
@@ -81,25 +78,25 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N,
 
 int RungeKutta4(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 y[N], dy[N], e1[N], e2[N], e3[N], e4, t = t0, half_dt = .5 * dt;
+    int i;
+    double y[N], dy[N], e1[N], e2[N], e3[N], e4, t, half_dt = .5 * dt;
 
     // Initialize with first values (specified)
-    for( j = 0; j < N; j++ )
-        y[j] = y1[j] = y0[j];
+    for( i = 0; i < N; i++ )
+        y1[i] = y0[i];
 
     logger(t0, N, y0);
 
     // Estimate for each interval
-    for( i = 0; i < intervals; i++ ) {
+    for( t = t0; t < t1; t = t + dt > t1 ? t1 : t + dt ) {
         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];
+        for( i = 0; i < N; i++ ) {
+            e1[i] = dy[i] * dt;
+            y[i] = y1[i] + .5 * e1[i];
         }
 
         if( f(t + half_dt, y1, dy, params) ) {
@@ -107,9 +104,9 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
             return -1;
         }
 
-        for( j = 0; j < N; j++ ) {
-            e2[j] = dy[j] * dt;
-            y[j] = y1[j] + .5 * e2[j];
+        for( i = 0; i < N; i++ ) {
+            e2[i] = dy[i] * dt;
+            y[i] = y1[i] + .5 * e2[i];
         }
 
         if( f(t + half_dt, y1, dy, params) ) {
@@ -117,9 +114,9 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
             return -1;
         }
 
-        for( j = 0; j < N; j++ ) {
-            e3[j] = dy[j] * dt;
-            y[j] = y0[j] + e3[j];
+        for( i = 0; i < N; i++ ) {
+            e3[i] = dy[i] * dt;
+            y[i] = y0[i] + e3[i];
         }
 
         if( f(t + dt, y1, dy, params) ) {
@@ -127,15 +124,13 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
             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;
+        for( i = 0; i < N; i++ ) {
+            e4 = dy[i] * dt;
+            y1[i] = y1[i] + (e1[i] + 2 * (e2[i] + e3[i]) + e4) / 6;
         }
 
         // Value can be printed/logged here
         logger(t, N, y1);
-
-        t = t0 + i * dt;
     }
 
     return 0;

modsim/ass3/logger.c

@@ -13,10 +13,10 @@ void logger(double t, int N, double *y) {
     if( !logfile )
         logfile = stdout;
 
-    fprintf(logfile, "%e", t);
+    fprintf(logfile, "%.12e", t);
 
     for( int i = 0; i < N; i++ )
-        fprintf(logfile, ",%e", y[i]);
+        fprintf(logfile, ",%.12e", y[i]);
 
     fprintf(logfile, "\n");
 }

modsim/ass3/main.c

@@ -5,6 +5,16 @@
 #include "models.h"
 #include "logger.h"
 
+#define DEFAULT_DT .001
+#define COMPARE_DT .5
+
+#define INTEGRATE(method, t0, t1, dt, N, FN) { \
+    if( method(t0, t1, dt, y0, y1, N, &f_ ## FN, (void *)args_ ## FN) ) \
+        puts("FAIL"); \
+}
+
+extern FILE *logfile;
+
 int main(int argc, char **argv) {
     if( argc < 2 ) {
         printf("Usage: %s ODE\n", argv[0]);
@@ -14,7 +24,7 @@ int main(int argc, char **argv) {
     double y0[3] = {.0, 5.0, .0},
            y1[3] = {.0, .0, .0};
 
-    const double args_osc[1]   = {.5},                         // {k}
+    const double args_osc[1]   = {1.0},                        // {k}
            args_forced_osc[3]  = {1.0, 0.1, 2.0},              // {k, D, w}
            args_lottka_volt[4] = {.5, 1.0, .1, .1},            // {a, b, c, d}
            args_modified_pp[5] = {.5, 1.0, .1, .1, 100.0},     // {a, b, c, d, y_m}
@@ -22,7 +32,19 @@ int main(int argc, char **argv) {
            *args_gilpin = NULL;
 
     switch( atoi(argv[1]) ) {
-        case 0: COMPARE(.0, 15.0, osc); break;
+        case 0:
+            logfile = fopen("report/osc_euler", "w");
+            puts("Integrating Euler method...");
+            INTEGRATE(Euler, .0, 40.0, COMPARE_DT, 2, osc);
+            logger_close();
+            logfile = fopen("report/osc_rk2", "w");
+            puts("Integrating Runge-Kutta 2 method...");
+            INTEGRATE(RungeKutta2, .0, 40.0, COMPARE_DT, 2, osc);
+            logger_close();
+            logfile = fopen("report/osc_rk4", "w");
+            puts("Integrating Runge-Kutta 4 method...");
+            INTEGRATE(RungeKutta4, .0, 40.0, COMPARE_DT, 2, osc);
+        break;
         case 1: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, forced_osc); break;
         case 2: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, lottka_volt); break;
         case 3: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, modified_pp); break;

modsim/ass3/models.h

 #ifndef __MODELS_H__
 #define __MODELS_H__
 
-#define DEFAULT_DT .001
-
-#define INTEGRATE(method, t0, t1, dt, N, FN) { \
-    if( method(t0, t1, dt, y0, y1, N, &f_ ## FN, (void *)args_ ## FN) ) \
-        puts("FAIL"); \
-}
-
-#define COMPARE(t0, t1, fn) { \
-    INTEGRATE(Euler, t0, t1, DEFAULT_DT, 2, fn); \
-    INTEGRATE(RungeKutta2, t0, t1, DEFAULT_DT, 2, fn); \
-    INTEGRATE(RungeKutta4, t0, t1, DEFAULT_DT, 2, fn); \
-}
-
 #define PARAM(n) ((const double *) params)[n]
 
 double *v_buf, *s_buf;

modsim/ass3/plot.py

@@ -3,6 +3,7 @@ from sys import argv, stdin
 from pylab import figure, clf, plot, show, savefig
 
 # Collect data
+lines = stdin.readlines()
 x = []
 y = [[] for i in range(lines[0].count(','))]
 for i in range(len(lines)):
@@ -18,7 +19,7 @@ for i in range(len(lines)):
 figure(1)
 clf()
 for i in range(len(y)):
-    plot(x, y[i], '-')
+    plot(x, y[i], 'x')
 if len(argv) == 2:
     savefig(argv[1])
 show()

modsim/ass3/report/report.tex

@@ -73,6 +73,25 @@ verassend, want RK4 voert de meeste operaties uit).
 
 \subsection{Harmonische oscillator}
 
+De vergelijking van de drie integratiemethoden kan worden uitgevoerd met:
+\begin{verbatim}
+$ sh ./compare_osc.sh
+\end{verbatim}
+
+In de volgende grafieken staan hiervan de resultaten (blauw
+
+\begin{figure}[H]
+    \includegraphics[scale=.5]{osc_euler.pdf}
+    \caption{Euler}
+\end{figure}
+\begin{figure}[H]
+    \includegraphics[scale=.5]{osc_rk2.pdf}
+    \caption{Runge-Kutta 2}
+\end{figure}
+\begin{figure}[H]
+    \includegraphics[scale=.5]{osc_rk4.pdf}
+    \caption{Runge-Kutta 4}
+\end{figure}
 
 \subsection{Damped, forced harmonic oscillator}
 

modsim/ass3/test.c

@@ -4,8 +4,6 @@
 #include "integration.h"
 #include "logger.h"
 
-#define DT .01
-
 extern FILE *logfile;
 
 int f_test_1(double t, double *y, double *dy, void *params) {
@@ -40,9 +38,20 @@ int f_test_4(double t, double *y, double *dy, void *params) {
     }
 }
 
+int f_cos(double t, double *y, double *dy, void *params) {
+    double result = cos(t);
+
+    if( isnan(result) ) {
+        return -1;
+    } else {
+        *dy = result;
+        return 0;
+    }
+}
+
 #define INTEGRATE(method, y_start, t0, t1, dy) { \
     *y0 = y_start; \
-    if( method(t0, t1, DT, y0, y1, 1, &dy, NULL) ) \
+    if( method(t0, t1, dt, y0, y1, 1, &dy, NULL) ) \
         puts("FAIL"); \
     else \
         printf("y1 for dy = " #dy " using %-12s %.10f\n", #method ":", *y1); \
@@ -53,9 +62,15 @@ int f_test_4(double t, double *y, double *dy, void *params) {
     INTEGRATE(RungeKutta4, y_start, t0, t1, dy); \
 }
 
-int main(void) {
-    double y0[1] = {.0}, y1[1] = {.0};
+int main(int argc, char **argv) {
+    double y0[1] = {.0}, y1[1] = {.0}, dt;
 
+    if( argc < 2 ) {
+        printf("Usage: %s DT\n", argv[0]);
+        exit(EXIT_FAILURE);
+    }
+
+    dt = atof(argv[1]);
     logfile = fopen("/dev/null", "w");
 
     COMPARE(.0, .0, 10.0, f_test_1);
@@ -67,6 +82,8 @@ int main(void) {
     COMPARE(1.0, -1.0, 1.0, f_test_4);
     puts("");
     COMPARE(-1.0, 1.0, 10.0, f_test_4);
+    puts("");
+    COMPARE(.0, .0, 5.0, f_cos);
 
-    return 0;
+    return EXIT_SUCCESS;
 }