# Tips on improving stability in numerical scheme for non-linear PDE

I am solving a non-linear second order system of PDEs in two variables. The equations are too complicated to write out here, but an essential feature is that there is a propagating wave which then bounces on a boundary.

The problem I have is that the numerics breaks down at the boundary, when the wave reaches this point. I have tried by "trial and error" to just change the way I compute derivatives at this point with different finite difference stenils, but with no luck (pseudospectral methods do not work btw), and I have no idea on how to systematically try to improve the stability and I have no intuition of what can work and what will not work.

Does anyone have any tips on what to try when one encounters such problems? How can I systematically move forward to try to make my numerical scheme stable at this boundary point? Basically, if someone just has a list of ideas to blindly try, that would also be great.

Edit: I should add that the "boundary" can also be interpreted as an origin in polar coordinates, so close to the boundary we have a cylindrical wave scattering at the origin.

Edit2: By request I have added a (reduced version) of the equations I am trying to solve

$\dot{c}=K'$

$\dot{K}=\frac{K^2f}{5r^2}+25(1-f)+\frac{5r}{4}(fc/r)'$

$\dot{f}=-\frac{f^2K'}{5r}$

The "polar coordinate origin" is at r=0, and the wave that propagates is in $c$ and $K$. I have removed many terms to obtain the above equations, and only kept the ones that I think are important. It seems that the instability is due to the $1-f$ term in $\dot{K}$ (if I remove it I get stable evolutions). The instability manifests itself by having $f$ and $K$ diverging with opposite signs at the points close to the origin. Around $r\sim 0$ we expect $K\sim r^4$ and $c\sim r^3$ and $f-1\sim r^2$.

The boundary conditions I put are $K=e^{-t^2}$ at some $r=r_0$, which induces a wave that should bounce at $r=0$. The initial conditions are $f=1$, $c=K=0$.

Some more details about my numerics: I have K and c on different grids (such that a point in $c$ lies between two points in $K$). This removed other instability issues I had. I use second/third order finite difference methods for computing derivatives and for interpolating between the grids. $c$ has a point at the origin, $K$ does not, and $f$ is on the same grid as $K$. At the origin I use the fact that the functions should be even when computing derivatives (by extending the radial coordinate to negative values), but this does not change much.

• Mehta method are you using to solve and what are the boundary conditions? Also knowing the actual equations would help. – boyfarrell Apr 30 '15 at 12:46
• I am using finite difference approximations for spatial derivatives, and an external library for the time evolution (the ode solver in scipy, which uses some backward differentation formula). Well, actually the "boundary" can be seen as the origin in polar coordinates, so it is sort of a cylindrical wave scattering at the origin, so the boundary condition is just regularity. – Jonathan Lindgren Apr 30 '15 at 13:47
• Coordinate singularities (like at the origin of a polar grid) can cause all sorts of weird behaviour. I agree with @boyfarrell, giving your equations and detailing what happens would really help. If the full set is too complicated, just give the relevant terms perhaps? Or name them so we can look them up at least. – Geoff Ryan Apr 30 '15 at 23:14
• @GeoffRyan Ok I have added some more details about the equations. – Jonathan Lindgren May 1 '15 at 23:39
• Sorry for not being clear, the only variables are time and r. All functions (c, K, f) depend on these two variables. Dot is derivative wrt time, prime is derivative wrt r – Jonathan Lindgren May 3 '15 at 20:50