# Can NDSolve handle discontinuous data?

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,
u' == 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,
u' == 0}, u, {x, 0, 1}] but in more complex scenarios I cannot hope to get an analytical answer.

• This is relevant: reference.wolfram.com/language/ref/DiscreteVariables.html Jul 27, 2014 at 15:46
• I think this is a regression: your code produces the expected answer with no problem in Mathematica version 8.
– Jens
Jul 27, 2014 at 18:36
• @Jens Could this be because DiscreteVariables were added in V9 -- that is, the error is the result of some advance in NDSolve that improves some area but is unreliable in another (hence the message)? (I do not know.) One can still achieve the V8 behavior via f[x_?NumericQ] := UnitBox[...] as Szabolcs has observed. Jul 27, 2014 at 20:34

You can force NDSolve to use the finite element method:

uif = NDSolveValue[{-u''[x] == UnitBox[(x - 0.5) 0.5/0.2], u == 0,
u' == 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, u' == 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,
u' == 0}, u, {x, 0, 1},
Method -> {"DiscontinuityProcessing" -> False}]


to get the old behavior.

• So FEM is not the default method, but after selecting it, we can handle discontinuous data, at least under Mathematica 10... Jul 29, 2014 at 20:33
• it is the default method in certain cases but not in this one, which may or may not be a bug/regression. Jul 30, 2014 at 5:56
• The fact that this is not the default in this case may or may not be a bug/regression. I'll have to check that. Jul 30, 2014 at 6:07

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' and very the parameter until u' has the desired value.)

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

FindRoot[fun[ud]', {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.

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, u' == 0}, u, {x, 0, 1}];

Plot[u[x] /. First[sol], {x, 0, 1}] • it's also possible to just make a numerical black box out of UnitBox, as in f[x_?NumericQ] := UnitBox[(x - 0.5) 0.5/0.2]. Then it works. But I'd like to understand why Mathematica thinks that it shouldn't attempt solving this equation when it knows that there are variables taking discrete values. Jul 27, 2014 at 15:50
• @Szabolcs I think this approach is the most general we can use to allow a discontinuous function inside an NDSolve, if we don't know in advance what discontinuos functions exactly can be involved (UnitStep, HeavisideTheta, Boole, ...). Thanks. There are some issue? I'm curious too to understand why Mathematica behave this way... Jul 27, 2014 at 16:59