Code not coming out of if statement.

1 次查看（过去 30 天）
Steve Avilesmejia
Steve Avilesmejia 2020-12-7
编辑： James Tursa 2020-12-14
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.

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

回答（4 个）

Walter Roberson
Walter Roberson 2020-12-7
Where do you assign to y(2)?
y(i+1) = yout
Likewisefor x
Pier Giorgio Petrolini
编辑：Image Analyst 2020-12-7
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 个评论
显示 1更早的评论
Steve Avilesmejia
Steve Avilesmejia 2020-12-7
whoa, umm the plot should be (xout,t,yout,t)
Walter Roberson
Walter Roberson 2020-12-7
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
Image Analyst 2020-12-7
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.
James Tursa
James Tursa 2020-12-11
编辑：James Tursa 2020-12-14
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.

类别

在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by