# Debugging an implemented numerical method: which term gives the drop in accuracy?

I have a numerical scheme to solve an hyperbolic system of equations of this type (this a more simple version just for clarity):

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\frac{\partial H}{\partial x}$$ $$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial x}=\frac{\partial H}{\partial y}$$ $$\frac{\partial H}{\partial t}+ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial x}=0$$

I have a steady state analytical solution, and I am checking the order of accuracy of the scheme. It is a direction splitting scheme, and all the terms should be second order accurate. I used 4 grids, with successive refinements $i=1\ldots 4$ (cartesian uniform grids with grid size $\Delta x_i=\Delta x_{i-1}/2$ and $\Delta x=\Delta y$). I saw that the order of accuracy is somehow around 1.5, therefore I decided to replace each derivative with the correspondent analytical term until I found the term responsible for the loss of accuracy. The culprits were: $\frac{\partial H}{\partial x}$ and $\frac{\partial H}{\partial y}$. Now my questions are:

1) Am I doing right by using a fixed CFL number, in the sense that every time I half the grid size I also half the time step? What if I use a fixed dt for all the grids? Note that the scheme is implicit and unconditionally stable.

2) Im using a staggered grid, i.e. H is cell centered and u and v are on the center of the cell edges. The terms $\frac{\partial H}{\partial x}$ and $\frac{\partial H}{\partial y}$ are discretized by second order central differences. When I print the error of $\frac{\partial H}{\partial x}$ after a single time step starting from the exact solution I see that it converges perfectly with second order accuracy (it is 4 times smaller at each refinement). However when I start from a second-order accurate approximation of H, $\frac{\partial H}{\partial x}$ has first order rate of convergence. I am trying to figure this out and it might make sense since $\frac{\partial H}{\partial x}= \frac{H_{m+1}-H_m}{\Delta x} + O(\Delta x^2)$ that is a second order accurate expression at the edge $m+1/2$ if $H_{m+1}$ and $H_m$ are exact. However, if I start from a second order accurate approximation of $H_{m+1}$ and $H_m$: $$H_{m+1} = H^{exact}_{m+1}+ O(\Delta x^2)$$ $$H_{m} = H^{exact}_{m}+ O(\Delta x^2)$$ (as it would be the second order accurate solution of the system) I get $$\frac{\partial H}{\partial x}= \frac{H^{exact}_{m+1}-H^{exact}_m}{\Delta x} + \frac{O(\Delta x^2)}{\Delta x} + O(\Delta x^2)$$ that is first order, and thats what I get (first order) if I print the numeric value of $\frac{\partial H}{\partial x}$ for the 4 different grids. However, if I print $H$ for the first half time step it is second order (it is the second half of the splitting that cause the order of accuracy). Am I doing something wrong here?

• This type of convergence error can show up if you're not handling your boundary conditions properly. What are your BCs, and how are they implemented? If you have 1st order approximations to $\nabla H$ at the boundary, this can kill the 2nd order convergence for the whole problem. – Tyler Olsen Jul 21 '15 at 14:45
• My boundary conditions are second order accurate. At least theoretically. I d love if somebody could kindly answer to point (1) and (2). – Millemila Jul 21 '15 at 16:24
• Deterministic equations should usually get integer orders (only methods on stochastic equations should get half-order convergences). You may need to check the convergence at a much lower $\Delta x$. But when it gets too small, you might have floating point errors causing inaccuracies as well. – Chris Rackauckas Feb 22 '17 at 19:37

1) This depends. If your analytical solution is a steady solution, then when you reach that solution, your time derivative term should be close to zero (or close to the machine error at least). Thus keeping a constant CFL condition is not necessary. Since your scheme is implicit and unconditionally stable, you do not need to worry about the influence of the time step of the error. For some solution we investigate, you will sometimes have an analytical solution that depends on time. In this case, your error, if you have a second order time discretization, will be: $$e \propto O(\Delta x^2) + O(\Delta t^2)$$
Thus if you keep the CFL constant $$e \propto O(\Delta x^2) + O (CFL \Delta x^2) \propto O(\Delta x^2)$$