# How to prevent instability blow up in NDSolve?

I have the following code to solve a PDE:

e = 2.5;
xmax = 5;
ymax = 5;
sol[x_, y_] = f[x, y] /. First@NDSolve[{
-D[f[x, y], x, x] - D[f[x, y], y, y] == e f[x, y],
Derivative[0, 1][f][x, -ymax] == Cos[\[Pi]/(2 xmax) x],
f[x, -ymax] == 0,
f[-xmax, y] == 0,
f[xmax, y] == 0
}, f[x, y], {x, -xmax, xmax}, {y, -ymax, ymax}]


I'd expect to get a Sin[]-like solution in y direction, but what I get instead is this:

Plot3D[sol[x, y], {x, -xmax, xmax}, {y, -ymax, ymax}, PlotRange -> {-1, 1}, AxesLabel -> {"x", "y", "f"},  MaxRecursion -> 3]


If I try using something like MaxStepSize -> 0.25, then the blow-up is just more frequent in x direction, and becomes visible at even smaller y values:

What can I do to prevent this blow-up and make NDSolve give the expected solution?

EDIT

In fact, the problem above is a minimal example showing the instability. What I'd actually like to solve is the equation with additional +U[x,y]f[x,y] on the LHS, where U[x,y]=-70Exp[-x^2-y^2]. I.e. the equation would look like

-D[f[x, y],x,x]-D[f[x,y],y,y]-70Exp[-x^2-y^2]f[x,y]==e f[x,y]


So, the answer I'm looking for should be extensible to this case. The answer by @xzczd solves the problem with minimal example, but unfortunately fails to extend to this one.

• Just a check. Do you mean to have both second derivative terms as negative? If you are trying to solve the wave equation then the first one should be positive. If both are negative then the solutions may blow up. Changing the sign of the first term gives a good solution. – Hugh Dec 29 '14 at 13:12
• @Hugh no, the equation is as intended, with both derivatives negative. I do understand that solutions could blow up in some circumstances (e.g. when e<0), but in this case the solution NDSolve gives blows up with the frequency of the grid, so in the limit of infinite grid density the function won't be differentiable, which doesn't look right. – Ruslan Dec 29 '14 at 13:44
• As bbgodfrey pointed out, many finite-difference-based methods actually become less stable when a finer grid is used, rather than more stable (this tends to surprise people, as intuition suggests using a finer grid will usually make things better, not worse). One solution is to use other methods, like Michael E2's suggestion. There might be other ways, too, but I'm not very knowledgeable here. – DumpsterDoofus Dec 30 '14 at 0:10
• Just a side note, today I came across this article and noticed your problem is actually ill-posed. Still, a quick search shows there seems to be several techniques to deal with it, but sadly I can't understand those materials I found… – xzczd Sep 25 '16 at 7:09

The Finite Element Method seems more stable for this type of problem.

e = 2.5;
xmax = 5;
ymax = 5;
sol = First@
NDSolve[{-D[f[x, y], x, x] - D[f[x, y], y, y] + NeumannValue[Cos[\[Pi]/(2 xmax) x], y == -ymax] == e f[x, y],
DirichletCondition[f[x, y] == 0, y > -ymax]},
f, {x, -xmax, xmax}, {y, -ymax, ymax},
Method -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> 0.01}}];

Plot3D[f[x, y] /. sol, {x, -xmax, xmax}, {y, -ymax, ymax},
PlotRange -> {-1, 1}, AxesLabel -> {"x", "y", "f"}, PlotPoints -> 101]


We can use "MeshOptions" to control the quality of the solution (as described in the linked tutorial above). Comparing the boundary condition along y == -ymax, we see that MaxCellMeasure -> 0.01 helps to meet the condition much better than the default of 0.25:

Plot[{Cos[π/(2 xmax) x], D[f[x, y], y] /. sol /. y -> -ymax} // Evaluate, {x, -xmax, xmax}]


I hope it works in your actual use case! :)

