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

4 views (last 30 days)
Andre Khosrow on 17 Mar 2023
Edited: Andre Khosrow on 18 Apr 2023
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
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.

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;
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 ;