Plot showing discontinuity where it shouldn't

Question

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:

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?

Answer 1

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

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

Answer 2

In newer versions (V11 or at least by V11.3), the gap shown at a point discontinuity is almost negligible, but present:

Plot[0.1 UnitStep[1/200 - a + H] + 
   1/2 (2 a H + (a - H)^2 UnitStep[-a + H]), {H, 4.9, 5.1}];
GraphicsRow[{%, Show[%, PlotRange -> {{4.995, 5.002}, Automatic}]}, ImageSize->Large]

Here's a general utility that removes exclusions that are continuities, but includes actual discontinuities (Exclusions -> None connects all segments, even the truly discontinuous ones):

ClearAll[includeContinuousSingularities];
SetAttributes[includeContinuousSingularities, HoldAll];
includeContinuousSingularities[tol_ /; NumericQ[tol]] := 
  Function[p, includeContinuousSingularities[p, tol], HoldAll];
includeContinuousSingularities[Plot[f_, {v_, a_, b_}, rest___], tol_: 1.*^-10] :=
 Module[{sing, disc},
  sing = Flatten@ Apply[And, DeleteCases[
    Simplify`FunctionSingularities[f, {v}, {"ALL", "IGNORE", "PWMINMAX"}], {}], {2}];
  sing = DeleteDuplicates@ Flatten@ Quiet@
       DeleteCases[NSolve[# && a < v < b, v] & /@ sing, _NSolve | {} | {{}}];
  disc = Equal @@@ Pick[sing,
     Chop[Block[{v = v /. #}, f] - Limit[f, #], tol] =!= 0 & /@ sing];
  Plot[f, {v, a, b}, rest, Exclusions -> disc]
  ]

Plot[
 0.1 UnitStep[1/200 - a + H] + 1/2 (2 a H + (a - H)^2 UnitStep[-a + H]),
 {H, 4.9, 5.1}] // includeContinuousSingularities;
GraphicsRow[{%, Show[%, PlotRange -> {{4.995, 5.002}, Automatic}]}, ImageSize->Large]

One could use Solve instead of NSolve, but NSolve over a bounded interval is pretty robust. One might also put a time constraint (TimeConstrained) on NSolve. With numerical methods, one has to be careful of comparing results; hence, I used Chop. The tolerance parameter tol is probably unnecessary, unless the function values are really small. See also Plot a piecewise function with black and white disks marking discontinuities.