From 96ecb30d32b34ea90360a265a19c128b98a234b6 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" <wking@tremily.us> Date: Tue, 12 Feb 2013 19:54:33 -0500 Subject: [PATCH] one_gaussian_bump.c: Add rk4 integrator (4th order Runge-Kutta) --- one_gaussian_bump.c | 63 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 62 insertions(+), 1 deletion(-) diff --git a/one_gaussian_bump.c b/one_gaussian_bump.c index a6739b7..1124022 100644 --- a/one_gaussian_bump.c +++ b/one_gaussian_bump.c @@ -39,10 +39,11 @@ static struct option long_options[] = { {"bump-width", required_argument, 0, 'b' }, {"euler", no_argument, 0, 'e' }, {"midpoint", no_argument, 0, 'M' }, + {"rk4", no_argument, 0, 'r' }, {} }; -static char *short_options = "t:d:x:y:X:E:Y:m:k:b:eM"; +static char *short_options = "t:d:x:y:X:E:Y:m:k:b:eMr"; /* dynamic system parameters */ typedef struct state_struct { @@ -172,6 +173,63 @@ void midpoint_step(state_t *state, system_t *system, double dt) free(dx_dt); } +void rk4_step(state_t *state, system_t *system, double dt) +{ + double *x = NULL, *dx_dt = NULL; + double *k1 = NULL, *k2 = NULL, *k3 = NULL, *k4 = NULL; + size_t i, x_size = sizeof(double) * state->n; + + /* allocate memory */ + x = malloc(x_size); + assert(x); /* check for successful malloc() */ + dx_dt = malloc(x_size); + assert(dx_dt); /* check for successful malloc() */ + k1 = malloc(x_size); + assert(k1); /* check for successful malloc() */ + k2 = malloc(x_size); + assert(k2); /* check for successful malloc() */ + k3 = malloc(x_size); + assert(k3); /* check for successful malloc() */ + k4 = malloc(x_size); + assert(k4); /* check for successful malloc() */ + + memcpy((void *)x, (void *)state->x, x_size); + /* k1 */ + (*system->dx_dt_fn)(state, (void *)system, dx_dt); + for (i = 0; i < state->n; i++) { + k1[i] = dx_dt[i] * dt; + state->x[i] = x[i] + 0.5*k1[i]; + } + state->t += dt/2; + /* k2 */ + (*system->dx_dt_fn)(state, (void *)system, dx_dt); + for (i = 0; i < state->n; i++) { + k2[i] = dx_dt[i] * dt; + state->x[i] = x[i] + 0.5*k2[i]; + } + /* k3 */ + (*system->dx_dt_fn)(state, (void *)system, dx_dt); + for (i = 0; i < state->n; i++) { + k3[i] = dx_dt[i] * dt; + state->x[i] = x[i] + k3[i]; + } + state->t += dt/2; + /* k4 */ + (*system->dx_dt_fn)(state, (void *)system, dx_dt); + for (i = 0; i < state->n; i++) { + k4[i] = dx_dt[i] * dt; + state->x[i] = x[i] + (k1[i] + 2*k2[i] + 2*k3[i] + k4[i])/6.0; + } + + /* release memory */ + free(x); + free(dx_dt); + free(k1); + free(k2); + free(k3); + free(k4); +} + /* prints "t\tx[0]\tx[1]\t...\tx[n-1]\n" * If system->energy_fn is defined, also print "\tenergy" before the * newline. @@ -238,6 +296,9 @@ void parse_args(int argc, char **argv, double *dt, double *t_max, case 'M': *step_fn = &midpoint_step; break; + case 'r': + *step_fn = &rk4_step; + break; case '?': exit(EXIT_FAILURE); default: -- 2.26.2