Clear Filters
Clear Filters

How do I fix the error of arrays having incompatible sizes?

5 views (last 30 days)
I'm trying to calculate 4th order runge kutta method to obtain a velocity and I need to do it by having KC in intervals. I tried using floor but not sure if that even changes my error. Should I use a for loop for the KC interval ?
Um = 1 ; % amplitude of oscillatory flow in m/s
D = 0.1 ; % diameter of the cylinder in m
m = 50 ; % in kg
Rho = 1024 ; % density of the water kg/m3
K = 200 ; % total spring stiffness in N/m
c = 100 ; % d6amping coefficient
Ca = 1 ; % added mass coefficient
Cd = 1.8 ; % drag coefficient
L = 1 ; % length of the cylinder
Ap = D*L ; % projected area
P = pi ;
md = Rho*P*(D*D)/(4*L) ; % added mass
% for part 4 KC has intervals
KC = 2:2:20 ;
T = (KC.*0.1)/1 ; % period due to KC
Om = (2*P)./T ; % omega
dt = T/40 ; % time step
ndt = (T/dt)*10 ; % total steps to be calculated
X(1:floor(ndt+1)) = 0 ;
V(1:floor(ndt+1)) = 0 ;
time(1:floor(ndt+1))=(0:floor(ndt)).*dt;
for n = 1:ndt
% to calculate kx1 and kv1
ta = time(n);
aw = Um*Om.*(cos(Om.*ta)-((2/3)*cos(2*Om.*ta))); % acceleration of the fluid
Vw = Um*sin(Om.*ta)-((Um/3)*sin(2*Om.*ta)) ; % velocity of the fluid
Va = V(n);
Xa = X(n);
t1 = Ca*md*aw;
t2 = 0.5*Rho*Cd*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx1 = V(n);
kv1 = (t1+t2+t3)/(m+Ca*md);
% to calculate kx2 and kv2
ta = time(n)+dt./2;
aw = Um*Om.*(cos(Om.*ta)-((2/3)*cos(2*Om.*ta))); % acceleration of the fluid
Vw = Um*sin(Om.*ta)-((Um/3)*sin(2*Om.*ta)) ; % velocity of the fluid
Va = V(n)+0.5*dt.*kv1;
Xa = X(n)+0.5*dt.*kx1;
t1 = Ca*md*aw;
t2 = 0.5*Rho*Cd*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx2 = V(n)+ 0.5*dt.*kv1;
kv2 = (t1+t2+t3)/(m+Ca*md);
% to calculate kx3 and kv3
ta = time(n)+dt./2;
aw = Um*Om*(cos(Om*ta)-((2/3)*cos(2*Om*ta))); % acceleration of the fluid
Vw = Um*sin(Om*ta)-((Um/3)*sin(2*Om*ta)) ; % velocity of the fluid
Va = V(n)+0.5*dt*kv2;
Xa = X(n)+0.5*dt*kx2;
t1 = Ca*md*aw;
t2 = 0.5*Rho*Cd*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx3 = V(n)+0.5*dt*kv2;
kv3 = (t1+t2+t3)/(m+Ca*md);
% to calculate kx4 and kv4
ta = time(n)+dt/2;
aw = Um*Om*(cos(Om*ta)-((2/3)*cos(2*Om*ta))); % acceleration of the fluid
Vw = Um*sin(Om*ta)-((Um/3)*sin(2*Om*ta)) ; % velocity of the fluid
Va = V(n)+0.5*dt*kv3;
Xa = X(n)+0.5*dt*kx3;
t1 = Ca*md*aw;
t2 = 0.5*Rho*Cd*Ap*abs(Vw-Va)*(Vw-Va);
t3 = -K*Xa-c*Va;
kx4 = V(n)+dt*kv3;
kv4 = (t1+t2+t3)/(m+Ca*md);
% next step values
X(n+1) = X(n)+(dt/6)*(kx1+2*kx2+2*kx3+kx4);
V(n+1) = V(n)+(dt/6)*(kv1+2*kv2+2*kv3+kv4);
end
  4 Comments
