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.

I was trying to integrate a continuous function with a kink and I did it two ways and both ways the plot of the result shows a discontinuity. I also later want to differentiate the Integrated function.

a=5;
DSolve[Y'[H] == Max[H, a], Y[H], H]
{{Y[H] -> C[1] + 1/2 (2 a H + (a - H)^2 UnitStep[-a + H])}}

or alternatively

Integrate[Max[H, a], H]

I get

$ \begin{array}{ll} \Big\{ & \begin{array}{ll} 5 H & H\leq 5 \\ \frac{25}{2}+\frac{H^2}{2} & \text{True} \end{array} \end{array} $

When I plot the output from the ODE or Plot what I get from the result of the Integration command, I get the same plot, of course, but with this discontinuity at $a$.

Plot[1/2 (2 a H + (a - H)^2 UnitStep[-a + H]), {H, 4.9, 5.1}]

I get:

Plot with discontinuity

Why is this happening? If I do other operations with this function, is Mathematica gonna treat it as a continuous function or a discontinuous one? I can't remember my math right now but this new function should be differentiable also right since it was the result of an Integration?

When I do this:

 FullSimplify[D[1/2 (2 a H + (a - H)^2 UnitStep[-a + H]), H]]

I get:

$ \begin{array}{ll} \{ & \begin{array}{ll} 5 & H<5 \\ H & H>5 \\ \text{Indeterminate} & \text{True} \end{array} \end{array} $

Is the fact that I am getting this Indeterminate in the middle something to do with the discontinuity?

share|improve this question
    
Try e.g. Plot[1/2 (2 a H + (a - H)^2 UnitStep[-a + H]), {H, 4.9, 5.1}, PlotPoints -> 300] –  Artes Oct 26 '12 at 0:56
    
@Artes ok that worked. Thanks! Do you have a sense of why I am getting an Intermediate when I differentiate? –  Amatya Oct 26 '12 at 1:08
    
Try a simpler case, e.g. : D[ UnitStep[x], x], UnitStep is a discontinuous function at 0. –  Artes Oct 26 '12 at 1:14
    
@Artes but 1/2 (2 a x + (a - x)^2 UnitStep[-a + x]) is a continuous and even differentiable function whose derivative should be Max[x,a]. I guess MMA is not able to figure this out and I need to communicate it to MMA somehow. –  Amatya Oct 26 '12 at 1:31
4  
The gap in the plot can also be removed using the option setting Exclusions->None. If you plot PiecewiseExpand[1/2 (2 a h + (a - h)^2 UnitStep[-a + h])] you also get the same gap at h=5 although D[PiecewiseExpand[1/2 (2 a h + (a - h)^2 UnitStep[-a + h])],h] does not have the Indeterminate piece; but Plot "sees" some discontinuity in both functions, and inreasing the setting for PlotPoints reduces the size of gap, and Exclusions->None removes it. (Mma V 8.04.0, Windows Vista 64bit) –  kguler Oct 26 '12 at 1:36

1 Answer 1

up vote 15 down vote accepted

The solution is to use Exclusions->None as option to Plot.

The gap happens exactly where UnitStep[-a+h] has its discontinuity

With[{a = 5}, 
  Plot[{1/2 (2 a H + (a - H)^2 UnitStep[-a + H]), 
    UnitStep[-a + H] + 25}, {H, 4.9, 5.1}]
 ]

Mathematica graphics

This behavior was introduced, when Wolfram decided, that discontinuities should be discontinuous displayed in Plot. When you look at the function, to see whether or not there is a crack, you should use Limit. Here you see, that the derivative is the same from both directions

With[{a = 5},
   Limit[D[1/2 (2 a H + (a - H)^2 UnitStep[-a + H]), H], H -> 5, 
    Direction -> #]
   ] & /@ {1, -1}  

(* {5, 5} *)

Therefore, it seems Plot internals work as if there is a discontinuity where there is none. The first two examples in the following suggest that forcing Plot to sample more points does not remove the "discontinuity". The third one (Exclusions->None) explicitly forces Plot to treat the line as continuous, and the last one (I think ... somehow) makes Plot internals "realize" that there are, in fact, no discontinuities in the first argument.

 Grid[{{Plot[1/2 (2 a h + (a - h)^2 UnitStep[-a + h]), {h, 4., 6},
        PlotPoints -> 50, ImageSize -> 300, 
        ExclusionsStyle -> Directive[AbsoluteThickness[5], Red],
        PlotLabel -> HoldForm[PlotPoints -> 50]],
   Plot[1/2 (2 a h + (a - h)^2 UnitStep[-a + h]), {h, 4., 6},
        PlotPoints -> 800, ImageSize -> 300, 
        ExclusionsStyle -> Directive[AbsoluteThickness[5], Red],
        PlotLabel -> HoldForm[PlotPoints -> 800]]},
  {Plot[1/2 (2 a h + (a - h)^2 UnitStep[-a + h]), {h, 4., 6},
       PlotPoints -> 10, ImageSize -> 300, Exclusions -> None,
       PlotLabel -> HoldForm[{Exclusions -> None, PlotPoints -> 10}]],
   Plot[1/2 (2 a hh + (a - hh)^2 UnitStep[-a + hh]) /. hh -> h, {h, 4., 6},
       PlotPoints -> 10, ImageSize -> 300, 
       ExclusionsStyle -> Directive[AbsoluteThickness[5], Red],
       PlotLabel -> HoldForm[{Plot[f[x] /. x -> h, _], PlotPoints -> 10}]]}}]

enter image description here

share|improve this answer
    
+1, it's not a comment ;-) Please review my additions. –  halirutan Oct 26 '12 at 8:57
    
@halirutan, very nice, thank you. –  kguler Oct 26 '12 at 13:32
    
@kguler Thank for for such an excellent explanation! I learnt about using limits in mathematica, especially inside a plot, as an additional bonus. :) It's interesting that you spotted that taking limits makes mathematica realize that the curve is continuous. cool beans! –  Amatya Nov 4 '12 at 20:46
    
@Amatya, thanks for the accept. –  kguler Nov 4 '12 at 23:29

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.