Commit 02f831b4 authored by Taddeüs Kroes's avatar Taddeüs Kroes

ModSim ass3: Added log function.

parent 6d2f51b7
...@@ -18,6 +18,11 @@ All functions return 0 upon success, or -1 when a computational error ...@@ -18,6 +18,11 @@ All functions return 0 upon success, or -1 when a computational error
occurs (overflow, NaN or such) 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, int Euler(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; int i, j, intervals = (t1 - t0) / dt;
...@@ -27,6 +32,8 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -27,6 +32,8 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
for( j = 0; j < N; j++ ) for( j = 0; j < N; j++ )
y1[j] = y0[j]; y1[j] = y0[j];
LOG(t0, y0)
// Estimate for each interval // Estimate for each interval
for( i = 0; i < intervals; i++ ) { for( i = 0; i < intervals; i++ ) {
if( f(t, y1, dy, params) ) { if( f(t, y1, dy, params) ) {
...@@ -34,10 +41,11 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -34,10 +41,11 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
return -1; return -1;
} }
for( j = 0; j < N; j++ ) { for( j = 0; j < N; j++ )
y1[j] = y1[j] + dy[j] * dt; y1[j] = y1[j] + dy[j] * dt;
// Value can be printed/logged here // Value can be printed/logged here
} LOG(t, y1);
t = t0 + i * dt; t = t0 + i * dt;
} }
...@@ -54,6 +62,8 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -54,6 +62,8 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N,
for( j = 0; j < N; j++ ) for( j = 0; j < N; j++ )
y1[j] = y0[j]; y1[j] = y0[j];
LOG(t0, y0)
// Estimate for each interval // Estimate for each interval
for( i = 0; i < intervals; i++ ) { for( i = 0; i < intervals; i++ ) {
if( f(t, y1, dy1, params) || f(t + dt, y1, dy2, params) ) { if( f(t, y1, dy1, params) || f(t + dt, y1, dy2, params) ) {
...@@ -61,10 +71,11 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -61,10 +71,11 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N,
return -1; return -1;
} }
for( j = 0; j < N; j++ ) { for( j = 0; j < N; j++ )
y1[j] = y1[j] + .5 * (dy1[j] + dy2[j]) * dt; y1[j] = y1[j] + .5 * (dy1[j] + dy2[j]) * dt;
// Value can be printed/logged here // Value can be printed/logged here
} LOG(t, y1);
t = t0 + i * dt; t = t0 + i * dt;
} }
...@@ -81,6 +92,8 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -81,6 +92,8 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
for( j = 0; j < N; j++ ) for( j = 0; j < N; j++ )
y[j] = y1[j] = y0[j]; y[j] = y1[j] = y0[j];
LOG(t0, y0)
// Estimate for each interval // Estimate for each interval
for( i = 0; i < intervals; i++ ) { for( i = 0; i < intervals; i++ ) {
if( f(t, y1, dy, params) ) { if( f(t, y1, dy, params) ) {
...@@ -121,9 +134,11 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -121,9 +134,11 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
for( j = 0; j < N; j++ ) { for( j = 0; j < N; j++ ) {
e4 = dy[j] * dt; e4 = dy[j] * dt;
y1[j] = y1[j] + (e1[j] + 2 * (e2[j] + e3[j]) + e4) / 6; y1[j] = y1[j] + (e1[j] + 2 * (e2[j] + e3[j]) + e4) / 6;
// Value can be printed/logged here
} }
// Value can be printed/logged here
LOG(t, y1);
t = t0 + i * dt; t = t0 + i * dt;
} }
......
typedef int (*f_ptr)(double t, double *y, double *dy, void *params); typedef int (*f_ptr)(double t, double *y, double *dy, void *params);
typedef void (*log_function)(double t, double *y);
log_function log_func = NULL;
int Euler(double t0, double t1, double dt, double *y0, double *y1, int N, int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
f_ptr f, void *params); f_ptr f, void *params);
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment