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

diff --git a/one_gaussian_bump.c b/one_gaussian_bump.c
index a6739b7ed91b3bace7e14db86a4462b0c12214ae..1124022144d069ec733bcdb6a8286e7c97b9ccdf 100644 (file)
--- 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: