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.

  • 2
    $\begingroup$ This is relevant: reference.wolfram.com/language/ref/DiscreteVariables.html $\endgroup$ – Szabolcs Jul 27 '14 at 15:46
  • $\begingroup$ I think this is a regression: your code produces the expected answer with no problem in Mathematica version 8. $\endgroup$ – Jens Jul 27 '14 at 18:36
  • 1
    $\begingroup$ @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. $\endgroup$ – Michael E2 Jul 27 '14 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] == 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}]

Mathematica graphics

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.

  • $\begingroup$ So FEM is not the default method, but after selecting it, we can handle discontinuous data, at least under Mathematica 10... $\endgroup$ – unlikely Jul 29 '14 at 20:33
  • $\begingroup$ it is the default method in certain cases but not in this one, which may or may not be a bug/regression. $\endgroup$ – user21 Jul 30 '14 at 5:56
  • $\begingroup$ 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. $\endgroup$ – user21 Jul 30 '14 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'[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}]

enter image description here

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

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

Mathematica graphics

  • 1
    $\begingroup$ 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. $\endgroup$ – Szabolcs Jul 27 '14 at 15:50
  • $\begingroup$ @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... $\endgroup$ – unlikely Jul 27 '14 at 16:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.