I am studying the behavior of the Robertson chemical reaction,
$$\begin{array}{rl} \dot{x} &= -0.04 x + 10^4 y z\\ \dot{y} &= 0.04 x - 10^4 y z - 3 \times 10^7 y^2\\ \dot{z} &= 3 \times 10^7 y^2\end{array}$$
(source: Wikipedia on stiff equations).
I have noticed that numerical schemes can integrate this equation forwards in time but not backwards; that is, for a given initial condition $(x_0, y_0, z_0)$, integrating forward to $t = t_1$ and then reversing the integration back to $t = t_0$ yields a solution not equal to the initial condition.
I am wondering about the mathematical reason for this, and if there is any remedy. Is this a feature of all stiff equations? Have I missed something obvious?
Thanks in advance.
More detail: I chose the initial condition $(3, 1, 2)$ as an example, but I have observed this effect with other initial conditions as well. Let $t_0 = 0$ and $t_1 = 1$. For the first integration:
$t = 1$, $(x, y, z) = (3.00609, 4.00023 \times 10^{-6}, 2.99391)$
And for the reverse integration:
$t = 0$, $(x, y, z) = ( 3.00657, 4.00151 \times 10^{-6}, 2.99343 )$
I use the BDF integration method and have tried multiple implementations; all of the results are consistent with each other.
The problem in the system you have written down is that the $y,z$ dynamics, especially the $y$ dynamics, are stiff, meaning that they demonstrate multiple time scales. At the beginning $y$ is rapidly converted into $z$, then the process mostly settles down. In the backward in time dynamics the integrator never sees $y$ start to get significantly larger and so it does not refine the step sizes to become small enough to resolve the transient effect.
This problem is apparently a benchmark for integrators for stiff ODEs, so you should not expect an "easy" fix to correct the problem.