Code not coming out of if statement.
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Respuestas (4)
Walter Roberson
el 7 de Dic. de 2020
Where do you assign to y(2)?
y(i+1) = yout
Likewisefor x
3 comentarios
Walter Roberson
el 7 de Dic. de 2020
After you assign to yout
x(i+1) = xout;
y(i+1) = yout;
Pier Giorgio Petrolini
el 7 de Dic. de 2020
Editada: Image Analyst
el 7 de Dic. de 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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/452698/image.png)
3 comentarios
Walter Roberson
el 7 de Dic. de 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.
Image Analyst
el 7 de Dic. de 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 comentarios
James Tursa
el 11 de Dic. de 2020
Editada: James Tursa
el 14 de Dic. de 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 comentarios
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!