Code not coming out of if statement.

clear all , clc, format compact
L = 30;
m = 68.1;
c_d = 0.25;
k = 40;
gamma = 8;
g = 9.81;
to = 0;
tf = 50;
h = 0.1;
n = ((tf-to)/h);
t = to:h:tf;
x1 = 0;
y1 = 0;
x=x1;
y=y1;

f1 = @(y) g - sign(y).*cd/m.*y.^2;
f2 = @(y) g - sign(y).*(cd/m).*y.^2-(k/m).*(x-L)-(gamma/m).*y;

for i = 1:n
    k1 = h*f1(y(i));
    z1 = h*f2(y(i));
    k2 = h*f1(y(i)+(z1/2));
    z2 = h*f2(y(i)+(z1/2));
    k3 = h*f1(y(i)+(z2/2));
    z3 = h*f2(y(i)+(z2/2));
    k4 = h*f1(y(i)+(z3));
    z4 = h*f2(y(i)+(z3));
    xout = x+((k1+(k2*2)+(k3*2)+k4)/6)*h
    yout = y+((z1+(z2*2)+(2*z3)+z4)/6)*h
 
end
xout
yout
plot(xout,t,yout,t)

It's giving an error :

"Index exceeds the number of array elements (1)."

and I've tried everything I canm but I can't seem to find what I need to fix it.

Answer 1

[Walter Roberson]

Where do you assign to y(2)?

 y(i+1) = yout

Likewisefor x

[Steve Avilesmejia]

im not sure what you are asking?

[Walter Roberson]

After you assign to yout

x(i+1) = xout;
y(i+1) = yout;

[Steve Avilesmejia]

are you putting this at the end in or out of the if statment.

Answer 2

[Pier Giorgio Petrolini]

Hello Steve,

Is it that what you were looking for ?

clear all , clc, format compact
L = 30;
m = 68.1;
c_d = 0.25;
k = 40;
gamma = 8;
g = 9.81;
to = 0;
tf = 50;
h = 0.1;
n = ((tf-to)/h);
t = to:h:tf;
x1 = 0;
y1 = 0;
x=x1;
y=y1;
f1 = @(y) g - sign(y).*cd/m.*y.^2;
f2 = @(y) g - sign(y).*(cd/m).*y.^2-(k/m).*(x-L)-(gamma/m).*y;

for y = 1:n
	k1 = h*f1(y);
	z1 = h*f2(y);
	k2 = h*f1(y+(z1/2));
	z2 = h*f2(y+(z1/2));
	k3 = h*f1(y+(z2/2));
	z3 = h*f2(y+(z2/2));
	k4 = h*f1(y+(z3));
	z4 = h*f2(y+(z3));
	xout(y,:) = x+((k1+(k2*2)+(k3*2)+k4)/6)*h;
	yout(y,:) = y+((z1+(z2*2)+(2*z3)+z4)/6)*h;	
end

% Starting Point
k1 = h*f1(y);
z1 = h*f2(y);
k2 = h*f1(y+(z1/2));
z2 = h*f2(y+(z1/2));
k3 = h*f1(y+(z2/2));
z3 = h*f2(y+(z2/2));
k4 = h*f1(y+(z3));
z4 = h*f2(y+(z3));
xout0 = x+((k1+(k2*2)+(k3*2)+k4)/6)*h;
yout0= y+((z1+(z2*2)+(2*z3)+z4)/6)*h;

xout = [xout0; xout]
yout = [yout0; yout]

plot(xout,t,yout,t)

[Walter Roberson]

No, after the loop, you always call f1(y) but y is not changing in the code there, so the value effectively becomes a constant when instead it needs to update to current position.

[Steve Avilesmejia]

whoa, umm the plot should be (xout,t,yout,t)

[Walter Roberson]

No it should not be. yout and xout are scalars in your code, you cannot plot them.

You are doing a Runge Kutta calculation. Each iteration estimates a new x y to feed back to the next iteration. When you use x(i) and y(i) you need to use the x and y that were calculated from the previous iteration, so you need to be storing into x and y. It is those recorded values that need to be plotted.

Answer 3

[Image Analyst]

From this code (before the loop):

y1 = 0;
x=x1;
y=y1;

we can see that y is a simple scalar, not an array, with a value of 0. Then you enter the loop. You have this line of code:

k1 = h*f1(y(i));

so you're trying to pass y(1), then y(2), then y(3), etc. into the f1 function. This works fine for the first iteration because y(1) still works when y is a scalar. However on the second iteration, you're trying to get y(2). Well....there is no second element of y, so you get the error.

Not sure what you want to do so not sure how to tell you to fix it. Why did you not put any comments in your code? If a programmer did that who worked for me, I'd "retrain" him. Comments are essential. In fact Steve Lord suggests that you FIRST write the whole program in comments and only then begin to fill in the lines below the comments with actual MATLAB code. Please comment your code and tell us what it is expected to do.

Answer 4

[James Tursa]

Sorry for the late reply ... I just saw this post.

These lines

f1 = @(y) g - sign(y).*cd/m.*y.^2;
f2 = @(y) g - sign(y).*(cd/m).*y.^2-(k/m).*(x-L)-(gamma/m).*y;

clearly indicate that the derivatives depend on both the current states of x and y, not just y. So you have a fundamental flaw in the construction of the derivative function handles and things will never work. Since you have two states, x and y, you should pass in both states to your function handles. To be consistent, I would do this for both f1 and f2 even though x is not in f1. E.g.,

f1 = @(x,y) g - sign(y).*cd/m.*y.^2;
f2 = @(x,y) g - sign(y).*(cd/m).*y.^2-(k/m).*(x-L)-(gamma/m).*y;

Having said that, it does lead me to wonder if you have these equations coded correctly. Why does f1 depend only on y but f2 depends on both x and y? You should double check this, and maybe post an image of the differential equations you are solving so we can verify this. I suspect that what you have for this might not be correct.

Your looping would then look something like this, using the current values of x and y (which are x(i) and y(i)) for your derivative calls:

n = numel(t);
x = zeros(size(t));
y = zeros(size(t));
x(1) = x1;
y(1) = y1;
for i = 1:n-1
    k1 = h * f1( x(i)       , y(i)        );
    z1 = h * f2( x(i)       , y(i)        );
    k2 = h * f1( x(i) + k1/2, y(i) + z1/2 );
    z2 = h * f2( x(i) + k1/2, y(i) + z1/2 );
    k3 = h * f1( x(i) + k2/2, y(i) + z2/2 );
    z3 = h * f2( x(i) + k2/2, y(i) + z2/2 );
    k4 = h * f1( x(i) + k3  , y(i) + z3   );
    z4 = h * f2( x(i) + k3  , y(i) + z3   );
    x(i+1) = x(i) + (k1 + 2*k2 + 2*k3 + k4)/6; % no h factor here because it is already part of the k's
    y(i+1) = y(i) + (z1 + 2*z2 + 2*z3 + z4)/6; % no h factor here because it is already part of the z's
 end

Note that by using spacing I have made the code much more readable and easier to debug.

The variables to plot will be the vectors x and y.

Finally, you should put comments on every line with a constant, stating the units of that constant. This makes it much easier to find unit mismatch errors and helps the reader decipher what the code is doing.