How to get value from InterpolatingFunction

Question

I have numerically solved a system of ODEs using NDSolve. I can plot them using ParametricPlot but I am having trouble evaluating them at specific points and getting the derivative. For example, if z0[s] is the solution to one of my equations, I want z0[2] or z0[3].

In the help file for InterpolatingFunction, their example seems trivial. They find the solution ifun of an ODE using NDSolve and then use the command ifun[2.5] to evaluate the function at 2.5. This is the functionality I want. I am having trouble seeing why this isn't working for my situation. Here is my code:

eq={(z0'')[s]==0,(z1'')[s]==0,(z2'')[s]==0,(z3'')[s]==0};
initCond={z1[0]==0,z0[0]==0,z2[0]==2,z3[0]==0,z0'[0]==1,z1'[0]==-1,z2'[0]==0,z3'[0]==0};
fullEq=Union[geoEq1,initCond1];
sol=NDSolve[fullEq,{z0[s],z1[s],z2[s],z3[s]},{s,0,20}]

Note that the systems of ODEs I use is very basic (the solution is just a straight line). I just want to learn how to do this for for a general function first.

This output

{{z0[s] -> InterpolatingFunction[{{0., 20.}}, <>][s], 
  z1[s] -> InterpolatingFunction[{{0., 20.}}, <>][s], 
  z2[s] -> InterpolatingFunction[{{0., 20.}}, <>][s], 
  z3[s] -> InterpolatingFunction[{{0., 20.}}, <>][s]}}

According to the help file, I think I should be able to do z0[2] to get the value and z0'[2] to get the value of the derivative. I try z0[2] and z0[2] /. sol but this doesn't yield the result.

For completeness, I have checked related questions on this SE but most seem to be for either a single InterpolatingFunction (i.e. not a system of 4 functions) and not for a parametrized solution z0[s]. I have not yet found any solution on the help file. I will keep looking and post back if I solve it.

Answer 1

There are basically, say, 2 x 2 = four forms of results one can extract from NDSolve, aside from some special-use cases. I think one of these four tends to be more useful than the others, although occasionally I find the other forms useful in particular cases. The one I like is based on a central idiom of Mathematica, the replacement Rule, and works well in Mathematica. The idea is not common in other languages, and so some may tend to avoid it. The one I like is also in terms of functions $f$ instead of an expression $f(x)$. Despite that everything is an expression in Mathematica, since we're doing calculus here, dealing in terms of functions plays to one of Mathematica's strengths.

The four forms can be viewed as the combinations of two binary divisions. The first division is the choice of the result as a rule or an expression:

NDSolve[eqns, u, {x, xmin, xmax}]        (* replacement rule *)
NDSolveValue[eqns, u, {x, xmin, xmax}]   (* function expression *)
  (* equiv. to u /. NDSolve[eqns, u, {x, xmin, xmax}] *)

The second is the form of u, either a function or the function-value at the dependent variable x; for instance

NDSolve[eqns, y, {x, xmin, xmax}]        (* (interpolating) function *)
NDSolve[eqns, y[x], {x, xmin, xmax}]     (* function-value at x *)

The one I like is

NDSolve[eqns, y, {x, xmin, xmax}]

which returns a rule of the form y -> InterpolatingFunction[<>] (inside lists, of course). This can be used to inject a function (the solution) into any expression. Mathematica will compute the rest, including any derivatives and integrals of y that appear in the expression. This is the principal reason I like it.

By contrast

NDSolve[eqns, y[x], {x, xmin, xmax}]

returns a rule of the form y[x] -> InterpolatingFunction[<>][x]. If it is applied to an expression, it will only replace y[x], literally. It will not replace y[2.5] or y'[x], and so forth. This is where the OP ran into trouble with his or her code.

Example

I'll use a slight modification of the OP's example. I'll reduce the number of functions to {z0, z1} and use slightly more interesting DEs. Often I find it is convenient to collect all the functions in the differential equation in a list such as vars below. If vars = {z0, z1}, then Through[vars[s]] equals {z0[s], z1[s]}. The solution solIF is in the form I prefer, {z0 -> InterpolationFunction[<>],...} and solS is in the form {z0 -> InterpolationFunction[<>][s],...}.

