Harm osc will write log to stdout.

parent 0db05df5
...@@ -7,23 +7,18 @@ ...@@ -7,23 +7,18 @@
#define DT .001 #define DT .001
#define INTEGRATE(method, t0, t1, FN) { \ #define INTEGRATE(method, t0, t1, FN) { \
if( method(t0, t1, DT, y0, y1, 2, &f_ ## FN, args_ ## FN) ) \ method(t0, t1, DT, y0, y1, 2, &f_ ## FN, args_ ## FN); \
puts("FAIL"); \ /*if( method(t0, t1, DT, y0, y1, 2, &f_ ## FN, args_ ## FN) ) puts("FAIL");*/ \
else \ /* else printf("y1 using %-12s s:%.10f\n", #method ":", *y1); */\
printf("y1 using %-12s s:%.10f\n", #method ":", *y1); \
} }
#define COMPARE(t0, t1, fn) { \ #define COMPARE(t0, t1, fn) { \
printf("Function %s\n", #fn); \ /*printf("Function %s\n", #fn);*/ \
INTEGRATE(Euler, t0, t1, fn); \ INTEGRATE(Euler, t0, t1, fn); \
INTEGRATE(RungeKutta2, t0, t1, fn); \ INTEGRATE(RungeKutta2, t0, t1, fn); \
INTEGRATE(RungeKutta4, t0, t1, fn); \ INTEGRATE(RungeKutta4, t0, t1, fn); \
} }
//void log(double t, double *y) {
// printf("%f,%f,%f\n", t, y[0], y[1]);
//}
//
#define PARAM(n) ((double *) params)[n] #define PARAM(n) ((double *) params)[n]
int f_osc(double t, double *y, double *dy, void *params) { int f_osc(double t, double *y, double *dy, void *params) {
...@@ -73,12 +68,14 @@ int main(void) { ...@@ -73,12 +68,14 @@ int main(void) {
COMPARE(.0, 10.0, forced_osc); COMPARE(.0, 10.0, forced_osc);
// TODO: this model is not yet completed. // TODO: this model is not yet completed.
printf("Function: lottka_volt\n"); //printf("Function: lottka_volt\n");
INTEGRATE(RungeKutta4, .0, 10.0, lottka_volt); INTEGRATE(RungeKutta4, .0, 10.0, lottka_volt);
// TODO: this model is not yet completed. // TODO: this model is not yet completed.
printf("Function: modified_pp\n"); //printf("Function: modified_pp\n");
INTEGRATE(RungeKutta4, .0, 10.0, modified_pp); INTEGRATE(RungeKutta4, .0, 10.0, modified_pp);
logger_close();
return 0; return 0;
} }
...@@ -2,12 +2,24 @@ ...@@ -2,12 +2,24 @@
FILE *logfile = NULL; FILE *logfile = NULL;
void logger(double t, double *y) { void logger_open(const char* filename) {
if( !logfile ) { if( logfile )
logfile = fopen("methods.log", "w"); logger_close();
} logfile = fopen(filename, "w");
}
fprintf(logfile, "%e,%e\n", t, y[0]); void logger(double t, int N, double *y) {
if( !logfile )
logfile = stdout;
fprintf(logfile, "%e", t);
if( N > 0 ) {
for( int i = 0; i < N - 1; i++ )
fprintf(logfile, ",%e", y[0]);
fprintf(logfile, "%e\n", y[N-1]);
}
else
fprintf(logfile, "\n");
} }
void logger_close() { void logger_close() {
......
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
//typedef void (*log_function)(double t, double *y); //typedef void (*log_function)(double t, double *y);
void logger(double t, double *y); void logger(double t, int N, double *y);
void logger_close(); void logger_close();
#endif // __LOG_H__ #endif // __LOG_H__
...@@ -28,7 +28,7 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -28,7 +28,7 @@ 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];
logger(t0, y0); logger(t0, N, y0);
// Estimate for each interval // Estimate for each interval
for( i = 0; i < intervals; i++ ) { for( i = 0; i < intervals; i++ ) {
...@@ -41,7 +41,7 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -41,7 +41,7 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
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
logger(t, y1); logger(t, N, y1);
t = t0 + i * dt; t = t0 + i * dt;
} }
...@@ -58,7 +58,7 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -58,7 +58,7 @@ 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];
logger(t0, y0); logger(t0, N, y0);
// Estimate for each interval // Estimate for each interval
for( i = 0; i < intervals; i++ ) { for( i = 0; i < intervals; i++ ) {
...@@ -71,7 +71,7 @@ int RungeKutta2(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -71,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; y1[j] = y1[j] + .5 * (dy1[j] + dy2[j]) * dt;
// Value can be printed/logged here // Value can be printed/logged here
logger(t, y1); logger(t, N, y1);
t = t0 + i * dt; t = t0 + i * dt;
} }
...@@ -88,7 +88,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -88,7 +88,7 @@ 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];
logger(t0, y0); logger(t0, N, y0);
// Estimate for each interval // Estimate for each interval
for( i = 0; i < intervals; i++ ) { for( i = 0; i < intervals; i++ ) {
...@@ -133,7 +133,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -133,7 +133,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
} }
// Value can be printed/logged here // Value can be printed/logged here
logger(t, y1); logger(t, N, y1);
t = t0 + i * dt; t = t0 + i * dt;
} }
......
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