• This is a great answer. Would you know why the default method does not check its convergence and falls back onto this solution? It makes using NDSolve hazardous that no such validation is carried out? – chris Dec 29 '14 at 19:38
• @MichaelE2 unfortunately it requires Mathematica 10. Still interesting anyway. What would be the way to add U[x,y]f[x,y] to the LHS of the equation with this awkward Neumann condition syntax? – Ruslan Dec 29 '14 at 19:46
• @Ruslan I think there might be a way to do FEM in V9....I'm not sure how much work it would take. See for instance Jens' answer here. – Michael E2 Dec 29 '14 at 20:20
• @chris Thanks. I don't know how NDSolve chooses the method. It may have to do with the OP's setup. Someone from WRI might see your comment and be able to respond. – Michael E2 Dec 29 '14 at 20:24
• It seems your formulation of the problem is not the same as the one in the OP. See wave equation section in the tutorial you linked. The problem in the OP basically only differs by a sign in a derivative (and other minor parts). But there's no NeumannValue[] used in the problem setup in that tutorial. And when I try to add a U[x,y]f[x,y] term to the LHS of the equation with your formulation, I get completely wrong results, which don't even satisfy initial conditions. (I used U[x,y]=-70Exp[-x^2-y^2]). Are you sure your formulation is equivalent to mine? – Ruslan Dec 30 '14 at 7:33

Your calculation is unstable, because the Courant limit is violated. For a square mesh, it is e > 4 dx^-2. Choosing e=27.5 gives

• So, what are the workarounds? Maybe some way to use rectangular mesh or some special Method of NDSolve? – Ruslan Dec 29 '14 at 14:20
• The equation can be solved analytically, although I have not been able to coax DSolve to do it yet.. – bbgodfrey Dec 29 '14 at 14:28
• Of course it can! But the equations I'm going to solve after fixing the stability problem can't! I just gave a minimal working example. – Ruslan Dec 29 '14 at 14:30
• What is actual equation you are trying to solve? Courant limits sometimes are sensitive to details. – bbgodfrey Dec 29 '14 at 14:32
• Add U[x,y]f[x,y] on the LHS of the equation where U[x,y] is some function, and e>U[x,y], then it's the equation I'm trying to solve. Restricting values of e beyond this is not a solution, of course. – Ruslan Dec 29 '14 at 14:33

Somewhat inspired by bbgodfrey's answer, after rereading the tutorial of "MethodOfLines", I found that a restriction for "MaxPoints" together with "DifferenceOrder" -> 2 solved the problem, no v10-feature is needed!:

e = 2.5;
xmax = 5;
ymax = 5;
sol[x_, y_] =
f[x, y] /.
With[{n = 18},
First@NDSolve[{-D[f[x, y], x, x] - D[f[x, y], y, y] == e f[x, y],
Derivative[0, 1][f][x, -ymax] == Cos[π/(2 xmax) x],
f[x, -ymax] == 0, f[-xmax, y] == 0, f[xmax, y] == 0},
f[x, y], {x, -xmax, xmax}, {y, -ymax, ymax},
InterpolationOrder -> All,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> n, "MinPoints" -> n,
"DifferenceOrder" -> 2}}]]

Plot3D[sol[x, y], {x, -xmax, xmax}, {y, -ymax, ymax}, AxesLabel -> {"x", "y", "f"}]
Plot[{Cos[π/(2 xmax) x] - D[sol[x, y], y] /. y -> -ymax} // Evaluate, {x, -xmax, xmax}]


Notice that "MinPoints" option isn't necessary here, but I found it lowers the error, at least at the boundary. And actually if you choose n = 10, "DifferenceOrder -> 2" can be eliminated, too, though a warning will be generated.

• Great! But I can't seem to extend this to the equation with +U[x,y]f[x,y] on the LHS, where U[x,y]=-70Exp[-x^2-y^2]. I try to increase n, increase difference order, but still the blow up appears at about y==0 or earlier. – Ruslan Dec 31 '14 at 8:32
• @Ruslan If this method would have any chance to work, it should be decreasing rather than increasing n or "DifferenceOrder", sadly it seems that decreasing doesn't work either 囧. – xzczd Dec 31 '14 at 14:29