Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}]

result

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

share|improve this question
1  
    
I think this is a regression: your code produces the expected answer with no problem in Mathematica version 8. –  Jens Jul 27 at 18:36
1  
@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. –  Michael E2 Jul 27 at 20:34

3 Answers 3

up vote 7 down vote accepted

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.

share|improve this answer
    
So FEM is not the default method, but after selecting it, we can handle discontinuous data, at least under Mathematica 10... –  unlikely Jul 29 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. –  user21 Jul 30 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. –  user21 Jul 30 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.

share|improve this answer

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

share|improve this answer
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. –  Szabolcs Jul 27 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... –  unlikely Jul 27 at 16:59

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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