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

Comments
Post a Comment