I was reasoning about the behaviour of the methods I'm using for my simulation and I noticed that, considering $h_s$ as the timestep over which I have unstable solutions and $h_a$ as the timestep where I start to have the expected accuracy, there's an interval of timesteps $h_a < h < h_s$ for which I don't have the expected accuracy. So I was wondering what could be the reason of that and my conclusions are:

  1. my $h_s$ is overestimated and it is actually smaller
  2. there's something wrong in the method
  3. it's something normal, meaning that the first timesteps smaller than $h_s$ doesn't guarantee the expected accuracy, depending on the method and the simulated problem I'm considering.

I'm more convinced by the latter since the higher is $h$ the worst is the evaluated solution and its derivatives and this in turn could generate an error which doesn't behave as expected, but I'm not sure about it. Moreover there could be some other reason of course.


In the following plot the y axis is $\lVert u_{4h} - u_{2h}\rVert / \lVert u_h - u_{2h}\rVert$ over $2^p$, where $u$ is the numerical solution, so near the $y = 1$ I have the method (verlet) behaving correctly. enter image description here

The following is what I get for Runge-Kutta 4 for the same system instead (after the last point on the left I reach the rounding error)

enter image description here

  • $\begingroup$ You are using some Runge-Kutta method with fixed time step? Your ODE system converges to some stable equilibrium point? Your system has some stiffness, and $h_s$ times the spectrum is inside the stability region with points on the boundary? You are saying that for $h<h_a$ the error behaves strongly like $c(t)h^p$? Did you produce a stepsize-error loglog plot? $\endgroup$ Apr 2, 2021 at 20:08
  • $\begingroup$ What is the method that you're using? How are you measuring the accuracy? $\endgroup$
    – nicoguaro
    Apr 2, 2021 at 20:16
  • $\begingroup$ I'm using a Runge-Kutta(4), a velocity Verlet and an Adams-Bashforth(4), and what I described happen for all of them. The ode system is a 3-body problem with a star at its center and the other bodies orbiting around it. I share a picture of the error behaviour in log-log plot for Verlet as an example. $\endgroup$
    – Zebx
    Apr 2, 2021 at 20:42
  • 1
    $\begingroup$ For larger values of $h$ the leading error term does not dominate the whole error, so distortions in direction of the next power, $h^4$ for Verlet, $h^5$ for RK4, are to be expected. $\endgroup$ Apr 3, 2021 at 4:20
  • $\begingroup$ All-in-all, your point 3. applies. If $L$ stands for a Lipschitz constant of the system, like the modulus of the largest eigenvalue (all having negative imaginary part), then the methods of moderate order in double precision floating point numbers will give behavior according to the order of the method, that is, the leading error term clearly dominates, in a range of $Lh=0.1$ to $10^{-4}$. The boundary of the stability region (if any) can usually be found for $Lh$ between $2$ and $3.5$, so indeed there will be a gap. $\endgroup$ Oct 3, 2021 at 14:25

1 Answer 1


Stability does not necessarily imply accuracy.

