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

Question

Piotr Haciuk 2019-1-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/441268-ode45-2dof-unable-to-perform-assignment-because-the-left-and-right-sides-have-a-different-number-of

编辑： David Goodmanson 2019-1-28

采纳的回答： David Goodmanson

在 MATLAB Online 中打开

Hi good people,

I'm trying to solve this:

clc
clear all
close all
x0=0;
xdot0=0;
IC=[x0 xdot0];
t=0:0.1:10;
[t,x]=ode45( @MDOF1, t, IC );

when the MDOF1 function looks like:

function dxdt=MDOF1(t,x)
m=[2 2];
M=diag(m);
C=[4 -2;-2 4];
K=[8 -4;-4 8];
F=10*sin(10*t);
dxdt(1)=x(2);
dxdt(2)=inv(M)*F-(inv(M)*C)*x(2)-(inv(M)*K)*x(1)
dxdt=[dxdt(1);dxdt(2)];
    end

and i get this error:

    Unable to perform assignment because the left and right sides have a different number of elements.
Error in MDOF1 (line 8)
dxdt(2)=inv(M)*F-(inv(M)*C)*x(2)-(inv(M)*K)*x(1)
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 mdof1_ode (line 8)
[t,x]=ode45( @MDOF1, t, IC );

Could anyone advise me what I'm doing wrong?

10 个评论
显示 8更早的评论隐藏 8更早的评论

Piotr Haciuk 2019-1-25

在 MATLAB Online 中打开

Your comments seem to be common sense, but I still don't know how to rewrite the code :-(

Any suggestions?

Help highly appreciated

P.S.

I've done the same thing for a one DOF problem and it works, but I really don't know how to translate it to 2 DOF, 1 DOF code below:

Function script:

function dx=mass(t,w)
global c1  M1 K1 f1 tm tm1 wn
f2=f1;%*sin(wn*t);
if t<tm ;
    
        f11=0;
        dx(1) = w(2);
        dx(2) = -(c1/M1)*w(2) - (K1/M1)*w(1) + f11/M1;
   elseif t>tm && t<tm1;
       
        f11=f2;
        dx(1) = w(2);
        dx(2) = -(c1/M1)*w(2) - (K1/M1)*w(1) + f11/M1;
   else t>tm1;
       
        f11=0;
        dx(1) = w(2);
        dx(2) = -(c1/M1)*w(2) - (K1/M1)*w(1) + f11/M1;
end
        dx=[dx(1); dx(2)];
    end

Main script:

clc 
clear all
global c1  M1 K1 f1 tm tm1 wn
%%%%%%%%%%%%%%%%%%%%%%%% BEAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%properties
%beam (concrete)
W=5;                                         % width[m]
L=200;                                       % length[m]
H=2;                                         % thickness[m]
rho=2400;                                    % density [kg/m3]
E=41e+9;                                     % Young's modulus[Pa]
A=W*H;                                       % cross section area [m^2]
K1=(E*A)/L;                                  % stiffness [N/m]
V=W*L*H;                                     % volumne [m^3]
M1=rho*V;                                    % mass [kg]
wn=sqrt(K1/M1);                              % natural frequency [rad/s]
cr1=2*sqrt(K1*M1);                           % critical damping [Ns/m]
%moving mass (avg car)
m=1500;                                        % [kg]
g=9.81;                                      % [m/s^2]
f1=m*g*(-1);
%%%%%%%%%%%%%%%%%%%%%%%% Initial conditions (t=0) %%%%%%%%%%%%%%%%%%%%%%%%%
t0=0;                                        % time [s]
t=0:0.1:10;                                  % time scale [s]    
%beam
w0=0;                                        % displacement [m]
wdot0=0;                                     % velocity [m/s]
IC1=[w0 wdot0];                  
%moving mass
xm0=0;                                       % position of the mass
% tm=5;                                        % time when mass on beam [s]
% tm1=tm+1;                                    % time when mass is gone [s]
%%%%%%%%%%%%%%%%%%%%%%%% Calculations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c1=190000000;                                % damping [Ns/m]
zeta=c1/cr1                                  % damping coeficient
tm=0;
c1=0;
for  c1 = 0: cr1/20 : cr1+0.1*cr1;
    tm=tm+1;
    tm1=tm+0.5;
    zeta=c1/cr1;
