Skip to main content

plotting - How does Plot work?


I considered the following function:


sin[x_] := Module[{},

Print["x=", x];
Sin[x]
]

in Mathematica. Next, I tried to plot it using:


Plot[sin[t], {t, 0, 2 Pi}]

Surprisingly, the first three lines of output are:


x=0.000128356
x=t

x=1.28228*10^-7

Can someone explain this behavior? In this case it doesn't cause a problem, but in my "real" case it does.


Summary


acl's answer below offers, at its very beginning a solution to the specific problem. In very-short, the reason that this x=t appears is hidden somewhere in the way Mathematica evaluates the functions. The answers below provide interesting insight into the way it works.


The interested reader should read all the answers and details below, they are invaluable, although might be behind the reach of some of the readers (like, partially, in my case).



Answer



If the problem is that a symbolic argument is passed, you can avoid it thus:


ClearAll[sin];
sin[x_?NumericQ] := Module[{},

Print[x];
Sin[x]
]

which simply defines sin so that it only matches for numeric arguments.


To see what it does, try sin[3.] and sin[x] and notice that the second evaluates to itself, as the definition above does not match.


You can also see what values of x are being evaluated by ClearAll[sin]; sin[x_] := Module[{}, Sow[x]; Sin[x] ] and then Plot[sin[x], {x, 0, 10}]; // Reap. x now appears.


However,


lst = {};
ClearAll[sin2];

sin2[x_] := Module[{},
AppendTo[lst, x];
Sin[x]
]

and then Plot[sin2[x], {x, 0, 10}]; and then we have no symbols in lst at the end.


EDIT: The explanation for this discrepancy between Sow/Reap and using a list is explained by Leonid in the comments. To test his proposal, I tried using a Bag instead of a list (this is undocumented, see Daniel Lichtblau's description) as follows:


AppendTo[$ContextPath, "Internal`"];
lst = Bag[];
sin3[x_] := Module[{},

StuffBag[lst, x];
Sin[x]
]

followed by Plot[sin3[x], {x, 0, 10}];. We now inspect the contents of the bag by BagPart[lst, All] and observe that there is indeed a symbol x in there.


Presumably it has to do with the way scoping constructs interact with the evaluations performed by AppendTo and StuffBag.


EDIT 2 (by Leonid Shifrin)


We can also demonstrate the same using more usual tools. In particular, instead of a Bag that has its own API, we can use any HoldAll wrapper (just not a list), and then the code for the function itself we need not change at all:


In[51]:= 
ClearAll[h];

SetAttributes[h,HoldAll];
lst=h[];
ClearAll[sin2];
sin2[x_]:=Module[{},AppendTo[lst,x];Sin[x]]

In[58]:=
Plot[sin2[x],{x,0,10}];
lst//Short

Out[59]//Short= h[0.000204286,x,<<1131>>,9.99657]


This clarifies what happens. The x inside List is substituted by the numerical value as a result of evaluation in AppendTo, roughly as follows:


In[60]:= 
Clear[x];
lst = {0.000204,x};
Block[{x = 2.04*10^(-7)},
AppendTo[lst,x]];
lst

Out[63]= {0.000204,2.04*10^-7,2.04*10^-7}


while the HoldAll attribute of h prevents the evaluation from happening (this will be more clear yet if we write AppendTo as lst = Append[lst,x]. It is the evaluation of the r.h.s. (Append), where lst is evaluated and x is substituted by its bound value). For h, x inside it does not evaluate, and is therefore kept symbolic. Similar thing happens with Reap-Sow, although the mechanism Reap-Sow is using to store the results is obviously different (but, whatever it is, it bypasses the main evaluation loop, and that is what matters).


EDIT 3 (acl):


There was a question in the comments as to why the numbers returned by Sow/Reap are not in ascending order. The reason is that Plot apparently uses an adaptive algorithm, in the same spirit as one does in adaptive integration (see en.wikipedia.org/wiki/Adaptive_quadrature, for instance).


Do Plot[sin[x], {x, 0, 10}]; // Reap // Last // Last // ListPlot to see it spend more effort at the turning points:


enter image description here


If you add the option MaxRecursion -> 0 to the Plot command, the algorithm does not subdivide steps that it deems inaccurate and the values are in order:


enter image description here


Maybe it is clearer to do it interactively. Let us play with MaxRecursion and PlotPoints:


ClearAll[sin];

sin[x_?NumericQ] := (Sow[x]; Sin[x])

Manipulate[
pts = ((plt = Plot[
sin[x],
{x, 0, 10},
PlotStyle \[Rule] {Red, Thin},
PlotPoints \[Rule] n,
MaxRecursion \[Rule] m
];) // Reap // Last // Last);


Show[
{
ListPlot[
Transpose@{pts, Sin[pts]},
PlotMarkers \[Rule] {Automatic, 3}
],
plt
}
],

{m, Range[0, 5]},
{{n, 10}, Range[1, 50]}
]

m is the value of MaxRecursion, n that of PlotPoints. The plot shows the resulting plot of Sin and, overlaid, the points that have been evaluated to produce it. Play with the numbers and it should become clear what is happening: PlotPoints tells Plot how many points to evaluate initially, MaxRecursion tells Plot how many times it may subdivide the regions thus defined if necessary (see here for a discussion of what "necessary" means).


enter image description here


Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...