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

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 

[Dyuman Joshi]

Are you sure the formula for calculating total steps (ndt) is correct?

You are dividing vector by another vector. The backlash operator performs matrix division.

KC  = 2:2:20;   
T   = (KC.*0.1);   % period due to KC 
Om  = (2*pi)./T;   % omega 
dt  = T/40;        % time step
ndt = (T/dt)*10;   % total steps to be calculated

[Andre Khosrow]

Yes. I am trying to calculate V with KC = 2,4,6,8,10....20. For instance when KC = 2 -> T = 0.2 -> Om = 31.42, dt = 0.005 (I calculated, via hand calculation, the ndt and figured that it doesnt change for any value of T so it will always be 400).

Then I want it to undergo runge kutta and get a value for V. After that is done, I want KC=4, then undergo the same steps and get a V value for KC=4. etc until it reaches KC=20 and have 10 values for V.

[Dyuman Joshi]

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]

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.

Answer 1

[Walter Roberson]

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;

[Walter Roberson]

Note:

If you are doing non-adaptive steps, if the steps are at fixed intervals, and you want to evaluate the effect of using different step sizes, then in the special case where the different candidate step sizes are all at integral multiples of the smallest step size, then you can start by evaluating at the smallest step size, recording the values in an array, and then after that instead of re-evaluating at larger step sizes, instead using appropriate indexing to pull back the already-calculated results.

Or.. you could use memoize

[Andre Khosrow]

I did as you suggested an ended up with this;

% Code for simple initial value problem of offshore energy extraction
% device. 
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
T = zeros(1,10)         ;
Om = zeros(1,10)        ;
dt = zeros(1,10)        ;
ndt= zeros(1,10)        ;
KC = 2:2:20             ;

for k= 1:10            
    T(k)   = (KC(k)*0.1)/1        ;  
    Om(k)  = (2*P)/T(k)           ;
    dt(k)  = T(k)/40              ;
    ndt(k) = (T(k)/dt(k))*10      ;
end  

X(1:ndt(k)+1) = 0                 ;
V(1:ndt(k)+1) = 0                 ;
time(1:ndt(k)+1)=(0:ndt(k))*dt(k) ;

