Kinematics while loop is infinite, plus other errors.

2020 3 19
2 个回答
更新时间：2020 3 21
3 次查看（30 天）

0 个投票

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.

4 个评论
显示 2更早的评论 隐藏 2更早的评论

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

请先登录，再进行评论。

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

回答（2 个）

Cris LaPierre
Cris LaPierre 2020-3-19
编辑：Cris LaPierre 2020-3-19

0 个投票

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
darova 2020-3-19

0 个投票

Here are some tips
You forgot about units: height should be in km for density

10 个评论
显示 8更早的评论 隐藏 8更早的评论

Sam Potter
Sam Potter 2020-3-19
编辑：Walter Roberson 2020-3-20
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
darova 2020-3-19
This should give you an answer
Let me know if it's clear
darova
darova 2020-3-19
Please use special button for code inserting
Sam Potter
Sam Potter 2020-3-20
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
darova 2020-3-20
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
darova 2020-3-20
  • Also does the change in h have to be constant for it to work?
No. All variable can be different (even dt)
Sam Potter
Sam Potter 2020-3-20
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
Sam Potter 2020-3-20
Actually I have just thought,is it because you differentiate displacement to get velocity and then differengiate that to get acceleration?
darova
darova 2020-3-20
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
Sam Potter 2020-3-21
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?

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by