Our equations are among those of physics that achieve 2nd order accuracy using a first order scheme. Variable Δt can forfeit this advantage if not done carefully. Perhaps the familiar coupled equations for the harmonic oscillator are the simplest to illustrate the opportunities and pitfalls here.

∂x/∂t = v
∂v/∂t = −x

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

xn+1 = xn + Δt vn
vn+1 = vn − Δt xn

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 Δt2/8. The corresponding difference equations are

xn+1 = xn + Δt vn
vn+1 = vn − Δt xn+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:

xn+1/2 = xn + Δt vn/2
vn+1 = vn − Δt xn+1/2
xn+1 = xn+1/2 + Δt vn+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 xn values from the computation. The new centered difference equations are

xn+1/2 = xn−1/2 + Δt vn
vn+1 = vn − Δt xn+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 x0 and x100 thus adding hair, but not in the inner loop:

x−1/2 = x0 − Δt v0/2
perform the centered equations for n from 0 thru 99
x100 = x100−1/2 + Δt v100/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 Δt2/12 or better. Serendipity!

### Changing Δt

Back to the original subject; if we merely change Δt between two performances of the centered equations, we introduce an error proportional to Δt thus forfeiting our Δt2 advantage. If, instead, we time align our x and v values at this point this error is avoided. In the above example if we wanted to double Δt before we proceed to the two radian point, we would adjust Δt while the variable x held the value of x100. The total adjustment between iterations of the centered equations is thus:

x100 = x100−1/2 + Δt v100/2
multiply Δt by adjustment factor
x100−1/2 = x100 − Δt v100/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.

### Crypto Connection

This scheme fits in the pattern of Feistel Structures which show up in crypto systems. Indeed if we implemented the dynamic variables above with data types for which additions was associative we could run the calculation backwards exactly. Truncation and rounding would be exactly reproduced in the opposite direction. The moral here is that we should not consider the ability to run backwards as evidence of accuracy.