Kinematics while loop is infinite, plus other errors.

I have posted similar things before, so apologies if you have already seen this.

My code is this:

clear
h(1)=100000;   %Initial Height
t(1)=0;
dt=0.005;  %Time Step
u=59.29;  %Initial Velocity
a(1)=0.03535;  %Initial Acceleration
v(1)=u+(a(1)*t(1));  %Velocity
p(1)=(((h(1))/71)+1.4);   %Air Density
g(1)=(40*10^7)/((6371+h(1))^2);   %Gravity
A=5;  %Area
c=0.7;
m=850;  %Mass
Fd(1)=0.5*c*(p(1))*A*(v(1))^2;   %Air Resistance
i=1;   %Loop counter
while h(end)>=0
    t(i+1)=t(i)+dt;
    h(i+1)=100000-(u*t(i+1))-(0.5*a(i)*t(i+1)^2);  %Suvat s=ut+0.5*a*t^2
    p(i+1)=(((h(i+1))/71)+1.4);
    g(i+1)=(40*10^7)/((6371+h(i+1))^2);
    Fd(i+1)=0.5*c*(p(i+1))*A*(v(i))^2;
    a(i+1)=g(i+1)-(Fd(i+1)/m);   %Acceleration=Gravity-(Fd/m)
    v(i+1)=u+(a(i)*t(i+1));   %Suvat v=u+at
    i=i+1;
end

The code is infinite. When I stop it running and plot a graph of (t,h) it appears that h drastically increases after the first time step. I cannot see why this is the case as the equation for height should mean h decreases as time goes in. Pehaps it is the order I have had to code my variables? I have tried various orders but coudnt get any to work without the error 'Index exceeds array elements (1)'. Any help would be appricieted.

[darova]

Can you show original equations?

[Sam Potter]

Hi. The initial height (h) of this part of the question is 100000. I have worked out the initial accelertion (a) at this height to be 0.03535.The initial velocity (u) is 59.29

Initial time (t) =0

The equations whilst 0>h>100000,. Mass (m) m=850.

Air constant (c )=0.7

Area = 5

Height (h) = 100000-(u*t)-(0.5*a*(t^2))

Gravity (g) = (40*10^7)/((6371+h)^2)

Air density (p)= (h/71)+1.4

Air resistance (Fd) =(0.5*c*p*A*(v^2)).

acceleration (a) = g-(Fd/m)

Velocity (v) = u+at

I think that is everything. If necessary I can supply the code of 150000>h>100000 if you want to see how I got my results for the first section of the flight. Thank You.

[darova]

Looks like differential equation. Can i see it on the paper or picture? Where did you get it?

[Sam Potter]

Hi, this is the original sheet of equations (not including suvat).

Answer 1

[Cris LaPierre]

Quick observation is that you are not being careful with your units. H is specified as height in km, but it looks like you are converting it to meters.

It also looks like you missed the negative sign in the calculation of rho (-H/71 + 1.4)

Answer 2

[darova]

Here are some tips

You forgot about units: height should be in km for density

[Sam Potter]

Hi, I have made the changes you said, my code now looks like this and runs well:

clear
h(1)=100000;   %Initial Height
t(1)=0;
dt=5;  %Time Step
u=59.29;  %Initial Velocity
a(1)=0.03535;  %Initial Acceleration
v(1)=u+(a(1)*t(1));  %Velocity
p(1)=-((h(1)/1000)/71)+1.4;   %Air Density
g(1)=(40*10^7)/((6371+h(1))^2);   %Gravity
A=5;  %Area
c=0.7;
m=850;  %Mass
Fd(1)=0.5*c*(p(1))*A*(v(1))^2;   %Air Resistance
i=1;   %Loop counter
while h(end)>=0
    t(i+1)=t(i)+dt;
    h(i+1)= h(i)-(v(i)*dt)  %Suvat s=ut+0.5*a*t^2
    p(i+1)=-((h(i+1)/1000)/71)+1.4; 
    g(i+1)=(40*10^7)/((6371+h(i+1))^2);
    Fd(i+1)=0.5*c*(p(i+1))*A*(v(i))^2;
    a(i+1)=g(i+1)-(Fd(i+1)/m);   %Acceleration=Gravity-(Fd/m)
    v(i+1)=v(i)+(a(i+1)*dt);   %Suvat v=u+at
    i=i+1;
