# Kinematics while loop is infinite, plus other errors.

6 views (last 30 days)
Sam Potter on 19 Mar 2020
Commented: Sam Potter on 21 Mar 2020
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.

Show 1 older comment
Sam Potter on 19 Mar 2020
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 on 19 Mar 2020
Looks like differential equation. Can i see it on the paper or picture? Where did you get it?
Sam Potter on 19 Mar 2020
Hi, this is the original sheet of equations (not including suvat).

Cris LaPierre on 19 Mar 2020
Edited: Cris LaPierre on 19 Mar 2020
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)

darova on 19 Mar 2020
Here are some tips
You forgot about units: height should be in km for density

Sam Potter on 20 Mar 2020
Actually I have just thought,is it because you differentiate displacement to get velocity and then differengiate that to get acceleration?
darova on 20 Mar 2020
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 on 21 Mar 2020
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?