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.

Answer 1

[Matt Fig]

[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];

[Grufff]

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]

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? :)

[Grufff]

Yes, it *is* familiar ;)

I have it know, and see how the diff command is then used to produce the derivative of my velocity values for acceleration data. I was expecting the reduction of order and two equation set up to output the accel data as part of the my first solution attempt though.

Looking at the "new variable, in and out vector" type method, I have an output vector which is suppsed to give me the derivative of the input vector terms, ie y' and v', which I read as velocity and accel. So I'm still a bit puzzled why I don't get that output (this is still my general unfamiliarity with this stuff I'm sure, I need to practice more).

Anyway, thanks again!

[Matt Tearle]

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]

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

Answer 3

[Grufff]

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.