# Order of numerical solver when calculating difference between forwards and backwards solution

I'm working in applied oceanography, where people are sometimes interested in calculating backwards trajectories'' of things floating on the ocean, i.e., going backwards in time to figure out where something came from. In trying out some toy examples, I came across the following point, which I thought was curious.

I'm solving an ODE

$$\dot{x} = f(x),$$

where $$x$$ is a 2-component vector, and $$f(x)$$ is a 2D vector field (constant in time). I'm using the 4th-order Runge-Kutta method, and I ran two tests:

• I calculated a trajectory (a solution) forwards in time, starting at $$x_0$$ and integrating over $$t = [0, T]$$, then I reversed the vector field (multiplying both components with $$-1$$), and calculated a trajectory backwards in time from $$x(T)$$ over the interval $$[T, 0]$$, ending up approximately back at $$x_0$$. Then I repeated the whole procedure for several timesteps, and plotted the error (final position - $$x_0$$), as a function of timestep.
• I calculated a forwards trajectory with a very short timestep, over the interval $$[0, T]$$, and called this $$x(T)$$ my reference solution. Then I repeated the exercise for several longer timesteps, and calculated the error relative to the reference solution, and plotted this as well as a function of timestep.

In the first case, I found that the error scaled as $$\Delta t^5$$, while in the second case it scaled as $$\Delta t^4$$, which is of course what I would expect from 4th-order Runge-Kutta. Results shown in the plot below.

My question is: What is the reason that I gain an extra order in the first case. I assume it must be something to do with the fact that going backwards along the same trajectory somehow cancels out the error in the pairwise opposing steps, but a more rigorous argument would be nice.

Update: Based on the suggestion from ConvexHull, I re-ran the analysis with a 3rd-order Runge-Kutta method (Kutta's 3rd-order method, https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Kutta's_third-order_method). Results in the plot below.

The paths of the forward and backward integration have a global error of $$O(h^4)$$, their distance is thus at most that large. In consequence, contributions of that difference to the coefficients of the error expansion do not occur in the first 3 coefficients.

The error expansion can thus be written as $$e(t,h)=c_5(t)h^5+c_6(t)h^6+...$$ for the local single step error for the step from $$t$$ to $$t+h$$.

For the same interval but backwards the error is $$e(t+h,-h)=-c_5(t+h)h^5+c_6(t+h)h^6+...$$ Forgetting for this exploration that the error accumulation is not fully additive, the combined error of forward and backward step on the interval $$[t,t+h]$$ is $$e(t,h)+e(t+h,-h)=(c_5(t)-c_5(t+h))h^5+(c_6(t)+c_6(t+h))h^6+...$$ As the coefficients are also smooth functions, the difference in the first term is proportional to $$h$$, making the first term $$O(h^6)$$. This means that the fully accumulated difference between the forward and backward integrations is $$O(h^5)$$ as observed.

If the parity of the degrees of the first error terms is reversed as in a third order method, then the first term does not cancel out, the order is preserved in the combined error and thus the distance is also of the same order.

What you observe should be an even-odd related problem. Recall that even- and odd functions are defined as

\begin{align} \text{even function:}\quad f(t)-f(-t)=0,\\ \text{odd function:}\quad f(t)+f(-t)=0.\\ \end{align}

Since you are using an even order Runge-Kutta method the error (a̲l̲l̲ e̲v̲e̲n̲ r̲e̲l̲a̲t̲e̲d̲ c̲o̲n̲t̲r̲i̲b̲u̲t̲i̲o̲n̲s̲) should cancel out when you reverse the time-stepping. In the same way you should observe the same order of convergence for both approaches if you are using a Runge-Kutta method of odd order.