From: W. Trevor King <wking@tremily.us>
Date: Wed, 13 Feb 2013 00:54:33 +0000 (-0500)
Subject: one_gaussian_bump.c: Add rk4 integrator (4th order Runge-Kutta)
X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=96ecb30d32b34ea90360a265a19c128b98a234b6;p=assignment-template.git

one_gaussian_bump.c: Add rk4 integrator (4th order Runge-Kutta)
---

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: