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'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}]
share|improve this question
    
I would first try to see how it works when nexf2 is linear –  andre Feb 22 '13 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. –  Jolyon Feb 22 '13 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}] –  andre Feb 22 '13 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. –  Jolyon Feb 22 '13 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 ? –  andre Feb 22 '13 at 17:18
show 1 more comment

1 Answer

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

share|improve this answer
    
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... –  Jolyon Feb 23 '13 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... –  Jolyon Feb 23 '13 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). –  Jolyon Feb 23 '13 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) –  andre Feb 23 '13 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? –  Jolyon Feb 24 '13 at 1:18
show 6 more comments

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.