∂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));
}

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));
}

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₀ and x₁₀₀ thus adding hair, but not in the inner loop:

x_−1/2 = x₀ − Δt v₀/2
perform the centered equations for n from 0 thru 99
x₁₀₀ = x_100−1/2 + Δt v₁₀₀/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));
}

Changing Δt

x₁₀₀ = x_100−1/2 + Δt v₁₀₀/2
multiply Δt by adjustment factor
x_100−1/2 = x₁₀₀ − Δt v₁₀₀/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));
}

Crypto Connection