# 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
Commented Dec 29, 2014 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. Commented Dec 29, 2014 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. Commented Dec 30, 2014 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… Commented Sep 25, 2016 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? Commented Dec 29, 2014 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? Commented Dec 29, 2014 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. Commented Dec 29, 2014 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. Commented Dec 29, 2014 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? Commented Dec 30, 2014 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? Commented Dec 29, 2014 at 14:20
• The equation can be solved analytically, although I have not been able to coax DSolve to do it yet.. Commented Dec 29, 2014 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. Commented Dec 29, 2014 at 14:30
• What is actual equation you are trying to solve? Courant limits sometimes are sensitive to details. Commented Dec 29, 2014 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. Commented Dec 29, 2014 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. Commented Dec 31, 2014 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 囧. Commented Dec 31, 2014 at 14:29