Added additional simulation models and cleaned Makefile.

modsim/ass3/Makefile

-CC=gcc
 CFLAGS=-Wall -Wextra -pedantic -std=c99 -D_GNU_SOURCE -g -O0
-LFLAGS=-lm
+LDFLAGS=-lm
 
-all: test harm_osc
+PROGS=test harm_osc
 
-test: test.o methods.o
-	$(CC) $(CFLAGS) $(LFLAGS) -o $@ $^
+.PHONY: all clean
 
-harm_osc: harm_osc.o methods.o
-	$(CC) $(CFLAGS) $(LFLAGS) -o $@ $^
+all: ${PROGS}
 
-clean:
-	rm *.o test harm_osc
+test: test.o logger.o methods.o
+harm_osc: harm_osc.o logger.o methods.o
+
+clean: ; rm -vf *.o ${PROGS}

modsim/ass3/harm_osc.c

@@ -2,40 +2,83 @@
 #include <stdio.h>
 #include <math.h>
 #include "methods.h"
+#include "logger.h"
 
 #define DT .001
 
-#define INTEGRATE(method, t0, t1) { \
-    if( method(t0, t1, DT, y0, y1, 2, &f_osc, &k) ) \
+#define INTEGRATE(method, t0, t1, FN) { \
+    if( method(t0, t1, DT, y0, y1, 2, &f_ ## FN, args_ ## FN) ) \
         puts("FAIL"); \
     else \
         printf("y1 using %-12s s:%.10f\n", #method ":", *y1); \
 }
-#define COMPARE(t0, t1) { \
-    INTEGRATE(Euler, t0, t1); \
-    INTEGRATE(RungeKutta2, t0, t1); \
-    INTEGRATE(RungeKutta4, t0, t1); \
+#define COMPARE(t0, t1, fn) { \
+    printf("Function %s\n", #fn); \
+    INTEGRATE(Euler, t0, t1, fn); \
+    INTEGRATE(RungeKutta2, t0, t1, fn); \
+    INTEGRATE(RungeKutta4, t0, t1, fn); \
 }
 
-extern log_function log_func;
+//void log(double t, double *y) {
+//    printf("%f,%f,%f\n", t, y[0], y[1]);
+//}
+//
 
-void log(double t, double *y) {
-    printf("%f,%f,%f\n", t, y[0], y[1]);
-}
+#define PARAM(n) ((double *) params)[n]
 
 int f_osc(double t, double *y, double *dy, void *params) {
+    // s' = v
     dy[0] = y[1];
-    dy[1] = -*(double *)(params) * y[0];
+    // v' = -k * s
+    dy[1] = - PARAM(0) * y[0];
+    return 0;
+}
 
+int f_forced_osc(double t, double *y, double *dy, void *params) {
+    // s' = v
+    dy[0] = y[1];
+    // v' = -k * s - D * v + sin(w * t)
+    dy[1] = - PARAM(0) * y[0] - PARAM(1) * y[1] + sin(PARAM(2) * t);
+    return 0;
+}
+
+int f_lottka_volt(double t, double *y, double *dy, void *params) {
+    // x' = -a * x + c * d * x * y
+    dy[0] = - PARAM(0) * y[0] + PARAM(2) * PARAM(3) * y[0] * y[1];
+    // y' = b * y - d * x * y
+    dy[1] = - PARAM(1) * y[1] - PARAM(3) * y[0] * y[1];
+    return 0;
+}
+
+int f_modified_pp(double t, double *y, double *dy, void *params) {
+    // x' = -a * x + c * d * x * y
+    dy[0] = - PARAM(0) * y[0] + PARAM(2) * PARAM(3) * y[0] * y[1];
+    // y' = b * ( 1 - y / y_m ) * y - d * x * y
+    //    = b * y - b * y * y / y_m - d * x * y
+    double by = PARAM(1) * y[1];
+    dy[1] = by - by * y[1] / PARAM(4) - PARAM(3) * y[0] * y[1];
     return 0;
 }
 
 int main(void) {
-    double k = .5, y0[2] = {.0, 5.0}, y1[2] = {.0, .0};
+    double y0[2] = {.0, 5.0},
+           y1[2] = {.0, .0};
+
+    double args_osc[1]         = {.5},                     // {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}
+    //printf("number of args: %d\n", sizeof(args_osc) / sizeof(double));
+    COMPARE(.0, 10.0, osc);
+    COMPARE(.0, 10.0, forced_osc);
 
-    log_func = &log;
+    // TODO: this model is not yet completed.
+    printf("Function: lottka_volt\n");
+    INTEGRATE(RungeKutta4, .0, 10.0, lottka_volt);
 
-    COMPARE(.0, 10.0)
+    // TODO: this model is not yet completed.
+    printf("Function: modified_pp\n");
+    INTEGRATE(RungeKutta4, .0, 10.0, modified_pp);
 
     return 0;
 }

modsim/ass3/logger.c

+#include "logger.h"
+
+FILE *logfile = NULL;
+
+void logger(double t, double *y) {
+    if( !logfile ) {
+        logfile = fopen("methods.log", "w");
+    }
+
+    fprintf(logfile, "%e,%e\n", t, y[0]);
+}
+
+void logger_close() {
+    if( logfile && fclose(logfile) ) {
+        printf("Failed to close logfile. Error no: %d", errno);
+    }
+}
+

modsim/ass3/logger.h

+#ifndef __LOG_H__
+#define __LOG_H__
+
+#include <stdio.h>
+#include <errno.h>
+
+//typedef void (*log_function)(double t, double *y);
+
+void logger(double t, double *y);
+void logger_close();
+
+#endif // __LOG_H__

modsim/ass3/methods.c

 #include <stdlib.h>
 #include <stdio.h>
 #include "methods.h"
+#include "logger.h"
 
 /*
 t0 – the start value t0 of t (input)
@@ -18,11 +19,6 @@ All functions return 0 upon success, or -1 when a computational error
 occurs (overflow, NaN or such)
 */
 
-#define LOG(t, y) { \
-    if( log_func != NULL ) \
-        log_func(t, y); \
-}
-
 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;
@@ -32,7 +28,7 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
     for( j = 0; j < N; j++ )
         y1[j] = y0[j];
 
-    LOG(t0, y0)
+    logger(t0, y0);
 
     // Estimate for each interval
     for( i = 0; i < intervals; i++ ) {
@@ -45,7 +41,7 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
             y1[j] = y1[j] + dy[j] * dt;
 
         // Value can be printed/logged here
-        LOG(t, y1);
+        logger(t, y1);
 
         t = t0 + i * dt;
     }
@@ -62,7 +58,7 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N,
     for( j = 0; j < N; j++ )
         y1[j] = y0[j];
 
-    LOG(t0, y0)
+    logger(t0, y0);
 
     // Estimate for each interval
     for( i = 0; i < intervals; i++ ) {
@@ -75,7 +71,7 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N,
             y1[j] = y1[j] + .5 * (dy1[j] + dy2[j]) * dt;
 
         // Value can be printed/logged here
-        LOG(t, y1);
+        logger(t, y1);
 
         t = t0 + i * dt;
     }
@@ -92,7 +88,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
     for( j = 0; j < N; j++ )
         y[j] = y1[j] = y0[j];
 
-    LOG(t0, y0)
+    logger(t0, y0);
 
     // Estimate for each interval
     for( i = 0; i < intervals; i++ ) {
@@ -137,7 +133,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
         }
 
         // Value can be printed/logged here
-        LOG(t, y1);
+        logger(t, y1);
 
         t = t0 + i * dt;
     }

modsim/ass3/methods.h

-typedef int (*f_ptr)(double t, double *y, double *dy, void *params);
-typedef void (*log_function)(double t, double *y);
+#ifndef __METHODS_H__
+#define __METHODS_H__
 
-log_function log_func = NULL;
+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);
@@ -10,3 +10,5 @@ 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);
+
+#endif // __METHODS_H__

modsim/ass3/plot.py

+#!/usr/bin/env python
 from sys import argv, exit
 from pylab import figure, clf, plot, show, savefig
 

statred/multivariate.py

-mean = [0,0]
-cov = [[1,0],[0,100]] # diagonal covariance, points lie on x or y-axis
-
-import matplotlib.pyplot as plt
-import numpy as np
-x,y = np.random.multivariate_normal(mean, cov, 5000).T
-print 'x:', x
-print 'y:', y
-plt.plot(x,y,'x'); plt.axis('equal'); plt.show()
-