for n = 1:ndt(k) 
    % to calculate kx1 and kv1
    ta = time(n);
    aw = Um*Om(k)*(cos(Om(k)*ta)-((2/3)*cos(2*Om(k)*ta))); % acceleration of the fluid
    Vw = Um*sin(Om(k)*ta)-((Um/3)*sin(2*Om(k)*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(k)/2;
    aw = Um*Om(k)*(cos(Om(k)*ta)-((2/3)*cos(2*Om(k)*ta))); % acceleration of the fluid
    Vw = Um*sin(Om(k)*ta)-((Um/3)*sin(2*Om(k)*ta))    ; % velocity of the fluid
    Va = V(n)+0.5*dt(k)*kv1;
    Xa = X(n)+0.5*dt(k)*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(k)*kv1;
    kv2 = (t1+t2+t3)/(m+Ca*md);
    % to calculate kx3 and kv3
    ta = time(n)+dt(k)/2;
    aw = Um*Om(k)*(cos(Om(k)*ta)-((2/3)*cos(2*Om(k)*ta))); % acceleration of the fluid
    Vw = Um*sin(Om(k)*ta)-((Um/3)*sin(2*Om(k)*ta))    ; % velocity of the fluid
    Va = V(n)+0.5*dt(k)*kv2;
    Xa = X(n)+0.5*dt(k)*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(k)*kv2;
    kv3 = (t1+t2+t3)/(m+Ca*md);
    % to calculate kx4 and kv4
    ta = time(n)+dt(k)/2;
    aw = Um*Om(k)*(cos(Om(k)*ta)-((2/3)*cos(2*Om(k)*ta))); % acceleration of the fluid
    Vw = Um*sin(Om(k)*ta)-((Um/3)*sin(2*Om(k)*ta))    ; % velocity of the fluid
    Va = V(n)+0.5*dt(k)*kv3;
    Xa = X(n)+0.5*dt(k)*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(k)*kv3;
    kv4 = (t1+t2+t3)/(m+Ca*md);
    % next step values
    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);
end 

I can see that the loop is working but I'm still getting 1 value of X(n+1) & V(n+1) instead of 10.

[Walter Roberson]

for k= 1:10

after the end of a for loop, the loop variable (k) will be left at the last value it was assigned, which would be 10 in this case. Your k loop is only a few lines, and after that you effectively have

k = 10;

in effect.

For your purposes you need to have a for k = 1:10 loop and you need to assign to

    X{k}(n+1) = X(n)+(dt(k)/6)*(kx1+2*kx2+2*kx3+kx4);
    V{k}(n+1) = V(n)+(dt(k)/6)*(kv1+2*kv2+2*kv3+kv4);

having initialized with

X{k}(1:ndt(k)+1) = 0;
V{k}(1:ndt(k)+1) = 0;

Possibly you might also need to make time into a cell array as well.

[Andre Khosrow]

Sorry for late reply, basically I got what I wanted (tweaked the original code) until the comment of calculating the next step values. I get the error "Unable to perform assignment because the left and right sides have a different number of elements". I managed to get the 10 values for kx1,kx2 .. kx4 but I just need next step values to calculate for each value of kx1,kx2 etc. I tried what you suggested but it didn't help with the original code.

% Code for simple initial value problem of offshore energy extraction
% device. 
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
T = zeros(1,10)         ;
Om = zeros(1,10)        ;
dt = zeros(1,10)        ;
ndt= zeros(1,10)        ;
KC = 2:2:20             ;

for k= 1:10            
    T(k)   = (KC(k)*0.1)/1        ;  
    Om(k)  = (2*P)/T(k)           ;
    dt(k)  = T(k)/40              ;
    ndt(k) = (T(k)/dt(k))*10      ;
end  

X(1:ndt(k)+1) = 0                 ;
V(1:ndt(k)+1) = 0                 ;
Power(1:ndt(k)+1)= 0              ;
time(1:ndt(k)+1)=(0:ndt(k))*dt(k) ;
for n = 1:ndt(k) 
    % 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(k)/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(k)*kv1;
    Xa = X(n)+0.5*dt(k)*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(k)*kv1;
    kv2 = (t1+t2+t3)/(m+Ca*md);
    % to calculate kx3 and kv3
    ta = time(n)+dt(k)/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(k)*kv2;
    Xa = X(n)+0.5*dt(k)*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(k)*kv2;
    kv3 = (t1+t2+t3)/(m+Ca*md);
    % to calculate kx4 and kv4
    ta = time(n)+dt(k)/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(k)*kv3;
    Xa = X(n)+0.5*dt(k)*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(k)*kv3;
    kv4 = (t1+t2+t3)/(m+Ca*md);
    
    % next step values trying to fix the error here

    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);
    Power(n+1)= c*V(n+1)^2                       ;
end 

Unable to perform assignment because the left and right sides have a different number of elements.

figure (01);
plot(time(1:ndt(k)),X(1:ndt(k)),'-r');
xlabel('time (s)'); %honrizontal axis's label is 'x'
ylabel('X (m)'); %vertical axis's label is 'T'
legend ('RK method' )

figure (02);
plot(time(1:ndt(k)),V(1:ndt(k)),'-r');
xlabel('time (s)'); %honrizontal axis's label is 'x'
ylabel('V (m)'); %vertical axis's label is 'T'
legend ('RK method' )

figure (03)
plot(time(1:ndt(k)),Power(1:ndt(k)),'-r');
xlabel('time (s)'); %honrizontal axis's label is 'x'
ylabel('Power'); 
legend ('Power')

[Walter Roberson]

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]

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                                           ;