Running a kinematics while loop containing an if loop. I get the error: Index exceeds number of array elements(1).

h(1)=150000;            %initial height
a(1)=(40*10^7)/(6371+h(1))^2;        %initial acceleration dependant on height
dt=0.005;        %time step
t(1)=0;         %initial time
v(1)=a(1)*t(1);         %velocity
g(1)=((40*10^7)/(6371+h(1))^2);     %Downward force h>100000, upward force = 0 above 100000
As= 5;      %Area
m=850;   %Mass
c=0.7;
p(1)=(100000/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*As*v^2;     %Downward force h<=100000
i=1; %loop counter
while h(end)>=0
    t(i+1)=t(i)+dt;
    h(i+1)=150000-(a(i+1)*(t(i+1))^2);
    if h(i+1)>100000
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=g(i+1);
    else
        p(i+1)=((h(i+1))/71)+1.4;
        Fd(i+1)=0.5*(p(i+1))*c*As*(v(i+1))^2;
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=(g(i+1))-((Fd(i+1))/m);
    end
    v(i+1)=(a(i+1))*(t(i+1));
    i=i+1;
end

I have posted this before and got some very useful advice but it did not fully fix the problem. Any help would be appreciated.

Answer 1

[Sriram Tadavarty]

Hi Sam,

In the second statement of while loop, the variable a is indexed with value of 2, before setting the value of it. This is the source of error.

h(i+1)=150000-(a(i+1)*(t(i+1))^2); % Before this execution, make sure the varaibles a and t, have a value at 'i+1' index

From the code it is not pretty clear as what you intended to do. More details on problem statement would help.

But, based on code comments and some assumptions, i can suggest some modifications to the code in loop:

• Update the height and velocity before the increment time
• Index the variable with one step before the increment time
• Then increment the time

while h(end)>=0
    v(i+1)=(a(i))*(t(i)); % Find the velocity of previous time increment
    h(i+1)=150000-(a(i)*(t(i))^2); % Find the height of previous time increment
    t(i+1)=t(i)+dt; % Update the time step
    if h(i+1)>100000
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=g(i+1);
        p(i+1)=p(1); % Set this value to some default value, chosen as the initial value (Since this is used directly with the indexing of 'i+1' in the else statement)
        Fd(i+1)=Fd(1); % Same reason as p
    else
        p(i+1)=((h(i+1))/71)+1.4;
        Fd(i+1)=0.5*(p(i+1))*c*As*(v(i+1))^2;
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=(g(i+1))-((Fd(i+1))/m);
    end
    i=i+1;
end

The udpate of your code with above code, will remove any errors that you see.

Hope this helps.

Regards,

Sriram

[Sam Potter]

I just saw your comment and came back to this thread. To simplify my code I got ride of the if loop and made it so it is a while loop from 100000<h>0.

>> clear

>> h(1)=100000; %initial height

a(1)=0.03535;

dt=0.005; %time step

t(1)=0; %initial time

v(1)=59.29; %Initial velocity

g(1)=((40*10^7)/(6371+h(1))^2); %Downward force h>100000, upward force = 0 above 100000

As= 5; %Area

m=850; %Mass

c=0.7;

p(1)=(100000/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*As*(v(1)^2); %Downward force h<=100000

i=1; %loop counter

while h(end)>=0

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

h(i+1)=10000-(v(i+1)*t(i+1))-(0.5*a(i+1)*(t(i+1))^2);

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

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

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

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

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

i=i+1;

end

However, I still get the error index exceed array elements. Any help would be apprecieted.

[Sriram Tadavarty]

May i know, if the previous while loop code worked?

If yes, the same can be applied here

[Sam Potter]

HI, the code did run but the plotted graphs of (t,v) and (a,v) didnt look correct.

[Sriram Tadavarty]

Thanks for the update Sam.

The errors comes from the first line of while loop, which index'es second element of variables 'a' and 't'. Like a(i+1) => a(2) which is not defined, similarly t(i+1) => t(2) also not defined. So, the error.

Update 'a' and 't' such that those values are known before. This should solve that error.

[Sam Potter]

Hi, my code already has a value for a(1) and t(1). Does this not mean the values are known before? If not how do i make it so that a and t are known before?

[Sriram Tadavarty]

Sam, your code only knows a(1) and t(1), but in the first line of while loop

v(i+1)=(a(i+1))*(t(i+1));
% Here you are accessing a(2) and t(2), which is not known (Since i is 1 before the while loop)

Hope you can see the issue now. As second element of matrix a is accessed without having that value, MATLAB errors.

So, i advice you to perform the loop from initial value and then update the time step. I could help more, if i know the exact equations.

[Sam Potter]

Hi, my code is currently this:

clear
h(1)=150000;            %initial height
a(1)=(40*10^7)/(6371+h(1))^2;        %initial acceleration dependant on height
dt=0.005;        %time step
t(1)=0;         %initial time
v(1)=a(1)*t(1);         %velocity
g(1)=((40*10^7)/(6371+h(1))^2);     %Downward force h>100000, upward force = 0 above 100000
As= 5;      %Area
m=850;   %Mass
c=0.7;
p(1)=(h(1)/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*As*v^2;     %Downward force h<=100000
i=1; %loop counter
while h(end)>=0
    t(i+1)=t(i)+dt;
    h(i+1)=150000-(a(i)*(t(i+1))^2); % Find the height of previous time increment
    if h(i+1)>100000
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=g(i+1);
        p(i+1)=p(1); % Set this value to some default value, chosen as the initial value (Since this is used directly with the indexing of 'i+1' in the else statement)
        Fd(i+1)=Fd(1); % Same reason as p
        v(i+1)=(a(i+1))*(t(i+1)); 
    else
        p(i+1)=((h(i+1))/71)+1.4;
        Fd(i+1)=0.5*(p(i+1))*c*As*(v(i+1))^2;
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=(g(i+1))-((Fd(i+1))/m);
        v(i+1)=59.29+(a(i+1))*(t(i+1)); 
    end
    i=i+1;
end

It works find until the error comes up: 'Index exceeds the number of array elements (237868).'. The number shown is when the loop switches from the if to the when loop. This means the error is with the while part of the loop (I think).

[Walter Roberson]

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

You do not assign to v(i+1) until 3 lines further down.

[Sriram Tadavarty]

Yes. If the update to v is placed ahead of update to Fd, the issue will go away.

% In the else statement
% you can try placing the update of v ahead with previous acceleration value
        p(i+1)=((h(i+1))/71)+1.4;
        v(i+1)=59.29+(a(i))*(t(i+1)); 
        Fd(i+1)=0.5*(p(i+1))*c*As*(v(i+1))^2;
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=(g(i+1))-((Fd(i+1))/m);
% Or you can just update the Fd value with previous v
Fd(i+1)=0.5*(p(i+1))*c*As*(v(i))^2;

Either of the two should work, but it depends on what is required there.

[Sam Potter]

clear
h(1)=150000;            %initial height
a(1)=(40*10^7)/(6371+h(1))^2;        %initial acceleration dependant on height
dt=0.005;        %time step
t(1)=0;         %initial time
v(1)=a(1)*t(1);         %velocity
g(1)=((40*10^7)/(6371+h(1))^2);     %Downward force h>100000, upward force = 0 above 100000
As= 5;      %Area
m=850;   %Mass
c=0.7;
p(1)=(h(1)/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*As*v^2;     %Downward force h<=100000
i=1; %loop counter
while h(end)>=0
    t(i+1)=t(i)+dt;
    h(i+1)=150000-(0.5*a(i)*(t(i+1))^2); % Find the height of previous time increment
    if h(i+1)>100000
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=g(i+1);
        p(i+1)=p(1)
        Fd(i+1)=Fd(1); % Same reason as p
        v(i+1)=(a(i+1))*(t(i+1)); 
    else
        v(i+1)=59.29+(a(i))*(t(i+1)); 
        p(i+1)=((h(i+1))/71)+1.4;
        Fd(i+1)=0.5*(p(i+1))*c*As*(v(i+1))^2;
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=(g(i+1))-((Fd(i+1))/m); 
    end
    i=i+1;
end

Hi, I have made he chnages you said, by placing v(i+1) before Fd. The code will run, but the graphs it creates of (t,v) and (t,a) look very wrong. Do you think this is due to the loop or a mistake in the equations?

[Sriram Tadavarty]

Hi Sam, It seems like the way equations are formulated to code, went wrong then. May I know what are the expected graphs for (t,v) and (t,a)? It would be helpful to easily find out what went wrong by writing the equations for the first three time steps in notes and see if the code is behaving as such, as a first step. Regards Sriram

[Sam Potter]

I have linked the code I used for the graphs of (t,v) an (t,a) while h>100000. The graphs should continue on smoothly from while h>0, with maybe a less steep graph.

h(1)=150000;    %Initial height
g(1)=(40*10^7)/(6371+h(1))^2;
dt=0.005;
t(1)=0;
v(1)=g(1)*t(1);

i=1;
while h(end)>=100000
    t(i+1)=t(i)+dt;
    h(i+1)=150000-(g(i)*(t(i+1))^2;
    g(i+1)=(40*10^7)/(6371+h(i+1))^2);
    v(i+1)=g(i+1)*t(i+1);
    i(i+1);
end

(I tried to lnk the graphs but an error arose. If you plot (t,v) and (t,g) for this code you will se the type of graphs I am expecting). It is worth noting that on the (t,v) graph I get for the while and if loop, the graph looks ok until h hits 100000.

[Sam Potter]

Just to add on, I have realised that when the code switches to the else section, the height dramatically increases, when it should keep on steadily decreasing. Also, when the code switches to the else section, the accelertion changes to a very large negative answer.

[Sriram Tadavarty]

Yes Sam. Even, i did observe it. But, I am not sure what your intended code should do.

[Sam Potter]

My code should give me the values of v and a when the object hits the ground.

Initial v=0, and mass = 850 Initial H=150000 Initial A=5

Whilst h>100000 there is no air resistance, meaning g, the acceleration = (40*10^7)/((6371+h)^2).

Whilst h>0, the acceleration is equal to g-(Fd/m). Fd=0.5*0.7*((H/71)+1.4)*A*(v^2).

I used the suvat equation for height, 150000-ut-(0.5*a*t^2). and velocity (v=u+at).

Hopefuly this might give you some context on the equations I have used in my code.

Do you think it would be easier to do a seperate while loop for h<10000?

[Sriram Tadavarty]

Thanks for the explanation of what is intended. Yes Sam, a new while loop would benefit. I see that your update to velocity is not updating properly, since u is ignored. For only first time step u is 0, for all other time steps, it must be the previous value. Right?

[Sam Potter]

No, (unless I am misunderstanding you), I believe u is constantly 0, whch is why I ignored it. One thing I forgot to mention is that I used the value of 59.29 in my else loop as that is the last value of v whilst h>100000.

This is the seperte while loop I have made:

>> clear

h(1)=100000;

t(1)=0;

dt=0.005;

u=59.29;

a(1)=0.03535;

v(1)=u+(a(1)*t(1)); %Initial velocity

p(1)=(((h(1))/71)+1.4); %Air Density

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

A=5;

c=0.7;

m=850;

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

i=1;

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);

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);

v(i+1)=u+(a(i+1)*t(i+1));

i=i+1;

end

The issue with this loop is that it seems to be infinite.