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 trying to workout how to solve a simple 2nd order ODE for y and y'. I can solve for y[t] as follows but how can i solve for y[t] and y'[t] ?

m = 1; g = 9.82;
sol = Flatten @ 
   NDSolve[{y''[t] == -m g - 0.3 y'[t], y'[0] == 10, y[0] == 0, 
     WhenEvent[y[t] < 0, tmax = t; "StopIntegration"]}, 
    y[t], {t, 0, 3}];
Plot[y[t] /. sol, {t, 0, tmax}]

Also the output from NDSolve is an interpolating function but how do i turn this into a function ? i.e i thought it would be something like this but that doesn't seem to work ?

g[t_]:= y[t]/. sol 

Thanks David.

share|improve this question
add comment

2 Answers

To solve for $y(t)$ and $y'(t)$ you simply add $y'(t)$ for the list of the dependent functions you want to solve for.

m = 1; g = 9.82;

sol = Flatten@
   NDSolve[{y''[t] == -m g - 0.3 y'[t], y'[0] == 10, y[0] == 0, 
     WhenEvent[y[t] < 0, tmax = t; "StopIntegration"]}, {y[t],y'[t]}, {t, 0, 3}];

Plot[{y[t] /. sol, y'[t] /. sol}, {t, 0, tmax}]

Mathematica graphics

share|improve this answer
add comment

This is a perfect use case for the (new as of version 9) NDSolveValue function, which can return pretty complicated expressions based on the dependent variables in a system of numerically integrated differential equations. The following, for example, will plot both y and y'.

With[{value = 
   NDSolveValue[{y''[t] == -m g - 0.3 y'[t], y'[0] == 10, y[0] == 0, 
    WhenEvent[y[t] < 0, tmax = t; "StopIntegration"]}, {y[t], 
    y'[t]}, {t, 0, 3}]},
 Plot[value, {t, 0, tmax}]]

Of course, Mathematica knows how to take the derivatives of interpolating functions returned by ordinary NDSolve, so the following will also work fine:

With[{solution = 
   NDSolve[{y''[t] == -m g - 0.3 y'[t], y'[0] == 10, y[0] == 0, 
     WhenEvent[y[t] < 0, tmax = t; "StopIntegration"]}, y, {t, 0, 3}]},
 Plot[{y[t], y'[t]} /. solution // Evaluate, {t, 0, tmax}]]

Also, if you want to avoid the semi-gross tmax = t side-effect in your WhenEvent statement, you can use this quasi-documented feature to get the domain of an interpolating function:

With[{if = 
   y /. First@
     NDSolve[{y''[t] == -m g - 0.3 y'[t], y'[0] == 10, y[0] == 0, 
       WhenEvent[y[t] < 0, tmax = t; "StopIntegration"]}, 
      y, {t, 0, 3}]},
 if["Domain"]]

This returns {{0., 1.86386}} on my machine. Complex objects taking strings like "Domain" as arguments as a way of providing more complex "methods" is a pretty common pattern.

share|improve this answer
    
Thanks for such a detailed response. Very neat. I was trying to differentiate the interpolating function but had trouble turning the interpolation function into a function i could use. –  David McHarg Nov 11 '13 at 23:14
add comment

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.