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

3 visualizaciones (últimos 30 días)
josh johnson
josh johnson el 17 de Nov. de 2022
Comentada: Torsten el 18 de Nov. de 2022
Ran in:
Im try to plot t,y. but i keep getting the effor.
Unable to perform assignment because the
left and right sides have a different
number of elements.
Error in testq3>MDOFODE45 (line 161)
yy(6) = (f(3) - k(9)*y(6) - k(6)*y(5) -
.2*k(6)*y(5) + .2*k(7)*y(4) - .2*k*(9)*y(6)
)/m(9);
Error in
testq3>@(t,y)MDOFODE45(t,y,m,k,Fd,FI,w,phiB,phiAC,phiD)
(line 116)
[t,y] = ode45(@(t,y) MDOFODE45(t, y, m, k,
Fd, FI, w , phiB, phiAC, phiD), tspan,y0) ;
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I
sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name,
ode, tspan, y0, options, varargin);
Error in testq3 (line 116)
[t,y] = ode45(@(t,y) MDOFODE45(t, y, m, k,
Fd, FI, w , phiB, phiAC, phiD), tspan,y0) ;
>>
clc;clear;
%---------------------------------variable -------------------
%distances to centre of gravity
I1 = 16; %distance to g side 1
I2 = 14; %distance to g side 1
I3 = 17; %distance to g side 2
I4 = 13; %distance to g side 2
M =9000*1000 ; %mass 9000 tonnes
I=1.35*10^9; %interia
Rout= 1.25 ; %outter radius of legs
Rin = Rout-.02; %inner raidus of legs
E=204*10^9; %youngs modulus
Im=(pi/4)*(Rout^4-Rin^4); %second moment of cylinder
dis1 = 41.42; %distance to pillar D
dis2 = dis1 / 2; %distance to pillar A & C
%------------------------------Case selection ------------------
% i = input('Enter input:');
i=1 % modified to run here
i = 1
switch i
case 1 %variables for sea state 1
w = 1.571; %frequency
T = 4; %Wave period
Wh = 2 ; %wave height
Fd = 3 ; %Max wave Drag (kN)
FI = 60 ; %Max interia wave (kN)
Lo = 24.956; %Deep water wave
case 2 %variables for sea state 2
w = 1.142;
T = 4;
Wh = 2 ;
Fd = 3 ;
FI = 60 ;
Lo = 47.183;
case 3 %variables for sea state 3
w = 0.898;
T = 4;
Wh = 2 ;
Fd = 3 ;
FI = 60 ;
Lo = 76.428;
case 4 %variables for sea state 2
w = 0.698;
T = 4;
Wh = 2 ;
Fd = 3 ;
FI = 60 ;
Lo = 126.341;
end
if (i == 1) %variables for sea state 1
wd = 10:5:30; % wave depth
L = Lo*sqrt(tanh(2*pi*wd/Lo)); % intermediate water waelength
wps = L / T; % wave proagation speed
keq=(3*E*Im)./(wd.^3); %lateral stiffness
elseif( i == 2 ) %variables for sea state 2
wd = 10:5:30;
L = Lo*sqrt(tanh(2*pi*wd/Lo));
wps = L / T;
keq=(3*E*Im)./(wd.^3);
elseif (i == 3) %variables for sea state 3
wd = 20:5:40;
L = Lo*sqrt(tanh(2*pi*wd/Lo));
wps = L / T;
keq=(3*E*Im)./(wd.^3);
elseif(i == 4) %variables for sea state 4
wd = 20:5:40;
L = Lo*sqrt(tanh(2*pi*wd/Lo));
wps = L / T;
keq=(3*E*Im)./(wd.^3);
end
%phi matix space allocation
phiB = zeros(1,5);
phiAC = zeros(1,5);
phiD = zeros(1,5);
for B = 1:5
pillart = 0/wps(B) ; %time to pillar B
pillartAC = dis2/wps(B) ; %time to pillar A & C
pillartD = dis1/wps(B) ; %time to pillar D
phiB(B) = (2*pi*pillart)/T; %phi at pillar B
phiAC(B) = (2*pi*pillartAC)/T; %phi at pillar A & C
phiD(B) = (2*pi*pillartD)/T; %phi at pillar D
end
%------------------------------------------matrices-------------
kmat=[];
for kx = 1:length(L)
keq=(3*E*Im)/(L(kx)^3); %lateral stiffness
kmat=[keq,kmat];
k=[4*keq 0 2*keq*(I1-I2); 0 4*keq 2*keq*(I3-I4); 2*keq*(I1-I2) 2*keq*(I3-I4) 2*keq*(I1^2+I2^2+I3^2+I4^2)];
end
c =0.2 * k;
m =[M 0 0;0 M 0;0 0 I ];
disp(k);
1.0e+09 * 0.0189 0 0.0189 0 0.0189 0.0377 0.0189 0.0377 8.5856
%ODE45 solver
y0=[0 0 0 0 0 0];% inital conditions
tspan=[0 10]; % time span
[t,y] = ode45(@(t,y) MDOFODE45(t, y, m, k, Fd, FI, w , phiB, phiAC, phiD), tspan,y0) ;
%results
% y_ = y(:, 4:end) ;
% ydot_ = y(:, 1:3) ;
figure(1)
plot (t,y(:,1)) % Plot the time
xlabel('Time (s)') % Label the plot
ylabel('Response X_1')
figure(2)
plot (t,y(:,2)) % Plot the time
xlabel('Time (s)') % Label the plot
ylabel('Response X_2')
% ODE Function
function yy = MDOFODE45(t, y, m, k, Fd, FI, w , phiB, phiAC, phiD)
I1 = 16;
I2 = 14;
I3 = 17;
I4 = 13;
FlegB = (3*Fd/4)*(sin(w*t+phiB)) - (3*sin(w*t+phiB)/3)+FI*cos(w*t+phiB); %forces at leg B
FlegAC = (3*Fd/4)*(sin(w*t+phiAC)) - (3*sin(w*t+phiAC)/3)+FI*cos(w*t+phiAC); %forces at leg A & C
FlegD = (3*Fd/4)*(sin(w*t+phiD)) - (3*sin(w*t+phiD)/3)+FI*cos(w*t+phiD); %forces at leg D
FlegBx = FlegB * sind(45); %forces at leg B in x direction
FlegBy = FlegB * cosd(45); %forces at leg B in y direction
FlegACx = FlegAC * sind(45); %forces at leg A & C in x direction
FlegACy = FlegAC * cosd(45); %forces at leg A & C in y direction
FlegDx = FlegD * sind(45); %forces at leg D in x direction
FlegDy = FlegD * cosd(45); %forces at leg D in y direction
%force matrix
f=[FlegACx+FlegBx+FlegACx+FlegDx ; FlegACy+FlegBy+FlegACy+FlegDy ; (FlegACx+FlegBx)*I1-(FlegACx+FlegDx)*I2-(FlegACy+FlegDy)*I3+(FlegBy+FlegACy)*I4];
yy = zeros(6,1);
yy(1) = y(4) ;
yy(2) = y(5) ;
yy(3) = y(6) ;
yy(4) = (f(1) - k(1)*y(1) - k(3)*y(6) - .2*k(1)*y(4) - .2*k(3)*y(6) )/m(1);
yy(5) = (f(2) - k(5)*y(2) - k(6)*y(3) - .2*k(5)*y(5) - .2*k(6)*y(6) )/m(5) ;
yy(6) = (f(3) - k(9)*y(6) - k(6)*y(5) - .2*k(6)*y(5) + .2*k(7)*y(4) - .2*k(9)*y(6) )/m(9);
end
I have also noticed that my k fuction doesnt doesnt full work because the k mat isnt fully statisfied by the keq values. any help would be greatly apprieated
  3 comentarios
