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.

Yesterday I had an issue that kept for quite a while. My code wasn't producing the right results. Now I know why and I am wondering how to avoid the issue in future. Consider the following code.

Clear[fAux];
fAux[i_] := fAux[i] = NDSolve[{y'[x] == Sin[y[x]^i], y[0] == 1}, y, {x, 0, Pi}][[1, 1, 2]];
f[i_, z_] := fAux[i][z];

Here, f returns the value of the solution of a ODE with parameter i at the value z. The solution of the ODE (an interpolating function) is memoized. You can convince yourself that it works fine by evaluating f[1,1.5] for example.

Now I want to evalute this:

Table[f[2, x], {x, 0, 3}]

and I get an error message

"NDSolve::dsvar: 0 cannot be used as a variable."

as well as unexpected results:

{Sin[y[0]^2][0], Sin[y[0]^2][1], Sin[y[0]^2][2], Sin[y[0]^2][3]}

Looking at the down values of fAux you can see some strange things:

DownValues[fAux]

{HoldPattern[fAux[1]] :> InterpolatingFunction[{{0.`, 3.141592653589793`}},"<>"],
HoldPattern[fAux[2]] :> Sin[y[0]^2],
HoldPattern[fAux[i_]] :> (fAux[i] = NDSolve[{(y^\[Prime])[x] == Sin[y[x]^i], y[0] == 1}, y, {x, 0, \[Pi]}][[1,1,2]])}

The culprit is the name of the iterator in the Table function. It has the same name, x, as the dummy variable in NDSolve. If you replace either with xx for instance, everything works as expected. I find this very confusing and dangerous, to keep track of what name you gave to which iterator variable throughout the entire notebook. I'm not even sure how this would behave across packages you call from within the notebook. Is this a wanted behavior? How can I make sure I won't make the same mistake again?

share|improve this question
add comment

3 Answers

up vote 8 down vote accepted

It may be helpful to understand the cause of the problem.

  1. Table acts much like Block in the way that it applies the values of its variables. Because of this it will affect things beyond explicit appearances of the variable in the body of Table.

  2. NDSolve, despite having syntax highlighting that indicates that x is localized, does in fact not have the Attributes required for this to take place.

You can observe that setting HoldAll on NDSolve fixes this problem, but I am not sure what doing so may break:

SetAttributes[NDSolve, HoldAll]

Table[f[2, x], {x, 0, 3}]
{1., 1.71208, 1.77067, 1.7724}

While it is arguably a pathological use of Formal Symbols in Table it is also worth noting that Formal Symbols do not protect against the Block-like scoping of Table:

ClearAttributes[NDSolve, HoldAll]
Clear[fAux];
fAux[i_] := 
  fAux[i] = NDSolve[{y'[\[FormalX]] == Sin[y[\[FormalX]]^i], y[0] == 1}, 
     y, {\[FormalX], 0, Pi}][[1, 1, 2]];
f[i_, z_] := fAux[i][z];

Table[f[2, \[FormalX]], {\[FormalX], 0, 3}]

NDSolve::dsvar: 0 cannot be used as a variable. >>

{Sin[y[0]^2][0], Sin[y[0]^2][1], Sin[y[0]^2][2], Sin[y[0]^2][3]}

Since globally setting HoldAll on NDSolve is probably inadvisable you could could manually localize x with Module while the definition is made:

Module[{x},
 fAux[i_] := 
   fAux[i] = NDSolve[{y'[x] == Sin[y[x]^i], y[0] == 1}, y, {x, 0, Pi}][[1, 1, 2]];
]

This calls Module only once at the time of definition rather than every time fAux is used.

share|improve this answer
    
Is there a technical reason for the different behavior of Table and NDSolve in this case? Having to remember to put Module around certain function definitions smells like an ugly workaround to me, and I wish Mathematica would just behave as I expect it in this case. –  Volker Apr 18 '13 at 9:55
    
@Volker Well, the technical reason is the lack of HoldAll on the NDSolve function. As to why it lacks that Attribute I cannot say; I have not adequately considered the situation. It almost seems like a bug given the fact that the syntax highlighter incorrectly marks that symbol as localized when in fact it is not. Sadly, existing code may require the non-held behavior of NDSolve therefore it is not a reasonable fix to simply apply that attribute. You could perhaps write a proxy function e.g. myNDSolve that does localize correctly; should I add an example of that? –  Mr.Wizard Apr 18 '13 at 10:03
    
You don't need to add an example for me, I think I'll manage. Thanks! –  Volker Apr 18 '13 at 10:05
    
@Volker Okay. Sorry I don't have a more satisfying solution for you. –  Mr.Wizard Apr 18 '13 at 10:06
add comment

You can use formal symbols instead for safety, like so:

fAux[i_] := fAux[i] = \[FormalY] /. First[NDSolve[
   {\[FormalY]'[\[FormalX]] == Sin[\[FormalY][\[FormalX]]^i], \[FormalY][0] == 1},
    \[FormalY], {\[FormalX], 0, Pi}]];
f[i_, z_] := fAux[i][z];

These formal symbols will not clash with your uses of normal x and y, so these should be considered. An alternative is to localize x and y with a scoping construct like Block[] or Module[]:

fAux[i_] := fAux[i] = Module[{x, y}, y /. First @
                      NDSolve[{y'[x] == Sin[y[x]^i], y[0] == 1}, y, {x, 0, Pi}]]

which was how it was done before the formal symbols came along.

share|improve this answer
    
The formal symbols suggestion only works if I use formal symbols in either. If I replace x by \[FormalX] in both the Table and the NDSolve command, I'm back at square 1. –  Volker Apr 18 '13 at 9:49
    
@Volker Correct. See my answer below. –  Mr.Wizard Apr 18 '13 at 9:50
add comment

Placing your function in a context or fully packaging it will turn your inner x into SomeContext`x and prevent the collision. It is then only a problem when you're writing or changing the package when you can take steps to avoid using the same symbol names.

Begin["SomeContext`"];
fAux[i_]:=fAux[i]=NDSolve[{y'[x]==Sin[y[x]^i],y[0]==1},y,{x,0,Pi}][[1,1,2]];
f[i_,z_]:=fAux[i][z];
End[];

Table[SomeContext`f[2,x],{x,0,3}]
{1.,1.71208,1.77067,1.7724}

This is presumably why such collisions do not happen with standard package functions all the time.

share|improve this answer
add comment

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.