∂x/∂t = v

∂v/∂t = −x

A naïve first order numerical implementation of these equations would be:

x_{n+1} = x_{n} + Δt v_{n}

v_{n+1} = v_{n} − Δt x_{n}

This is symmetrical between x and v. The matching C code:

#include <stdio.h> #include <math.h> int main(){ double x = 0, v = 1, dt = 0.01; int j; for(j = 0; j<100; ++j) {double tmp = x; x = x + dt*v; v = v - dt*tmp;} printf("%13.9f %13.9f, %13.9f %13.9f\n", x, v, x-sin(1), v-cos(1)); }yields: 0.845670565 0.543038634, 0.004199580 0.002736328

suggesting an error of about Δt/2 which is consistent with a first order scheme.

The lazy or prescient programmer might be tempted to write instead

#include <stdio.h> #include <math.h> int main(){ double x = 0, v = 1, dt = 0.01; int j; for(j = 0; j<100; ++j) { x = x + dt*v; v = v - dt*x;} printf("%13.9f %13.9f, %13.9f %13.9f\n", x, v, x-sin(1), v-cos(1)); }which yields: 0.841483755 0.536091381, 0.000012770 −0.004210925

and this suggests a 2nd order scheme yielding an error in x of Δt

x_{n+1} = x_{n} + Δt v_{n}

v_{n+1} = v_{n} − Δt x_{n+1}

This has probably been discovered more frequently than designed. We must understand what is going on here if we are to hold to this advantage in more complex cases. And what about that Δt/2 error in v?

Here is one way to make sense of this luck, and thus to pretend that we invented it. We modify the difference equations above to:

x_{n+1/2} = x_{n} + Δt v_{n}/2

v_{n+1} = v_{n} − Δt x_{n+1/2}

x_{n+1} = x_{n+1/2} + Δt v_{n+1}/2

In words, we evaluate the influence of v on the rate of change of x at a time halfway thru the time interval between n and n+1. The trick works on differential equations of the form.

∂x/∂t = f(v)

∂v/∂t = g(x)

The diffusion equation, or Coriolis force do not fit this form.
The last of these equations can be combined with the first equation of the next time step and thus entirely removing the x_{n} values from the computation.
The new **centered** difference equations are

x_{n+1/2} = x_{n−1/2} + Δt v_{n}

v_{n+1} = v_{n} − Δt x_{n+1/2}

We are then back to the simple discovered equations, having renamed the x values, but there is now a theory.

The assimilation trick is not proper for the beginning, where we established initial conditions,
nor at the ending, where we report the results.
To be fussy we reintroduce x_{0} and x_{100} thus adding hair, but not in the
inner loop:

x_{−1/2} = x_{0} − Δt v_{0}/2

perform the centered equations for n from 0 thru 99

x_{100} = x_{100−1/2} + Δt v_{100}/2

The matching code is

#include <stdio.h> #include <math.h> int main(){ double x = 0, v = 1, dt = 0.01; int j; x -= dt*v/2; for(j = 0; j<100; ++j) { x = x + dt*v; v = v - dt*x;} x += dt*v/2; printf("%13.9f %13.9f, %13.9f %13.9f\n", x, v, x-sin(1), v-cos(1)); }which yields: 0.841462718 0.540298800, −0.000008267 −0.000003506

suggesting an error of Δt

x_{100} = x_{100−1/2} + Δt v_{100}/2

multiply Δt by adjustment factor

x_{100−1/2} = x_{100} − Δt v_{100}/2

At this point it is best to add t += Δt to the centered equations to avoid confusion. Thus t, not iteration count, tracks progress.

To integrate to two radians, changing Δt midway, we now have the code:

#include <stdio.h> #include <math.h> int main(){ double x = 0, v = 1, t=0, dt = 0.01; int j; x -= dt*v/2; for(j = 0; j<100; ++j) { x += dt*v; v -= dt*x; t += dt;} {x += dt*v/2; dt *= 2; x -= dt*v/2;} // adjust dt. for(j = 0; j<50; ++j) { x += dt*v; v -= dt*x; t += dt;} x += dt*v/2; printf("%13.9f %13.9f %13.9f, %13.9f %13.9f\n", t, x, v, x-sin(2), v-cos(2)); }This yields 2.000000000 0.909260340 −0.416192336, −0.000037086 −0.000045500

which is consonant with the longer interval of integration and the double Δt.