6
$\begingroup$

I want a plot with this function

-((H^3 T)/(200 Pi)) + (H^2 T^2)/96 + H^4 l + H^2 m 

with

rule = {l -> 27/(40000 Pi^2), m -> -85, T -> 106}

So if I give as input

f = -((H^3 T)/(200 Pi)) + (H^2 T^2)/96 + H^4 l + 
  H^2 m   /. {l -> 27/(40000 Pi^2), m -> -85, T -> 106}

and plot

Plot[Log[Abs[f]] Sign[f], {H, 0, 4000}, PlotRange -> All]

the plot line is not continuous.

Instead, if I didn't set f and I explicitly write

Plot[Log[Abs[-((H^3 T)/(200 Pi)) + (H^2 T^2)/96 + H^4 l + 
  H^2 m   /. {l -> 27/(40000 Pi^2), m -> -85, T -> 106}
    ]] Sign[-((H^3 T)/(200 Pi)) + (H^2 T^2)/96 + H^4 l + 
  H^2 m   /. {l -> 27/(40000 Pi^2), m -> -85, T -> 106}
   ], {H, 0, 4000}, PlotRange -> All]

the plot line is continuous.

Why there is this difference? Is there a way to obtain a continuous line labelling the function as f?

$\endgroup$
  • 1
    $\begingroup$ Plot[Log[Abs[f]] Sign[f], {H, 0, 4000}, PlotRange -> All, Exclusions -> None] $\endgroup$ – corey979 Oct 9 '16 at 8:14
  • 1
    $\begingroup$ I do think it's still relevant to know why the difference in exclusions. $\endgroup$ – Feyre Oct 9 '16 at 8:54
  • $\begingroup$ @Feyre Indeed it's intriguing, but I'm kinda stumped. With f as in the OP (so with /.rule already) consider this: Plot[Log[Abs[f]] Sign[f] /. rule, {H, 0, 4000}, PlotRange -> All] (so I basically used the same rule twice!). $\endgroup$ – corey979 Oct 9 '16 at 9:50
6
$\begingroup$

Why the difference arises

The difference arises because Log[f] (as well as Sign[f]) is discontinuous where f == 0, which Plot computes symbolically and gets different results in the two Plot codes. There should be vertical asymptotes at those points. (So Exclusions -> None is not a fix, as it gives a continuous graph.)

Here are the exclusions computed by Plot in each of the OP's cases:

Visualization`ExpandExclusions[Log[f] Sign[f], {H}, Automatic]
(*
  {{Im[(769 H^2)/24 + (27 H^4)/(40000 π^2) - (53 H^3)/(100 π)] == 0, 
    Re[(769 H^2)/24 + (27 H^4)/(40000 π^2) - (53 H^3)/(100 π)] <= 0},
   {(769 H^2)/24 + (27 H^4)/(40000 π^2) - (53 H^3)/(100 π) == 0, True},
   {(769 H^2)/24 + (27 H^4)/(40000 π^2) - (53 H^3)/(100 π) == 0, True}}
*)

Visualization`ExpandExclusions[
 Log[Abs[-((H^3 T)/(200 Pi)) + (H^2 T^2)/96 + H^4 l + H^2 m /.
     {l -> 27/(40000 Pi^2), m -> -85, T -> 106}]] *
  Sign[-((H^3 T)/(200 Pi)) + (H^2 T^2)/96 + H^4 l + H^2 m /.
     {l -> 27/(40000 Pi^2), m -> -85, T -> 106}],
 {H}, Automatic]
(*
  {}
*)

One can see these calls with

Trace[
 (* plot command *),
 _Visualization`ExpandExclusions,
 TraceInternal -> True]

Update: The reason no exclusions are produced is that the unevaluated expression being plotted contains the member ReplaceAll of the following blacklist, which prevents further analysis of potential exclusions:

Visualization`DiscontinuityDump`$BlackList
(*
  {CompiledFunction, InterpolatingFunction, LinearSolveFunction, 
   NearestFunction, TransformationFunction, NIntegrate, NSum, NDSolve, 
   FindRoot, FindMinimum, FindMaximum, NMinimize, NMaximize, FixedPoint, 
   FixedPointList, Nest, NestList, Fold, FoldList, NestWhile, 
   NestWhileList, Apply, Map, Table, Do, For, While, Set, SetDelayed, 
   Decrement, PreDecrement, Increment, PreIncrement, Rule, RuleDelayed, 
   ReplaceAll, ReplaceRepeated, Replace, Nearest}
*)

Getting an accurate plot

The order of growth of f at the asymptotes is very small, and to get large enough values of f that look like asymptotes, it takes values of H that are closer to the singularities than MachinePrecision will allow.

For instance, to get a magnitude greater than 40, we have to be within roughly 10^-24 of the singularities:

Block[ (* H near singularities *)
 {H = H + {-1*^-24, 1*^-24} /. 
    NSolve[f == 0 && 0 < H < 4000, H, WorkingPrecision -> 32]},
 Log[Abs[f]] Sign[f]
 ]
(*  {{-46.5568, 46.5568}, {41.780, -41.780}}  *)

It is impractical to do this with Plot, since it automatically controls selection of the sample points. I will use ListLinePlot instead. If we include the points above, as well as one close to H == 0 where there is another asymptote, we can get this picture:

ListLinePlot[
 Table[{H, Log[Abs[f]] Sign[f]},
   {H, Union[Range[1.`32*^-20, 4000, 5], 
        Flatten[H + {-1*^-20, 1*^-20} /. 
         NSolve[f == 0 && 0 < H < 4000, H, WorkingPrecision -> 32]]]}], 
 PlotRange -> {-25, 25}]

Mathematica graphics

It's virtually impossible to mark up the asymptotes because the graph is so close to them. (One can try dashing or transparency, but the result is not very good, imo.)

$\endgroup$
  • $\begingroup$ Excellent answer (+1). Does documentation exist for 21 forms of Visualization? $\endgroup$ – bbgodfrey Oct 9 '16 at 14:37
  • 3
    $\begingroup$ @bbgodfrey Thanks! No, not that I know of, and most seem not to have a definition. They may be internal option names or constants set at runtime. In the case of LegendReap, I would guess it might be used to set up code in which it is changed to Reap at the appropriate execution time with Block[]. I've used ExpandExclusions on this site before. The FindDate* symbols seem to be the only other functions, and they might be useful if you're plotting dates. $\endgroup$ – Michael E2 Oct 9 '16 at 14:48
  • $\begingroup$ Hi, thanks for your answer! $\endgroup$ – Giancarlo Oct 11 '16 at 1:48
  • $\begingroup$ Hi, thanks for your answer! It's still not clear to me why Plot gets different results in the two Plot codes. You said that is because Plot computes computes symbolically the Log. But it seems to me that Plot should compute in the same way both expressions; so why the different output? $\endgroup$ – Giancarlo Oct 11 '16 at 1:54
  • 1
    $\begingroup$ @Giancarlo Plot has the attribute HoldAll, which means, at least in this case, that it symbolically analyzes the expression in the form you pass to Plot. The first form gets analyzed and the graph is separated at the discontinuities. Because of the /. (= ReplaceAll) in the second form, Plot automatically does not analyze the discontinuities (because of the blacklist). As a result, it connects the pieces. $\endgroup$ – Michael E2 Oct 11 '16 at 2:00

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.