25 views (last 30 days)

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.

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)

Walter Roberson
on 7 Dec 2020

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.

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.

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.

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

Start Hunting!
## 0 Comments

Sign in to comment.