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

1 vue (au cours des 30 derniers jours)
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 commentaires
josh johnson
josh johnson le 18 Nov 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 le 18 Nov 2022
Make a loop over k, call ODE45 in the loop for each value of k separately and save the results.

Connectez-vous pour commenter.

Réponse acceptée

Cris LaPierre
Cris LaPierre le 17 Nov 2022
Modifié(e) : Cris LaPierre le 17 Nov 2022
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 commentaire
Torsten
Torsten le 17 Nov 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

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by