Mostrar 1 comentario más antiguo
josh johnson
josh johnson el 18 de Nov. de 2022
Much apprecitaed !!
is there a way to get the k to run for each wave depth (wd) then have figures that show the response at each depth ?
Torsten
Torsten el 18 de Nov. de 2022
Make a loop over k, call ODE45 in the loop for each value of k separately and save the results.

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuesta aceptada

Cris LaPierre
Cris LaPierre el 17 de Nov. de 2022
Editada: Cris LaPierre el 17 de Nov. de 2022
Ran in:
The problem is that your calculation for yy(6) returns a 3x3 matrtix, and you are assigning it to a 1x1 space. Hence, the left (1x1) and the right (3x3) have a different number of elements.
I think this is your mistake
(f(3) - k(9)*y(6) - k(6)*y(5) - .2*k(6)*y(5) + .2*k(7)*y(4) - .2*k*(9)*y(6) )/m(9)
% ^
I think you want .2*k(9)*y(6) instead
clc;clear;
%---------------------------------variable -------------------
%distances to centre of gravity
I1 = 16; %distance to g side 1
I2 = 14; %distance to g side 1
I3 = 17; %distance to g side 2
I4 = 13; %distance to g side 2
M =9000*1000 ; %mass 9000 tonnes
I=1.35*10^9; %interia
Rout= 1.25 ; %outter radius of legs
Rin = Rout-.02; %inner raidus of legs
E=204*10^9; %youngs modulus
Im=(pi/4)*(Rout^4-Rin^4); %second moment of cylinder
dis1 = 41.42; %distance to pillar D
dis2 = dis1 / 2; %distance to pillar A & C
%------------------------------Case selection ------------------
% i = input('Enter input:');
i=1 % modified to run here
i = 1
switch i
case 1 %variables for sea state 1
w = 1.571; %frequency
T = 4; %Wave period
Wh = 2 ; %wave height
Fd = 3 ; %Max wave Drag (kN)
FI = 60 ; %Max interia wave (kN)
Lo = 24.956; %Deep water wave
case 2 %variables for sea state 2
w = 1.142;
T = 4;
Wh = 2 ;
Fd = 3 ;
FI = 60 ;
Lo = 47.183;
case 3 %variables for sea state 3
w = 0.898;
T = 4;
Wh = 2 ;
Fd = 3 ;
FI = 60 ;
Lo = 76.428;
case 4 %variables for sea state 2
w = 0.698;
T = 4;
Wh = 2 ;
Fd = 3 ;
FI = 60 ;
Lo = 126.341;
end
if (i == 1) %variables for sea state 1
wd = 10:5:30; % wave depth
L = Lo*sqrt(tanh(2*pi*wd/Lo)); % intermediate water waelength
wps = L / T; % wave proagation speed
keq=(3*E*Im)./(wd.^3); %lateral stiffness
elseif( i == 2 ) %variables for sea state 2
wd = 10:5:30;
L = Lo*sqrt(tanh(2*pi*wd/Lo));
wps = L / T;
keq=(3*E*Im)./(wd.^3);
elseif (i == 3) %variables for sea state 3
wd = 20:5:40;
L = Lo*sqrt(tanh(2*pi*wd/Lo));
wps = L / T;
keq=(3*E*Im)./(wd.^3);
elseif(i == 4) %variables for sea state 4
wd = 20:5:40;
L = Lo*sqrt(tanh(2*pi*wd/Lo));
wps = L / T;
keq=(3*E*Im)./(wd.^3);
end
%phi matix space allocation
phiB = zeros(1,5);
phiAC = zeros(1,5);
phiD = zeros(1,5);
for B = 1:5
pillart = 0/wps(B) ; %time to pillar B
pillartAC = dis2/wps(B) ; %time to pillar A & C
pillartD = dis1/wps(B) ; %time to pillar D
phiB(B) = (2*pi*pillart)/T; %phi at pillar B
phiAC(B) = (2*pi*pillartAC)/T; %phi at pillar A & C
phiD(B) = (2*pi*pillartD)/T; %phi at pillar D
end
%------------------------------------------matrices-------------
kmat=[];
for kx = 1:length(L)
keq=(3*E*Im)/(L(kx)^3); %lateral stiffness
kmat=[keq,kmat];
k=[4*keq 0 2*keq*(I1-I2); 0 4*keq 2*keq*(I3-I4); 2*keq*(I1-I2) 2*keq*(I3-I4) 2*keq*(I1^2+I2^2+I3^2+I4^2)];
end
c =0.2 * k;
m =[M 0 0;0 M 0;0 0 I ];
disp(k);
1.0e+09 * 0.0189 0 0.0189 0 0.0189 0.0377 0.0189 0.0377 8.5856
%ODE45 solver
y0=[0 0 0 0 0 0];% inital conditions
tspan=[0 500]; % time span
[t,y] = ode45(@(t,y) MDOFODE45(t, y, m, k, Fd, FI, w , phiB, phiAC, phiD), tspan,y0) ;
%results
% y_ = y(:, 4:end) ;
% ydot_ = y(:, 1:3) ;
figure(1) % ##################### <--------------- modified this line
plot (t,y(:,1)) % Plot the time
xlabel('Time (s)') % Label the plot
ylabel('Response X_1')
figure(2) % ##################### <--------------- modified this line
plot (t,y(:,2)) % Plot the time
xlabel('Time (s)') % Label the plot
ylabel('Response X_2')
% ODE Function
function yy = MDOFODE45(t, y, m, k, Fd, FI, w , phiB, phiAC, phiD)
I1 = 16;
I2 = 14;
I3 = 17;
I4 = 13;
FlegB = (3*Fd/4)*(sin(w*t+phiB)) - (3*sin(w*t+phiB)/3)+FI*cos(w*t+phiB); %forces at leg B
FlegAC = (3*Fd/4)*(sin(w*t+phiAC)) - (3*sin(w*t+phiAC)/3)+FI*cos(w*t+phiAC); %forces at leg A & C
FlegD = (3*Fd/4)*(sin(w*t+phiD)) - (3*sin(w*t+phiD)/3)+FI*cos(w*t+phiD); %forces at leg D
FlegBx = FlegB * sin(45); %forces at leg B in x direction
FlegBy = FlegB * cos(45); %forces at leg B in y direction
FlegACx = FlegAC * sin(45); %forces at leg A & C in x direction
FlegACy = FlegAC * cos(45); %forces at leg A & C in y direction
FlegDx = FlegD * sin(45); %forces at leg D in x direction
FlegDy = FlegD * cos(45); %forces at leg D in y direction
%force matrix
f=[FlegACx+FlegBx+FlegACx+FlegDx ; FlegACy+FlegBy+FlegACy+FlegDy ; (FlegACx+FlegBx)*I1-(FlegACx+FlegDx)*I2-(FlegACy+FlegDy)*I3+(FlegBy+FlegACy)*I4];
yy = zeros(6,1);
yy(1) = y(4) ;
yy(2) = y(5) ;
yy(3) = y(6) ;
yy(4) = (f(1) - k(1)*y(1) - k(3)*y(6) - .2*k(1)*y(4) - .2*k(3)*y(6) )/m(1);
yy(5) = (f(2) - k(5)*y(2) - k(6)*y(3) - .2*k(5)*y(5) - .2*k(6)*y(6) )/m(5) ;
yy(6) = (f(3) - k(9)*y(6) - k(6)*y(5) - .2*k(6)*y(5) + .2*k(7)*y(4) - .2*k(9)*y(6) )/m(9);
end
  1 comentario
Torsten
Torsten el 17 de Nov. de 2022
And most probably
FlegBx = FlegB * sind(45); %forces at leg B in x direction
FlegBy = FlegB * cosd(45); %forces at leg B in y direction
FlegACx = FlegAC * sind(45); %forces at leg A & C in x direction
FlegACy = FlegAC * cosd(45); %forces at leg A & C in y direction
FlegDx = FlegD * sind(45); %forces at leg D in x direction
FlegDy = FlegD * cosd(45); %forces at leg D in y direction
instead of
FlegBx = FlegB * sin(45); %forces at leg B in x direction
FlegBy = FlegB * cos(45); %forces at leg B in y direction
FlegACx = FlegAC * sin(45); %forces at leg A & C in x direction
FlegACy = FlegAC * cos(45); %forces at leg A & C in y direction
FlegDx = FlegD * sin(45); %forces at leg D in x direction
FlegDy = FlegD * cos(45); %forces at leg D in y direction

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by