In another side note, the assignment {solIF} =... works if eqs are quasilinear system, because the solution will consist of a single set of functions for the dependent variables. If there is an error and NDSolve does not return a solution, the assignment will fail. One advantage over using solIF = NDSolve[...] is in code that uses solIF immediately: On the one hand, if NDSolve fails, solIF will equal the original NDSolve[] call, and every time solIF is evaluated, it leads to another failed evaluation of NDSolve[]. On the other hand, if the system is not quasilinear (and I expected it to be so) and NDSolve returns more than one set of solutions, then the assignment will fail. This lets me knowing I wasn't thinking straight or that the system was not constructed properly. (It's a minor point, but I thought I should explain why I do it.)

(* damped/harmonic oscillators *)
eqs = {z0''[s] + z0'[s]/8 + z0[s] == 0, z1''[s] + 4 z1[s] == 0};
ics = {z0[0] == 0, z0'[0] == 1, z1[0] == 0, z1'[0] == -1};
vars = {z0, z1};
{solIF} = NDSolve[{eqs, ics}, vars, {s, 0, 20}]             (* z0 -> ifn, etc *)
{solS} = NDSolve[{eqs, ics}, Through[vars[s]], {s, 0, 20}]  (* z0 -> ifn[s], etc *)

To evaluate z0 at 2.5, say, here's how to use each form. It is somewhat more straightforward with the z0 -> InterpolatingFunction[<>] form in solIF.

z0[2.5] /. solIF
z0[s] /. solS /. s -> 2.5

Computing the derivative z1'[s] is even more natural with solIF. In the function-value expression form (like the OP's), only the last form of the derivative of z1 is computed correctly, since it is the only code in which the replacement rules are applied to z1[s].

z1'[s] /. solS
D[z1[s], s] /. solS
D[z1[s] /. solS, s]
(*
  Derivative[1][z1][s]
  Derivative[1][z1][s]
  InterpolatingFunction[{{0., 20.}}, <>][s]
*)

By contrast, all three form work with the function form:

z1'[s] /. solIF
D[z1[s], s] /. solIF
D[z1[s] /. solIF, s]

Often, we use the solutions to perform other calculations. For instance, we can examine the energy in the oscillator systems represented by {eqs, ics}. Here is some code for computing the energy:

ClearAll[energy];

SetAttributes[energy, Listable];
energy[harmonicODE_, s0_] := Module[{z, s, m, γ, k},
   {{z, s}} = Cases[harmonicODE, Derivative[2][z_][s_] :> {z, s}, Infinity, 1];
   {k, γ, m} = Coefficient[harmonicODE /. Equal -> Subtract, {z[s], z'[s], z''[s]}];
   k z[s0]^2/2 + m z'[s0]^2/2
   ];

energy[eqs, s]
(*
  {z0[s]^2/2 + 1/2 z0'[s]^2, 
   2 z1[s]^2 + 1/2 z1'[s]^2}
*)

These are symbolic expressions for the energies. We can evaluate them by injecting the solution functions in solIF into the energies and plugging in values for s, say, via Plot[]:

Plot[energy[eqs, s] /. solIF // Evaluate, {s, 0, 20}]

I don't know a simple way to do that with solS.

I won't say much about NDSolveValue[]. It returns InterpolatingFunction[<>] or InterpolatingFunction[<>][s]. It can be convenient if you just want to evaluate the function. One can inject them in an expression by writing your own rules, but that can be a bit inconvenient if you have several variable {z0, z1,...} as in the OP's case. You're probably better off using NDSolve[] in that case.

Answer 2

Try this, hopefully I have correctly fixed what I think where typos

geoEq1 = {z0''[s] == 0, z1''[s] == 0, z2''[s] == 0, z3''[s] == 0};
initCond1 = {z1[0] == 0, z0[0] == 0, z2[0] == 2, z3[0] == 0, 
  z0'[0] == 1, z1'[0] == -1, z2'[0] == 0, z3'[0] == 0};
fullEq = Union[geoEq1, initCond1];
{o0, o1, o2, o3} = {z0, z1, z2, z3} /. 
  NDSolve[fullEq, {z0, z1, z2, z3}, {s, 0, 20}][[1]];
Plot[{o0[s], o1[s], o2[s], o3[s]}, {s, 0, 20}]
{o0[1], o1[3], o2[2], o3[4]}

and you should see the four overlaid plots and four evaluations of your interpolating functions.

I used o0,o1,o2,o3 because I wasn't sure if you were going to need to maintain the original values of your zn assignments.

Note the [[1]] trailing the NDSolve to eliminate one extra layer of {}. Note the change in form of your zn inside the argument list to NDSolve. Note the way that each on is initialized and later used.

There are at least a dozen other ways of doing anything in Mathematica. Pick one that you can master and remember for next time.