one_gaussian_bump.c: Add a midpoint integrator

diff --git a/one_gaussian_bump.c b/one_gaussian_bump.c
index f3a4646e2ca3942c970a2c008006a30f384890bc..a6739b7ed91b3bace7e14db86a4462b0c12214ae 100644 (file)
--- a/one_gaussian_bump.c
+++ b/one_gaussian_bump.c
@@ -19,6 +19,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #include <stdlib.h>       /* for exit(), EXIT_*, malloc(), free() */
 #include <stdio.h>        /* for stderr, *printf() */
+#include <string.h>       /* for memcpy() */
 #include <getopt.h>       /* for getopt_long() */
 #include <assert.h>       /* for assert() */
 #include <math.h>         /* for exp(), sqrt() */
@@ -36,10 +37,12 @@ static struct option long_options[] = {
        {"mass",           required_argument, 0, 'm' },
        {"bump-magnitude", required_argument, 0, 'k' },
        {"bump-width",     required_argument, 0, 'b' },
+       {"euler",          no_argument,       0, 'e' },
+       {"midpoint",       no_argument,       0, 'M' },
        {}
 };
 
-static char *short_options = "t:d:x:y:X:E:Y:m:k:b:";
+static char *short_options = "t:d:x:y:X:E:Y:m:k:b:eM";
 
 /* dynamic system parameters */
 typedef struct state_struct {
@@ -140,6 +143,35 @@ void euler_step(state_t *state, system_t *system, double dt)
        free(dx_dt);
 }
 
+void midpoint_step(state_t *state, system_t *system, double dt)
+{
+       double *x = NULL, *dx_dt = 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() */
+
+       memcpy((void *)x, (void *)state->x, x_size);
+       (*system->dx_dt_fn)(state, (void *)system, dx_dt);
+       /* take a half step */
+       for (i = 0; i < state->n; i++)
+               state->x[i] += dx_dt[i] * dt / 2;
+       state->t += dt/2;
+       /* evaluate slopes at the midpoint */
+       (*system->dx_dt_fn)(state, (void *)system, dx_dt);
+       /* take a full step using the midpoint slopes */
+       for (i = 0; i < state->n; i++)
+               state->x[i] = x[i] + dx_dt[i] * dt;
+       state->t += dt/2;       
+
+       /* release memory */
+       free(x);
+       free(dx_dt);
+}
+
 /* prints "t\tx[0]\tx[1]\t...\tx[n-1]\n"
  * If system->energy_fn is defined, also print "\tenergy" before the
  * newline.
@@ -160,7 +192,7 @@ void print_state(state_t *state, system_t *system)
 }
 
 void parse_args(int argc, char **argv, double *dt, double *t_max,
-                state_t *state, system_t *system)
+                state_t *state, system_t *system, step_fn_t *step_fn)
 {
        int ch;
        gaussian_state_t *gstate = (gaussian_state_t *)system->dx_dt_state;
@@ -200,6 +232,12 @@ void parse_args(int argc, char **argv, double *dt, double *t_max,
                case 'b':
                        gstate->b = safe_strtod(optarg, "bump-width");
                        break;
+               case 'e':
+                       *step_fn = &euler_step;
+                       break;
+               case 'M':
+                       *step_fn = &midpoint_step;
+                       break;
                case '?':
                        exit(EXIT_FAILURE);
                default:
@@ -228,7 +266,7 @@ int main(int argc, char **argv) {
        gaussian.b = 1;
        state.x[2] = -sqrt(2*0.5/system.m);  /* v = -sqrt(2E/m),  with E = 0.5 */
 
-       parse_args(argc, argv, &dt, &t_max, &state, &system);
+       parse_args(argc, argv, &dt, &t_max, &state, &system, &step_fn);
 
        print_state(&state, &system);
        while (state.t < t_max) {