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.

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?

share|improve this question
1  
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 '13 at 10:33
    
The (VectorQ[#] &) is not needed, VectorQ will do the job in this case. –  mmal Jun 24 '13 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. –  Michael E2 Jun 24 '13 at 14:11
    
@ssch, this worked perfectly, thank you. –  cosmoguy Jun 24 '13 at 21:41

1 Answer 1

up vote 1 down vote accepted

I think the following could help:

Use in your code

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.

share|improve this answer
    
This works, cheers! –  cosmoguy Jun 24 '13 at 21:43
    
Great to hear! Btw. First[expr] (as ssch writes) and expr[[1]], are equivalent. So, chose the notation you like more ;-) –  partial81 Jun 25 '13 at 6:19

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.