# 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? 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. 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. 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. 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. 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}$. Dec 29 '15 at 17:12