# Intermediate variables in NDSolve

I'm interested in using NDSolve as an integrator for a system of differential equations that looks like the following:

x'[t] == f[x[t], p[t]]
p'[t] == g[x[t], p[t]]


Both f and g are expensive functions to compute, but they both need to perform similar computations, so it is cheaper to write

vec'[t] == newf[vec[t]]


where vec is a list whose elements are {x,p}. Now, I have a quantity computed in newf that just depends on x and p at a given time (call it y), and I would like this quantity to be included in the solution. I could compute it separately as

vec'[t] == newf[vec[t]]
y[t] == calcy[vec[t]]


but this would be expensive, as I'd be computing it twice. I'd like to include it in the vector and just include it with newf's output, but as I'm specifying y and not y', this doesn't work. I tried making a list out of everything as

{vec'[t], y[t]} == newf2[vec[t]]


but Mathematica spits back

NDSolve::underdet: There are more dependent variables, {vec[t], y[t]}, than equations, so the system is underdetermined.

Any suggestions?

Here is a minimal working example, where I would like to have y return in the solution as well as the vector solution.

f[vec_List] := Module[{y},
y = Norm[vec];
{1, y}]

NDSolve[{vec'[t] == f[vec[t]], vec[0] == {0, 0}},
{vec},
{t, 0, 1}]

Plot[vec[t] /. %, {t, 0, 1}]

• I would first try to see how it works when nexf2 is linear Commented Feb 22, 2013 at 16:34
• I tried it with {vec'[t],y[t]}=={{0,0},0}, which is as simple as it gets, and received the same error message. Commented Feb 22, 2013 at 16:35
• {vec'[t],y[t]}=={{0,0},0} is not a correct syntax. The correct syntax, with initial condition, is NDSolve[{vec'[t] == {0, 0} , y[t] == 0, vec[0] == {0, 0}}, {vec, y}, {t, 0, 1}] Commented Feb 22, 2013 at 16:54
• Sure, I know this, but it doesn't help with the my question - how to specify the driving term vec'[t] and the auxiliary function y[t] in the same function. Commented Feb 22, 2013 at 16:56
• I don't understand very well the motivation behind making the same function. I think this is not possible, or if it is possible, it doesn't bring anything interessing. You say "computing newf[] and calcy []" would be "twice computing". Do you mean that newf and caly are similar function ? Commented Feb 22, 2013 at 17:18

Very specific to the example in the question, you can do :

NDSolve[{vec'[t] == {1, y[t]}, y[t] == Norm[vec[t]],
vec[0] == {0, 0}}, {vec, y}, {t, 0, 1}]


Here is another solution that corresponds better to what you describe in your comments:

f01[vec_List] := Module[{y}, y = Norm[Take[vec, 2]]; {1, y^2, y}];

NDSolve[{
vec'[t] == driving[t],
driving[t] == f01[vec[t]],
vec[0] == {0, 0, 0},
driving[0] == {1, 0, 0}
},
{vec, driving}, {t, 0, 1}]


encountered difficulties :

• NDSolve understand that driving is of length 3 only after addition of driving[0] == {1, 0, 0}
(which is accepted without error !)

• one can not modify the length of vec and driving in vec'[t] == driving[t]
my solution (efficient ?) was to add a dummy element in vec so that Length[vec]=3

• on the contrary one can extract a part of vec in the function f01 (see y = Norm[Take[vec, 2]])

Note :
With this code NDSolve gives a answer without complaining.
I have not verified that this answer is OK.

Hope this helps

• The minimal working example was just that - a minimal example. What I'm actually wanting to do is to use a long and complicated module, where y is an intermediate result used in calculating vec'[t]. I'm trying to return vec'[t] as well as y. For the trivial example that I gave above, you can do what you suggest, but for the more complicated scenarios, it's not so simple. Although, I guess I could have all of the driving terms stored along with the solution in a similar manner... Commented Feb 23, 2013 at 2:23
• I tried generalizing your suggestion to the following: f[vec_List] := Module[{y}, y = Norm[vec]; {1, y^2, y}]; NDSolve[{vec'[t] == Drop[driving[t], -1], driving[t] == f[vec[t]], vec[0] == {0, 0}}, {vec, driving}, {t, 0, 1}]. However, this gives me an error of a different sort... Commented Feb 23, 2013 at 2:41
• Ok, this is much closer to what I want. I found that as you say, you cannot use Drop[driving[t],-1] in the equation in NDSolve. What you can do, however, is to call a function droplast[driving[t]], where droplast[x_List]:=Drop[x,-1]. I'm really not sure what the difference is, but it seems to work. The importance of this is that the integrator step size will then not be influenced by the dummy function (and it doesn't spend effort integrating the dummy function). Commented Feb 23, 2013 at 20:14
• I have tried something of the kind of your droplast, but it didn't seem to work. Can you post the code in a comment, please ? (I will see that later) Commented Feb 23, 2013 at 20:35
• f01[vec_List] := Module[{y}, y = Norm[vec]; {1, y^2, y}]; g[vec_List] := Drop[vec, -1]; NDSolve[{vec'[t] == g[driving[t]], driving[t] == f01[vec[t]], vec[0] == {0, 0}, driving[0] == {1, 0, 0}}, {vec, driving}, {t, 0, 1}] was what I was using. Perhaps the _List condition is important? Commented Feb 24, 2013 at 1:18