I'm developing a solver for solving linear hyperbolic equations of first order with respect to time and spatial derivatives. The formal order of accuracy of the solver must be 5 because I use 5th-order accurate Runge--Kutta method and 5th-order finite-difference WENO for Hamilton--Jacobi equations.

To verify my solver, I do a convergence study. I have to cases: in one case I have an exponentially growing solution, and in another case I have an exponentially decaying solution. I do not know an exact solution, so I employ the solution on the finest grid as an exact one to measure the discretization error.

After measuring the error in $L_2$-norm, I obtain the following results for the case of an exponentially growing solution:

Resolution  L2-error  L2-order
  199       3.36e-05  nan
  399       7.65e-07  5.45
  799       2.36e-08  5.01
 1599       5.05e-09  2.22
 3199       3.01e-09  0.74
 6399       2.15e-09  0.48
51199       nan       nan

I have 'weird' numbers in the 'Resolution' column, because spatial step is computed as $$\Delta x = \frac{1}{N+1},$$ where $N$ is the resolution from the 'Resolution' column. This way I can be sure that, for example, each grid point in the solution on the coarsest grid matches each 256th grid point in the solution on the finest grid (52000 / 200 = 256). Then the discretization error is measured for grid functions as

$$ e = \left( \frac{1}{N+1} \sum u^{(199)}_{i} - u^{(51999)}_{256i} \right) $$

Data from the table above demonstrate that initially I have the 5th order of accuracy, although I do not get asymptotic convergence to value 5.00. Then, the order of the solution drops significantly, which I can explain with floating-point errors.

However, the real issue is that when I do convergence study for the case of an exponentially decaying solution, I do not get the 5th order of accuracy at all:

Resolution  L2-error  L2-order
199         1.73E-08  nan
399         1.14E-11  10.56
799         4.92E-12   1.22
1599        2.18E-12   1.17
3199        1.00E-12   1.12
6399        5.46E-13   0.87
51199       nan       nan

Initially, I get ridiculous 10th order of accuracy, which immediately degrade to the 1st order. I should note that in this case the solution varies in range 1e-5 to 1e-12 (it decays). Can the small amplitude of the solution be a reason for this strange behavior of the order of accuracy?

So, I have the following questions:

  1. Is self-convergence test reliable? Do I measure discretization error correctly by subtracting two solutions pointwise?

  2. What is the reason for drastically different behavior of convergence for the same solver with two different problems?

  • $\begingroup$ This sort of behavior can be caused by incorrect implementation of the operator split. How did you verify the two schemes? $\endgroup$ Dec 28, 2015 at 19:04
  • $\begingroup$ Try if you can construct a simple example to verify your code first. Some thoughs: 1) if you converge, but to a wrong solution, you will never find out in this way. 2) the rates will inevitably drop at some point if you compare to a discrete solution 3) your errors are pretty low, especially for your second example. make sure you don't have roundoff issues involved or "accidentally" compute an exact solution. $\endgroup$ Dec 28, 2015 at 19:35
  • $\begingroup$ @Biswajit, I don't do any operator splitting here. It's a 1-D problem. I just approximate spatial derivative and compute local source functions. $\endgroup$ Dec 28, 2015 at 23:24

1 Answer 1


For all practical purposes, an L2 error of 1e-9 is zero, due to accumulated roundoff. So, for your first example, everything with more than 1600 resolution is exact. For your second example, everything with better than 400 resolution is exact. At these points, you can no longer expect a significant improvement in accuracy by using more points, and the convergence "rates" you compute make no sense.

To test your code better, you need to find a problem where the error is larger on coarse meshes. (You should make sure that it's the relative, not absolute error that's large, so just multiplying the right hand side, initial and boundary values by 1e6 does not count.)

  • $\begingroup$ It's not clear to me, why the threshold for L2-error must be 1e-9. How then ODE solvers with adaptive step sizes can work to the rel_tol=1e-15? Regarding finding a problem - do you mean that I should use the method of manufactured solutions? $\endgroup$ Dec 28, 2015 at 23:23
  • $\begingroup$ I've assumed that the L2 norm of your solution (or $L_\infty$ norm) is around 1, so round-off would be around 1e-16. But you accrue round-off in many steps, and the larger your linear system becomes, the more you accrue. Whether that accumulated round-off is around 1e-11 or 1e-9 I don't know, and it depends on your implementation, your problem, etc. I just see from your data that you level off at around 1e-9, and so that's where I assume your accumulated round-off is. $\endgroup$ Dec 29, 2015 at 17:08
  • $\begingroup$ As for finding a problem -- there is really nothing inherently wrong in comparing with your finest solution. All I wanted to say is that you should make your problem more complicated for your solver. For example, I can typically approximate a solution of the form $\sin(x)e^{-t}$ with just a very small number of elements and time steps down to accuracy $10^{-8}$ (assuming a domain $[0,1]$ and time interval $[0,1]$, as an example). But I can't do this with a solution of the form $\sin(50x)\sin(50t)$ which requires many more elements and time steps. $\endgroup$ Dec 29, 2015 at 17:10
  • $\begingroup$ Ah, and on the question of ODE solvers: the linear systems to be solved in every time step are tiny, and their size is independent of the number of time steps. You don't accrue any significant round-off there, so you can get down to $10^{-12}$ or $10^{-14}$. $\endgroup$ Dec 29, 2015 at 17:12

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.