Skip to content
Snippets Groups Projects
Commit d33225f4 authored by Taddeüs Kroes's avatar Taddeüs Kroes
Browse files

Implemented Gilpin's model.

parent b6fa77bc
No related branches found
No related tags found
No related merge requests found
...@@ -8,6 +8,6 @@ PROGS=test main ...@@ -8,6 +8,6 @@ PROGS=test main
all: ${PROGS} all: ${PROGS}
test: test.o logger.o integration.o test: test.o logger.o integration.o
harm_osc: main.o models.o logger.o integration.o main: main.o models.o logger.o integration.o
clean: ; rm -vf *.o ${PROGS} clean: ; rm -vf *.o ${PROGS}
#include <stdlib.h> #include <stdlib.h>
#include <stdio.h> #include <stdio.h>
#include "methods.h" #include "integration.h"
#include "logger.h" #include "logger.h"
/* /*
......
File moved
#include <stdlib.h> #include <stdlib.h>
#include <stdio.h> #include <stdio.h>
#include <math.h> #include <math.h>
#include "methods.h" #include "integration.h"
#include "models.h" #include "models.h"
#include "logger.h" #include "logger.h"
...@@ -11,25 +11,25 @@ int main(int argc, char **argv) { ...@@ -11,25 +11,25 @@ int main(int argc, char **argv) {
return EXIT_FAILURE; return EXIT_FAILURE;
} }
double y0[2] = {.0, 5.0}, double y0[3] = {.0, 5.0, .0},
y1[2] = {.0, .0}; y1[3] = {.0, .0, .0};
const double args_osc[1] = {.5}, // {k} const double args_osc[1] = {.5}, // {k}
args_forced_osc[3] = {1.0, 0.1, 2.0}, // {k, D, w} 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_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} args_modified_pp[5] = {.5, 1.0, .1, .1, 100.0}, // {a, b, c, d, y_m}
args_pid_control[5] = {-2.0, -1.0, .1, 10.0, .001}; // {P, D, a, M, d} args_pid_control[5] = {-2.0, -1.0, .1, 10.0, .001}, // {P, D, a, M, d}
//printf("number of args: %d\n", sizeof(args_osc) / sizeof(double)); *args_gilpin = NULL;
switch( atoi(argv[1]) ) { switch( atoi(argv[1]) ) {
case 0: COMPARE(.0, 10.0, osc); break; case 0: COMPARE(.0, 15.0, osc); break;
case 1: INTEGRATE(RungeKutta4, .0, 20.0, forced_osc); break; case 1: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, forced_osc); break;
case 2: INTEGRATE(RungeKutta4, .0, 20.0, lottka_volt); break; case 2: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, lottka_volt); break;
case 3: INTEGRATE(RungeKutta4, .0, 20.0, modified_pp); break; case 3: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, modified_pp); break;
case 4: case 4:
y0[0] = .0; // = v(0) y0[0] = .0; // = v(0)
y0[1] = 3.0; // = s(0) y0[1] = 3.0; // = s(0)
buf_size = ceil(args_pid_control[4] / DT); buf_size = ceil(args_pid_control[4] / DEFAULT_DT);
v_buf = malloc(buf_size * sizeof(double)); v_buf = malloc(buf_size * sizeof(double));
s_buf = malloc(buf_size * sizeof(double)); s_buf = malloc(buf_size * sizeof(double));
...@@ -38,10 +38,14 @@ int main(int argc, char **argv) { ...@@ -38,10 +38,14 @@ int main(int argc, char **argv) {
s_buf[i] = y0[1]; s_buf[i] = y0[1];
} }
INTEGRATE(RungeKutta4, .0, 100.0, pid_control); INTEGRATE(RungeKutta4, .0, 100.0, DEFAULT_DT, 2, pid_control);
free(v_buf); free(v_buf);
free(s_buf); free(s_buf);
break; break;
case 5:
y0[0] = y0[1] = y0[2] = 10.0;
INTEGRATE(RungeKutta4, .0, 1000.0, .1, 3, gilpin);
break;
default: puts("ODE does not exist!"); return EXIT_FAILURE; default: puts("ODE does not exist!"); return EXIT_FAILURE;
} }
......
#include <math.h> #include <math.h>
#include "models.h" #include "models.h"
int pos = 0;
int f_osc(double t, double *y, double *dy, void *params) { int f_osc(double t, double *y, double *dy, void *params) {
// s' = v // s' = v
dy[0] = y[1]; dy[0] = y[1];
// v' = -k * s // v' = -k * s
dy[1] = - PARAM(0) * y[0]; dy[1] = - PARAM(0) * y[0];
return 0; return 0;
} }
...@@ -16,6 +15,7 @@ int f_forced_osc(double t, double *y, double *dy, void *params) { ...@@ -16,6 +15,7 @@ int f_forced_osc(double t, double *y, double *dy, void *params) {
dy[0] = y[1]; dy[0] = y[1];
// v' = -k * s - D * v + sin(w * t) // v' = -k * s - D * v + sin(w * t)
dy[1] = - PARAM(0) * y[0] - PARAM(1) * y[1] + sin(PARAM(2) * t); dy[1] = - PARAM(0) * y[0] - PARAM(1) * y[1] + sin(PARAM(2) * t);
return 0; return 0;
} }
...@@ -24,6 +24,7 @@ int f_lottka_volt(double t, double *y, double *dy, void *params) { ...@@ -24,6 +24,7 @@ int f_lottka_volt(double t, double *y, double *dy, void *params) {
dy[0] = - PARAM(0) * y[0] + PARAM(2) * PARAM(3) * y[0] * y[1]; dy[0] = - PARAM(0) * y[0] + PARAM(2) * PARAM(3) * y[0] * y[1];
// y' = b * y - d * x * y // y' = b * y - d * x * y
dy[1] = PARAM(1) * y[1] - PARAM(3) * y[0] * y[1]; dy[1] = PARAM(1) * y[1] - PARAM(3) * y[0] * y[1];
return 0; return 0;
} }
...@@ -34,9 +35,12 @@ int f_modified_pp(double t, double *y, double *dy, void *params) { ...@@ -34,9 +35,12 @@ int f_modified_pp(double t, double *y, double *dy, void *params) {
// = b * y - b * y * y / y_m - d * x * y // = b * y - b * y * y / y_m - d * x * y
double by = PARAM(1) * y[1]; double by = PARAM(1) * y[1];
dy[1] = by - by * y[1] / PARAM(4) - PARAM(3) * y[0] * y[1]; dy[1] = by - by * y[1] / PARAM(4) - PARAM(3) * y[0] * y[1];
return 0; return 0;
} }
int pos = 0;
int f_pid_control(double t, double *y, double *dy, void *params) { int f_pid_control(double t, double *y, double *dy, void *params) {
int buf_begin = (pos + 1) % buf_size; int buf_begin = (pos + 1) % buf_size;
...@@ -55,3 +59,14 @@ int f_pid_control(double t, double *y, double *dy, void *params) { ...@@ -55,3 +59,14 @@ int f_pid_control(double t, double *y, double *dy, void *params) {
return 0; return 0;
} }
int f_gilpin(double t, double *y, double *dy, void *params) {
// x' = x(1 - 0.001x - 0.001y - 0.01z)
dy[0] = y[0] * (1.0 - .001 * y[0] - .001 * y[1] - .01 * y[2]);
// y' = y(1 - 0.001y - 0.0015x - 0.001z)
dy[1] = y[1] * (1.0 - .001 * y[1] - .0015 * y[0] - .001 * y[2]);
// z' = z(0.005x - 0.0005y - 1)
dy[2] = y[2] * (.005 * y[0] + .0005 * y[1] - 1.0);
return 0;
}
#ifndef __MODELS_H__ #ifndef __MODELS_H__
#define __MODELS_H__ #define __MODELS_H__
#define DT .001 #define DEFAULT_DT .001
#define INTEGRATE(method, t0, t1, FN) { \ #define INTEGRATE(method, t0, t1, dt, N, FN) { \
method(t0, t1, DT, y0, y1, 2, &f_ ## FN, (void *) args_ ## FN); \ if( method(t0, t1, dt, y0, y1, N, &f_ ## FN, (void *)args_ ## FN) ) \
/*if( method(t0, t1, DT, y0, y1, 2, &f_ ## FN, args_ ## FN) ) puts("FAIL");*/ \ puts("FAIL"); \
/* else printf("y1 using %-12s s:%.10f\n", #method ":", *y1); */\
} }
#define COMPARE(t0, t1, fn) { \ #define COMPARE(t0, t1, fn) { \
/*printf("Function %s\n", #fn);*/ \ INTEGRATE(Euler, t0, t1, DEFAULT_DT, 2, fn); \
INTEGRATE(Euler, t0, t1, fn); \ INTEGRATE(RungeKutta2, t0, t1, DEFAULT_DT, 2, fn); \
INTEGRATE(RungeKutta2, t0, t1, fn); \ INTEGRATE(RungeKutta4, t0, t1, DEFAULT_DT, 2, fn); \
INTEGRATE(RungeKutta4, t0, t1, fn); \
} }
#define PARAM(n) ((const double *) params)[n] #define PARAM(n) ((const double *) params)[n]
...@@ -26,5 +24,6 @@ int f_forced_osc(double t, double *y, double *dy, void *params); ...@@ -26,5 +24,6 @@ int f_forced_osc(double t, double *y, double *dy, void *params);
int f_lottka_volt(double t, double *y, double *dy, void *params); int f_lottka_volt(double t, double *y, double *dy, void *params);
int f_modified_pp(double t, double *y, double *dy, void *params); int f_modified_pp(double t, double *y, double *dy, void *params);
int f_pid_control(double t, double *y, double *dy, void *params); int f_pid_control(double t, double *y, double *dy, void *params);
int f_gilpin(double t, double *y, double *dy, void *params);
#endif // __MODELS_H__ #endif // __MODELS_H__
...@@ -2,10 +2,10 @@ ...@@ -2,10 +2,10 @@
from sys import argv, stdin from sys import argv, stdin
from pylab import figure, clf, plot, show, savefig from pylab import figure, clf, plot, show, savefig
try: #try:
f = open(argv[1], 'r') # f = open(argv[1], 'r')
except IndexError: #except IndexError:
f = stdin f = stdin
lines = f.readlines() lines = f.readlines()
f.close() f.close()
x = [] x = []
...@@ -22,6 +22,6 @@ figure(1) ...@@ -22,6 +22,6 @@ figure(1)
clf() clf()
for i in range(len(y)): for i in range(len(y)):
plot(x, y[i], '-') plot(x, y[i], '-')
if len(argv) == 3: if len(argv) == 2:
savefig(argv[2]) savefig(argv[1])
show() show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment