ode45 - solving 2nd Order ODE IVP problem

Question

Grufff 2011-5-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/7332-ode45-solving-2nd-order-ode-ivp-problem

I hope someone can help, I thought I understood this but I can't make it work and it's driving me bananas.

I've been given a simple 2nd order ODE which is supposed to describe the motion of a mars lander type vehicle in the last 6 seconds of its descent to a planet surface - this period is where the retro rockets are firing and hence slowing the vehicle.

I'm told that the rockets are fired 20m from the surface, when the velocity of the lander is 67.056m/s (and I take this to be at time zero, I need to solve for time span 0 - 6secs).

Equation: y'' = g - (k/m)*(y')^2

g=3.885m/s^2

k=1.2

m=150kg

(these are g, gravity of the planet, k, air resistance co-eff, and m, mass of the lander)

So, it's a 2nd order equation, before passing it to ode45 I need to reduce it to a system of 2 1st order ODEs and then call my function giving time span and initial values. This is my attempt (apologies if my variable names etc don't follow usual convention, I'm new to this):

y'' = g - (k/m)*y'^2

introduce new variable:

v = y'

v' = y''

input and output vectors:

input = [y;v]

output =[y';v']

(I mean y prime and v prime here, not matlab transpose operator of course).

Giving my system:

output(1) = input(2) output(2) = g - (k/m)*(input(2)^2)

So, I'm writing my function as:

********

function output = lander (t, input)
g = 3.688;
output(1) = input(2);
output(2) = g - (1.2/150)*(input(2)*input(2));
output = output(:);

**********

and calling with command:

[t,output] = ode45(@lander, [0:0.05:6], [20,67.056]);

I'm not getting the results I'm supposed to, and can't figure this out.

Behaviour should be the velocity reduces from 67m/s down to close to zero, and acceleration reduces from -30 (ish?) toward zero. COrrect result should be this:

My velocity seems correct but acceleration is keeps climbing. It looks like I'm messing up my conversion to 1st order, or using my IVP's incorrectly. Mine result here:

I'be ever so grateful for a pointer or two in the right direction. Please point out my errors!

The solution I've been given, which produces the data trying to match, uses a different method:

basically uses a function containing one equation

DV = G - K/M * V.^2;

it's used once to calculate velocity against time, then run through the same equation a second time using the v data to produce accel data.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt Fig 2011-5-13

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/7332-ode45-solving-2nd-order-ode-ivp-problem#answer_10077

在 MATLAB Online 中打开

[t,y] = ode45(@lander, 0:.005:6, [20,67.056]);
% Calculate the acceleration from the velocity...
acc = diff(y(:,2))/(t(2)-t(1)); % Accel. is der. of velocity!
plot(t,y(:,2),'r',t(2:end),acc,'b')
legend({'Velocity [m/s]';'Acceleration [m^2/s]'})
xlabel('Time [s]')

I used this LANDER:

function ydot = lander (t,y)
g = 3.885;  % Force of gravity
k = 1.2;  % Air resistance
m = 150;  % lander mass
ydot = [y(2);g - (k/m)*y(2).^2];

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Grufff 2011-5-14

Yes, your function is much neater than mine. I was pretty sure I could do the element by element ^2 squaring using y(2).^2 but I was trying to reduce my function down to it's simplest terms to cut out any errors which might be messing up my results.

I like your last line, it didn't occur to me to combine both of my equations in one line in that way, very neat.

Thanks again for your help.

请先登录，再进行评论。

Answer 2

Matt Tearle 2011-5-13

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/7332-ode45-solving-2nd-order-ode-ivp-problem#answer_10078

在 MATLAB Online 中打开

Nothing's wrong, except your interpretation of the output. The columns of output correspond to y and y' (v), not v and v'. To check:

plot(t(2:end),diff(output)./[diff(t),diff(t)])

This plots an approximation to the derivatives of the columns of output. Look familiar? :)

3 个评论
显示 1更早的评论隐藏 1更早的评论

Matt Tearle 2011-5-16

Your rate equation function (lander) takes y & v as inputs and returns y' & v' as outputs. However, the ODE solver (ode45) takes these equations and initial conditions as inputs, then, after integrating, returns the solution y & v.*

So another way to plot the derivatives would be to stick y & v (the output from ode45) into lander, and plot the output. This will be messy, though, because of the way the dimensions work. You'd probably have to use a loop.

*

[If it helps, don't think in terms of y & v, but of a vector z. Your equations are z' = f(z). lander implements the function f. ode45 integrates these equations to get z(t).]

Grufff 2011-5-19

Thanks for the further explanation, it's very helpful and I appreciate your taking the time.

请先登录，再进行评论。

Answer 3

Grufff 2011-5-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/7332-ode45-solving-2nd-order-ode-ivp-problem#answer_10079

Thank you both so much!

I'm ever so glad to know I wasn't getting it badly wrong, just some mis-understanding of my output. I get it now you point out my mistake.

For some reason I thought my function was outputting v and v', but as you say , it's y and v. Thank you, both.

I don't think I can accept two answers, so I will accept Matt Fig's, only because he was first to reply. You were both super speedy and very helpful, thank you again.