end

I understand the changes to air density and the time step, but please can you explain the chanages to h and v in the loop? Especially using dt instead of t(i+1), as I dont understand how using a constant value for time can work the way it does. Sorry if its a stupid question. Thank You.

[darova]

This should give you an answer

Let me know if it's clear

[darova]

Please use special button for code inserting

[Sam Potter]

Hi,thank you for your help. Just wondering,is there a general formula that shows your point? Also does the change in h have to be constant for it to work?

[darova]

It's not that easy to explain it on your case. You have complicated example

Your equation is:

can be re-written as:

and g are depend on v and H which depend on time

I'll give you some links

Euler method - simplest method to solve numericaly diff. equation (used in your case)

ODE - what is differential equation

ode45 - main MATLAB solver of ordinary diff equations

Feel free to ask

[darova]

• Also does the change in h have to be constant for it to work?

No. All variable can be different (even dt)

[Sam Potter]

I think I understand rulers method,the one thing I am not sure on is why we use v(i) for the h(i+1) equation. Is it because for each time step v*t is the distance travelled?

[Sam Potter]

Actually I have just thought,is it because you differentiate displacement to get velocity and then differengiate that to get acceleration?

[darova]

You can use v*t to calculate distance only when you have v=constant (no changing)

But velocity in this case changes all the time

• Actually I have just thought,is it because you differentiate displacement to get velocity and then differengiate that to get acceleration?

There is no differentiation here. Only integration: h(i+1) = h(i) + v(i)*dt - integration (summation)

I can't explain it to you here. That is not that simple. You should dig into this by yourself. Practice

[Sam Potter]

Hi, whilst this issue is not the same one as the thread title, it is the same block of code. Im the question, when the object hits h=3000, a 150m area parachute opens. I have tried tos how this in my code:

h(1)=150000; %initial height

a(1)=(40*10^7)/(6371+h(1))^2; %initial acceleration dependant on height

dt=5; %time step

t(1)=0; %initial time

v(1)=a(1)*t(1); %velocity

g(1)=((40*10^7)/(6371+h(1))^2);

A= 5; %Area of spaceship

m=850; %Mass

c=0.7;

p(1)=-(100/71)+1.4; % Initial Air Density (Air density occurs at h=100000, from then p=(h/71)+1.4)

Fd(1)=0.5*p(1)*c*A*v^2; %Downward force h<=100000

i=1; %loop counter

while h(end)>=0

t(i+1)=t(i)+dt;

h(i+1)=h(i)-(v(i)*dt); % Find the height of previous time increment

g(i+1)=(40*10^7)/((6371+h(i+1))^2);

if h(i+1)<=3000

A=5;

else

A=150;

end

if h(i+1)>100000

a(i+1)= g(i+1) %Acceleration=Gravity-(Fd/m)

v(i+1)=v(i)+(a(i+1)*dt);

elseif h(i+1)>=0

p(i+1)=-((h(i+1)/1000)/71)+1.4;

Fd(i+1)=0.5*c*(p(i+1))*A*(v(i))^2;

a(i+1)=g(i+1)-(Fd(i+1)/m); %Acceleration=Gravity-(Fd/m)

v(i+1)=v(i)+(a(i+1)*dt);

end

i=i+1;

end

I showed this with the lines

if h(i+1)<=3000

A=5;

else

A=150;

end

However, when i run it, what happens it acceleration dramitically decreases until it is a huge negative, also causing velocity to be negative.

What should happen is a drop in acceleration, and a drop in velocity, but it should then level out. Any help oon this would be appreciated. Is the issue to do with where on my code I have inserted the new if loop?