# Are self-convergence tests reliable?

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?

• This sort of behavior can be caused by incorrect implementation of the operator split. How did you verify the two schemes? – Biswajit Banerjee Dec 28 '15 at 19:04
• 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. – Christian Waluga Dec 28 '15 at 19:35
• @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. – Dmitry Kabanov Dec 28 '15 at 23:24

• 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. – Wolfgang Bangerth Dec 29 '15 at 17:08
• 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. – Wolfgang Bangerth Dec 29 '15 at 17:10
• 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}$. – Wolfgang Bangerth Dec 29 '15 at 17:12