MATLAB Answers

Code not coming out of if statement.

25 views (last 30 days)
Steve Avilesmejia
Steve Avilesmejia on 7 Dec 2020
Edited: James Tursa on 14 Dec 2020
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.

  0 Comments

Sign in to comment.

Answers (4)

Walter Roberson
Walter Roberson on 7 Dec 2020
Where do you assign to y(2)?
y(i+1) = yout
Likewisefor x

  3 Comments

Steve Avilesmejia
Steve Avilesmejia on 7 Dec 2020
im not sure what you are asking?
Walter Roberson
Walter Roberson on 7 Dec 2020
After you assign to yout
x(i+1) = xout;
y(i+1) = yout;
Steve Avilesmejia
Steve Avilesmejia on 7 Dec 2020
are you putting this at the end in or out of the if statment.

Sign in to comment.


Pier Giorgio Petrolini
Pier Giorgio Petrolini on 7 Dec 2020
Edited: Image Analyst on 7 Dec 2020
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)

  3 Comments

Walter Roberson
Walter Roberson on 7 Dec 2020
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
Steve Avilesmejia on 7 Dec 2020
whoa, umm the plot should be (xout,t,yout,t)
Walter Roberson
Walter Roberson on 7 Dec 2020
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.

Sign in to comment.


Image Analyst
Image Analyst on 7 Dec 2020
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.

  0 Comments

Sign in to comment.


James Tursa
James Tursa on 11 Dec 2020
Edited: James Tursa on 14 Dec 2020
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.

  0 Comments

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!