# Using Dot Product with InterpolatingFunctions

I am using NDSolve on a vector function $\mathbf y'(x)=f[\mathbf y(x)]$ with initial condition $\mathbf y_0$, where the dimension of the vector should be user-defined.

Clear[y, y0]
y0 = {1, 2}
sol = NDSolve[{y'[x] + y[x] == 0, y[0] == y0}, {y}, {x, 0, 10}]


I'm attempting to use the resulting InterpolatingFunction as an input into a matrix ODE, $\mathbf z'(x) = g[\mathbf z(x)]$, with initial condition $\mathbf z_0$.

Clear[z,z0]
z0 = {{1., 2.},{3.,4.}}
failing = NDSolve[{z'[x] + foo[y[x] /. sol, {1, 1}] z[x] == 0, z[0] == z0}, {z}, {x, 0, 10}]


Where

foo[vectx_?(VectorQ[#] &), vecty_?(VectorQ[#] &)] := vectx.vecty


The function foo[x,y] is failing to compute the dot product correctly because the dimensions don't match up (Length[y[x] /. sol,] is 1), although VectorQ[y[x] /. sol] yields True. Naively, I expected that the result of NDSolve would be something like

$\mathbf y(x) = (\text{InterpolatingFunction}[y_1(x)], \text{InterpolatingFunction}[y_2(x)],...)$

So that the dot product

$\mathbf y \cdot (a,b,...)=a \, \text{InterpolatingFunction}[y_1(x)] + b \, \text{InterpolatingFunction}[y_2(x)]+...$

Rather, it's something more like {InterpolatingFunction[{{0.,10.}},<>][t]}, which is trivially a vector since it's enclosed by {...}.

How do I express my "dot product" logic above in Mathematica-ese? Or is there a more elegant way to get the code flow that I want?

• You can change sol = NDSolve... to sol = First@NDSolve.. to get y[x]/.sol to have the form InterpolatingFunction[...] instead of {InterpolatingFunction[...]}
– ssch
Jun 24, 2013 at 10:33
• The (VectorQ[#] &) is not needed, VectorQ will do the job in this case.
– mmal
Jun 24, 2013 at 12:20
• Alternatively (VectorQ[#, NumericQ]&) is more precise. (Note that VectorQ[y[x] /. sol] is True, but VectorQ[y[0.4] /. sol] is False.) One could also add a condition, foo[..] /; Length[vectx] == Length[vecty] := .., but that's not really necessary. The major problem has been pointed out by ssch and partial81. Jun 24, 2013 at 14:11
• @ssch, this worked perfectly, thank you. Jun 24, 2013 at 21:41

I think the following could help:

failing = NDSolve[{z'[x] +
foo[Evaluate[y[x] /. sol][[1]], {1, 1}] z[x] == 0, z[0] == z0}, {z}, {x, 0, 10}]


Please tell me if this works for you.

Btw. You can also skip the Evaluate and use only foo[y[x] /. sol[[1]], {1, 1}]. The result is the same.

• This works, cheers! Jun 24, 2013 at 21:43
• Great to hear! Btw. First[expr] (as ssch writes) and expr[[1]], are equivalent. So, chose the notation you like more ;-) Jun 25, 2013 at 6:19

You asked if there's a more elegant way to achieve your output. I think the simplest approach is to combine your 2 ODEs into a single ODE, something like:

y0 = {1, 2};
z0 = {{1, 2}, {3, 4}};

{ysol, zsol} = NDSolveValue[
{
y'[x] + y[x] == 0, y[0] == y0,
z'[x] + y[x].{1,1} z[x] == 0, z[0] == z0
},
{y, z},
{x, 0, 10}
];


Note that I used NDSolveValue instead of NDSolve, and just used Dot instead of your auxiliary function foo.

Visualization:

Plot[zsol[x], {x, 0, 10}]