# HoldAll Attributes and Delayed Evaluation: ContourPlot with solutions to NDSolve

This is a general question on evaluating functions that are defined to be evaluated with delay.

In particular, I am having difficulties working with the InterpolationFunction from NDSolve when I define the family of solutions as a delayed evaluation function. The end goal is to produce a contour plot dependent on different parameters in the differential equations. I provide a MWE below.

Take a simple system of differential equations dependent on some parameter b (one could use ParametricNDSolve also) :

test[b_] := NDSolve[{y''[x] - x Sin[x] y[x] == 0, y[0] == 0, y'[0] == b},
y, {x, -10, 10}]


I then define a function which is the solution of this system for a given parameter.

y[b_, x_] := y[x] /. test[b]


Since I have not chosen a parameter b this function must be a delayed evaluation function.

This function gives me no problems when I work with it. However, when I try and perform operations with graphics or symbolic functions with the attribute HoldAll such as Table or ContourPlot I am unable to get Mathematica to evaluate the function and give me the result. For example, with ContourPlot:

ContourPlot[y[b, x], {b, -10, 10}, {x, 0, 10}]


returns

NDSolve::dsvar: 0.0007142857142857143 cannot be used as a variable.
ReplaceAll::reps: {NDSolve[{-5.10204038248*10^-7 y[0.000714285714286]+(y^\[Prime]\[Prime])[0.000714285714286]==0,y[0]==0,(y^\[Prime])[0]==-9.99857142857},y,{0.000714285714286,-10,10}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.


However, even with Evaluate the problem persists:

ContourPlot[Evaluate[y[b, x]], {b, -10, 10}, {x, 0, 10}]

NDSolve::ndinnt: Initial condition b is not a number or a rectangular array of numbers.
ReplaceAll::reps: {NDSolve[{-x Sin[x] y[x]+(y^\[Prime]\[Prime])[x]==0,y[0]==0,(y^\[Prime])[0]==b},y,{x,-10,10}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.


Similar behaviour can be observed with Table

Table[y[i, x], {i, 0, 5}, {x, 0, 10}]
Table[Evaluate[y[i, x]], {i, 0, 5}, {x, 0, 10}]


so I suspect the attribute HoldAll is preventing the evaluation. But somehow, the use of Evaluate and the delayed evaluation also prevents Mathematica from calculating a solution.

I have searched other, similar, questions in Mathematica.SE but I have not found any that reproduce the problem or provide a solution. How do I get Mathematica to work with these sorts of evaluations?

The error messsage actually originates from NDSolve and is realted to x being set to numeric value before the evaluation.

Workaround is to limit the scope of x you use in NDSolve:

test[b_] :=
Block[{x},
NDSolve[{y''[x] - x Sin[x] y[x] == 0, y[0] == 0, y'[0] == b},
y, {x, -10, 10}]]
y[b_, x_] := y[x] /. test[b]
ContourPlot[y[b, x], {b, -10, 10}, {x, 0, 10}]


Second workaround, which is slightly less nice, is to use a different name for a symbol in the plot:

test[b_] :=
NDSolve[{y''[x] - x Sin[x] y[x] == 0, y[0] == 0, y'[0] == b},
y, {x, -10, 10}]
y[b_, x_] := y[x] /. test[b]
ContourPlot[y[b, x2], {b, -10, 10}, {x2, 0, 10}]


which gives the same plot.

Edit

More detailes on what went wrong. From the documentation of ContourPlot:

ContourPlot has attribute HoldAll, and evaluates the $f_i$ and $g_i$ only after assigning specific numerical values to $x$ and $y$.

Let's mimic the effect of HoldAll, by making a function f, which replaces the values as does plot:

In[180]:= f[a_] := a /. {b -> 0, x -> 0}
SetAttributes[f, HoldAll]
f[y[b, x]]

During evaluation of In[180]:= NDSolve::ndinnt: Initial condition b is not a number or a rectangular array of numbers.

During evaluation of In[180]:= ReplaceAll::reps: {NDSolve[{-x Sin[x] y[x]+(y^\[Prime]\[Prime])[x]==0,y[0]==0,(y^\[Prime])[0]==b},y,{x,-10,10}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

During evaluation of In[180]:= NDSolve::dsvar: 0 cannot be used as a variable.

During evaluation of In[180]:= ReplaceAll::reps: {NDSolve[{(y^\[Prime]\[Prime])[0]==0,y[0]==0,(y^\[Prime])[0]==0},y,{0,-10,10}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

Out[182]=
y[0] /. NDSolve[{(y^\[Prime]\[Prime])[0] == 0, y[0] == 0,
Derivative[1][y][0] == 0}, y, {0, -10, 10}]


The probelm is that Pattern f[a_] is already solved before the numeric values are inserted and then they get inserted too agressvily.

ParametricNDSolve takes care of localizing the symbolic x for solving for you and for your application, it is actually the nicesest solution.

• Why exactly do we have this behaviour? Surely, the choice of the symbols for the integration and the plotting shouldn't matter? – OldTomMorris Sep 11 '18 at 9:24
• See the updated answer. – Johu Sep 11 '18 at 13:17

For your example, using ParametricNDSolveValue directly would be simplest:

yp = ParametricNDSolveValue[
{y''[x]-x Sin[x] y[x]==0, y[0]==0, y'[0]==b},
y,
{x, -10, 10},
b
];

ContourPlot[yp[b][x], {b, -10, 10}, {x, 0, 10}]


(The OP asks how ParametricNDSolveValue avoids the issues in using SetDelayed and NDSolve)

The key difference is that ParametricNDSolveValue encapsulates the process of solving the ODE for different values of the parameter, there is no need to worry about evaluation order or conflicting arguments. When a symbolic parameter is given, the ParametricFunction doesn't evaluate:

yp[b]


When a numeric argument is given, the ParametricFunction evaluates to an InterpolatingFunction with no dependence on the ODE independent variable:

yp[1]
FreeQ[yp[1], x]


True

Here are some variations on using a ParametricFunction:

yp[b][x]
yp[b][1]
yp[1][x]
yp[1][1]


1.04674

Unlike the OP approach, no errors occur for any of the above inputs.

• Thanks for the solution. Why does ParametricNDSolve work and how does it avoid the problems @Johu mentions in his answer? I didn't think the difference would be so large between a delayed evaluatin NDSolve and ParametricNDSolve`. – OldTomMorris Sep 11 '18 at 9:26