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

Question

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.

Answer 1

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.

Answer 2

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.