Can NDSolve handle discontinuous data?

Question

Backslide introduced in 9.0, persisting through 11.0.1

---

It is possible to numerically solve a differential equation if not-smooth data are involved?

For example the following instruction return the error NDSolve::bvdisc:

NDSolve[{-u''[x] == UnitBox[(x - 0.5) 0.5/0.2], u[0] == 0, 
  u'[1] == 0}, u, {x, 0, 1}]

NDSolve::bvdisc: NDSolve is not currently able to solve boundary value problems with discrete variables. >>

Here I suppose discrete variable means discontinuous functions but maybe the problem is completely different.

This error is reported in Mathematica 10, where the Finite Element Framework should be able to return weak, not classical, solutions and should be definitely able to deal with such a simple equation, but maybe we need to enable some option or methods.

Please note I know I can do

DSolve[{-u''[x] == UnitBox[(x - 0.5) 0.5/0.2], u[0] == 0, 
  u'[1] == 0}, u, {x, 0, 1}]

but in more complex scenarios I cannot hope to get an analytical answer.

Answer 1

You can force NDSolve to use the finite element method:

uif = NDSolveValue[{-u''[x] == UnitBox[(x - 0.5) 0.5/0.2], u[0] == 0, 
    u'[1] == 0}, u, {x, 0, 1}, 
   Method -> {"PDEDiscretization" -> "FiniteElement"}];

Possibly, this might be the default.

In any case if you compare to the analytical solution:

aif = DSolveValue[{Rationalize[-u''[x] == UnitBox[(x - 0.5) 0.5/0.2]],
    u[0] == 0, u'[1] == 0}, u, {x, 0, 1}];

The result looks good:

Plot[aif[x] - uif[x], {x, 0, 1}]

For a general case one should keep in mind, is that the quality of the solution will depend on the mesh having nodes or internal boundaries at the discontinuity. There is some documentation about this in the SolvingPDEwithFEM tutorial in the variable coefficients section and probably the mesh generation tutorial is also of interest.

Another option is to switch off "DiscontinuityProcessing"

NDSolve[{-u''[x] == UnitBox[(x - 0.5) 0.5/0.2], u[0] == 0, 
  u'[1] == 0}, u, {x, 0, 1}, 
 Method -> {"DiscontinuityProcessing" -> False}]

to get the old behavior.

Answer 2

This is just an extended comment.

I'm not quite certain what's going on as we can easily implement a shooting method manually here.

(Shooting method, in short: parametrize in terms of u'[0] and very the parameter until u'[1] has the desired value.)

fun = ParametricNDSolveValue[{-u''[x] == UnitBox[(x - 0.5) 0.5/0.2], u[0] == 0, u'[0] == ud}, u, {x, 0, 1}, ud]

FindRoot[fun[ud]'[1], {ud, 1}]
(* {ud -> 0.4} *)

Plot[{UnitBox[(x - 0.5) 0.5/0.2], fun[0.4][x]}, {x, 0, 1}]

This gives the solution to you equation.

Answer 3

I've had difficulty before with UnitBox, although I've forgotten the exact context. UnitStep usually works better.

sol = NDSolve[{-u''[x] == UnitStep[x - 0.3] - UnitStep[x - 0.7], 
    u[0] == 0, u'[1] == 0}, u, {x, 0, 1}];

Plot[u[x] /. First[sol], {x, 0, 1}]