Dyuman Joshi
Dyuman Joshi on 20 Mar 2023
The problem is here -
%no need to use floor as ndt is an integer
ndt = (T/dt)*10 ; %400
time(1:ndt+1)=(0:ndt).*dt;
0:400 has 401 elements, where as dt has 10 elements, so you can not do element-wise multiplicaiton.
Which is why I asked you if ndt=(T/dt)*10 is correct or not.
You can not change dt as KC is fixed (as you mentioned the range in your comment above), and dt depends on KT. And you mentioned ndt is always 400, so 0:ndt is fixed as well.
Andre Khosrow
Andre Khosrow on 21 Mar 2023
Edited: Andre Khosrow on 21 Mar 2023
I see, so is there anyway to calculate it with a for loop that eliminates the array imcompatibility ? Or is it easier to just calculate it separately ? [i.e have KC = 2 which would make dt fixed with 1 element, obtain the V and then manually change KC=4 and obtain V etc.

Sign in to comment.

Answers (1)

Walter Roberson
Walter Roberson on 21 Mar 2023
Yes, you need to loop over KC values.
Your ndt is a vector with values calculated in terms of KC. And that means that
X(1:floor(ndt+1)) = 0 ;
is not going to work the way you want. With ndt being a vector, that is not going to create an row of X for each different value in ndt, with the rows being different lengths -- not going to happen. Instead, when you have a : operator that has a non-scalar value as one of the parameters to the operator, then the first element is going to be used and the rest ignored, as if you had written
X(1:floor(ndt(1)+1)) = 0;
  6 Comments
Walter Roberson
Walter Roberson on 17 Apr 2023
Om = zeros(1,10) ;
Om is a vector of length 10.
aw = Um*Om.*(cos(Om.*ta)-((2/3)*cos(2*Om.*ta))); % acceleration of the fluid
Vw = Um*sin(Om.*ta)-((Um/3)*sin(2*Om.*ta)) ; % velocity of the fluid
Va = V(n);
aw is defined in terms of the vector Om, so aw is a vector. Likewise for Vw. But at this point, Va is assigned a scalar.
t1 = Ca*md*aw;
t2 = 0.5*Rho*Cd*Ap.*abs(Vw-Va).*(Vw-Va);
aw and Vw are vectors so t1 and t2 are vectors.
kv1 = (t1+t2+t3)/(m+Ca*md);
and that makes kv1 a vector as well.
aw = Um*Om.*(cos(Om.*ta)-((2/3)*cos(2*Om.*ta))); % acceleration of the fluid
Vw = Um*sin(Om.*ta)-((Um/3)*sin(2*Om.*ta)) ; % velocity of the fluid
For unknown reasons, aw and Vw are recalculated. But for the same reasons as before they become vectors.
Va = V(n)+0.5*dt(k)*kv1;
Va is also recalculated, but this time it is defined in terms of vector kv1, so Va now becomes a vector.
t1 = Ca*md*aw;
t2 = 0.5*Rho*Cd*Ap.*abs(Vw-Va).*(Vw-Va);
t1 and t2 are recalculated, but for the same reasons as before they become vectors. The differences is that before Va was a scalar and now Va is a vector so the resulting t1 might be different than the first time.
Anyhow, you go on and recalculate variables several times, but t1 and t2 continue to be vectors so
kv3 = (t1+t2+t3)/(m+Ca*md);
is going to be a vector, and several kx* variables are going to be vectors.
X(n+1) = X(n)+(dt(k)/6)*(kx1+2*kx2+2*kx3+kx4);
V(n+1) = V(n)+(dt(k)/6)*(kv1+2*kv2+2*kv3+kv4);
but with kx* and kv* variables being vectors, the right hand side of those two computations are going to be vectors, but you are assigning to a scalar location in each case.
Andre Khosrow
Andre Khosrow on 18 Apr 2023
Edited: Andre Khosrow on 18 Apr 2023
Ohh I see. How would I fix it then? So I just took the first value of and ran it through. I'll do the same for the other 9 values.
X(n+1) = X(n)+(dt(k)/6)*(kx1+2*kx2(1:1)+2*kx3(1:1)+kx4(1:1)) ;
V(n+1) = V(n)+(dt(k)/6)*(kv1(1:1)+2*kv2(1:1)+2*kv3(1:1)+kv4(1:1));
Power(n+1)= c*V(n+1)^2 ;

Sign in to comment.

Categories

Find more on Fluid Mechanics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!