Intermediate variables in NDSolve

Question

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

Answer 1

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