[t,w]=ode45( @mass, t, IC1 );                % solver
%%%%%%%%%%%%%%%%%%%%%%%% Results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
plot(t,w(:,1));
xlabel('t [s]'); ylabel('w [m]');
legend('displacement')%, 'velocity')
title(sprintf('Critical damping=%4.1f, damping =%4.1f, zeta =%4.3f, time=%4.1f ',cr1,c1,zeta,tm));
% ylim([-2 2]);
    drawnow;

Piotr Haciuk 2019-1-26

在 MATLAB Online 中打开

David,

I made a slight modification to the function script:

function dxdt=MDOF1(t,x)
m=[10000 10000];
M=diag(m);
C=[800 -400;-400 800];
K=[800 -400;-400 800];
F1=100*sin(10*t);
F=[F1; F1];
dxdt = zeros(4,1);
dxdt(1:2)=x(3:4);
dxdt(3:4)=inv(M)*F-(inv(M)*C)*x(3:4)-(inv(M)*K)*x(1:2);
end

The equation I'm aiming to solve looks like:

[M]x''+[C]x'+[K]x=F

Dividing through by [M] and leaving x'' on the LHS gives:

x''=F/[M]-[C]x'-[K]x

I assume the function script is written correctly now

Next question is how to extract x(1),x(2) and x'(1),x'(2) and plot them as a total response?

What i have so far is this:

clc
clear all
close all
x10=0;
x1dot0=0;
x20=0;
x2dot0=0;
IC=[x10 x20 x1dot0  x2dot0];
t=0:0.01:500;
[t,x]=ode45( @MDOF1, t, IC );
x1 = x(:,1);
x2 = x(:,2);
v1 = x(:,3);%x'(1)
v2 = x(:,4);%x'(2)
xfinal=x1+v1;
vfinal=x2+v2;

Piotr Haciuk 2019-1-26

在 MATLAB Online 中打开

untitled.jpg

Guess I got it!

function dxdt=MDOF1(t,x)
m=[100000 100000];
M=diag(m);
C=[80000 -40000;-40000 80000];
K=[80000 -40000;-40000 80000];
F1=100*sin(10*t);
F2=10*sin(5*t);
F=[F1; F2];
dxdt = zeros(4,1);
dxdt(1:2)=x(3:4);
dxdt(3:4)=inv(M)*F-(inv(M)*C)*x(3:4)-(inv(M)*K)*x(1:2);
end

Main:

clc
clear all
close all
x10=1;
x1dot0=0.1;
x20=0.5;
x2dot0=0.05;
IC=[x10 x20 x1dot0   x2dot0];
t=0:0.01:50;
[t,x]=ode45( @MDOF1, t, IC );
plot(t,x(:,1))
xlabel('t'); ylabel('x');
hold on
plot(t,x(:,2))
hold on
plot(t,x(:,3))
hold on
plot(t,x(:,4))
legend('x1','x2','v1','v2')
hold off

And the result attached

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Goodmanson 2019-1-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/441268-ode45-2dof-unable-to-perform-assignment-because-the-left-and-right-sides-have-a-different-number-of#answer_358165

HI Piotr,

I decided to post one last observation as an answer, which you may respond to as you wish.

Since you have a couple of driving sine functions, the system must eventually settle down to an equilibrium state of constant amplitudes for x's and v's. My first reaction on seeing your plot was that there must be something wrong, because the solution was decaying away. However, it turns out that your initial conditions are quite large compared to equilibrium. If you run the time out to 100 seconds and zoom in, you can see equilibrium ampliudes around 10^-4, considerably smaller than your initial conditions. Or if you use IC = [0 0 0 0] you can see the system settle down quite a bit more quickly. Overall it appears that things are looking good.

2 个评论
显示无隐藏无

Piotr Haciuk 2019-1-27

I'd like to thank you for your contribution towards the outcome, without your help I'd struggle for the next few weeks. so thank you very much! now i have to built on that my main target. Respond to the masses movements on the surface of a bridge modelled using Finite element method.

David Goodmanson 2019-1-28

编辑：David Goodmanson 2019-1-28

Your'e welcome, I was glad to contribute. At this point expanding to more elements is simple, at least in principle, and hopefully you will get good results.

请先登录，再进行评论。

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

10 个评论
显示 8更早的评论隐藏 8更早的评论

采纳的回答

2 个评论
显示无隐藏无

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

10 个评论 显示 8更早的评论隐藏 8更早的评论

采纳的回答

2 个评论 显示 无隐藏 无

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

10 个评论
显示 8更早的评论隐藏 8更早的评论

2 个评论
显示无隐藏无