5
$\begingroup$

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?

$\endgroup$
4
$\begingroup$

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

enter image description here

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.

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

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

enter image description here

Addendum

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

enter image description here

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]

enter image description here

True

Here are some variations on using a ParametricFunction:

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

enter image description here

1.04674

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – OldTomMorris Sep 11 '18 at 9:26

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.