I'll demonstrate this with a simple scalar ODE $y' = \lambda y$ (known as Dahlquist's test equation). This simple ODE is generally interpreted as a linearisation of a generic ODE $y'=f(y)$ around a given solution point, i.e. it reproduces, at least locally, the behaviour of the more complex ODE you might be considering. Here, $\lambda$ is the eigenvalue of the system.

Let's now apply forward-Euler to that equation: $$y_{n+1} = y_{n} + \lambda \Delta t y_{n}$$ Let us introduce $z=\lambda \Delta t$. We have: $$y_{n+1} = \underbrace{(1+z)}_{R(z)} y_n \tag{1}$$

The stability function is: $$R(z)=\frac{y_{n+1}}{y_n} = 1 + z$$

To ensure stability, you have to ensure that $|R(z)| \leq 1$, i.e. $z$ is (in the complex plane), in the disc centered on $-1 + 0i$ with radius $1$.

Now the true solution of the simple ODE is an exponential $y(t)=y(0) \exp(\lambda t)$. If we assume the initial condition $y_n$ for the current time step is equal to the true solution at that time, then the true solution at the next time step is:

$$y(t_{n+1}) = y_n \exp(\lambda\Delta t) = y_n \exp(z) \tag{2}$$

Compare Equations $(1)$ and $(2)$: we see that $R(z)$ is an approximation of $\exp(z)$. The actual Taylor expansion around 0 of the exponential is: $$\exp(z) = \sum\limits_{i=0}^{+\infty} \frac{z^n}{n!} = \underbrace{1 + z}_{R(z)} + \frac{z^2}{2} + ...$$ Therefore, we directly see that $R(z)$ from the forward Euler method has a local error of order 2, hence its order of convergence is 1 (for the global error).

Now we can study the error of the forward Euler solution (for one single time step)
relative to the true solution for different values of $z=\lambda \Delta t$: enter image description here

A quick explanation for this graph: The red curve is the limit of the stability domain, the contour field represent the log10 of the absolute relative error $|\frac{R(z)-\exp(z)}{\exp(z)}|$, and the black isolines represents the same error level but in percent. So if you want to have a precise solution, you would typically reduce the time steps so that $\lambda\Delta t$ is within the 1% precision zone (first black circle around 0) for all the eigenvalues $\lambda$ of your ODE system. The zone with black hatches correspond to thz one where your numerical approximation has more than 100% relative error to the true solution for a single time step.

Here is the same plot with RK4: enter image description here

We see that the zone of high-precision (typically look at the 1% isoline) is much larger, because $R(z)$ is then a much better approximation of the exponential. We also see that the solution may be precise even outside the stability domain. For example for eigenvalues with a positive real part (for example that may occur with chemical reactions), the true solution diverges, hence stability is not always a good indicator of the precision of your solution. Your numerical solution may be precise yet unstable, or it may be stable yet very imprecise...

So, coming back to what you observed (intermediate range of $\Delta t$ where the error is large), this corresponds to the zone where $\Delta t$ diminshes, therefore $\lambda \Delta t$ progressively nears 0, improving gradually the accuracy of the approximation $R(z)$ of the exponential.

I'll be happy to complement my answer if needed !

EDIT1: To see why the above test equation is indeed relevant, let s consider the general system of $N$ ODEs, written in vector the form $y' = f(y)$. We can linearize f around $y_0$: $$ y' \approx f(y_0) + \nabla f(y_0) (y-y_0) = \nabla f(y_0) y + \underbrace{(f(y_0) - \nabla f(y_0))}_{=A=constant} = \nabla f(y_0) + A$$ The constant term $A$ yields a linear term (wrt time) in the solution, which is trivially integrated exactly by any numerical integrator. The true challenge only comes from the non-constant part, which is not integrated exactly (except by exponential integrators, but I don't know much about these). Therefore the interesting equation is $\hat{y}'= \nabla f(y_0) \hat{y}$. Now, assuming you can diagonalize $\nabla f(y_0)$ (eigenvalues $\lambda_i$), you get an equivalent system of ODEs that is $\bar{y}_i' = \lambda_i \bar{y}_i, i=1..N$.

  • $\begingroup$ Thank you, your answer's pretty clear to me, a part from the linearization. I mean when I tried to do the same in the past for my ode system I found something wrong,so I've just focused on findng an empirical stability region for specific integration times and a set of timesteps. The problem is: I have a 6 equations system of ode, 2 for each dimension. As far as I understood I should evaluate the Jacobian of the function (a 6x6 matrix), then I find the eigenvalues and, if $Re(\lambda_k)<0$, I search for the set of $\Delta t$ such that (considering Runge-Kutta 4) $\lvert R(h\Delta t) \rvert<1$. $\endgroup$
    – Zebx
    Apr 3, 2021 at 11:23
  • $\begingroup$ So I should find the contours of $R(z)$. But some eigenvalue didn't have $Re(\lambda)<0$ and the condition $R(z)<1$ is never satisfied. I know it should be better to provide more details so I can make another question eventually. $\endgroup$
    – Zebx
    Apr 3, 2021 at 11:23
  • 1
    $\begingroup$ @Zebx I have added an edit to address the linearization. Regarding your second comment, finding $\Delta t$ such that all eigenvalues are within the stability domain only guarantees stability, but beware that, in a nonlinear problem, the eigenvalues change with the solution, so that choice of $\Delta t$ may become unstable after some time. It does not guarantee a given level of accuracy. $\endgroup$
    – Laurent90
    Apr 3, 2021 at 20:23
  • 1
    $\begingroup$ @Zebx for eigenvalues with a positive real part, this means that some components of the exact solution are growing (at least locally around the point of linearisation): requiring stability from your numerical solution means you are not describing the dynamics accurately... For eigenvalues with a positive real part, I guess this happens for a nonlinear system which may have eigenvalues moving from the left half complex plane to the right (for example eigenvalues associated with chemical reactions). If you have an eigenvalue that always has a positive real part, the true solution diverges. $\endgroup$
    – Laurent90
    Apr 3, 2021 at 20:29
  • $\begingroup$ The $y_0$ for the linearization of the problem can be any point of the solution (supposing it is sufficiently smooth to allow me to do that) or it can be just an equilibrium point of the problem? $\endgroup$
    – Zebx
    Apr 4, 2021 at 21:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.