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

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?

[Stephan]

the result of dxdt(1) is a scalar. the result of dxdt(2) is a 2x2 matrix. is that the desired dimension of the result?

[Piotr Haciuk]

the result should be matrix 2 by n, where n are the time steps ode45 is making, than i want to plot first column with respect to the time.

[David Goodmanson]

Hi Piotr,

As Stephan pointed out, dxdt(1) is a scalar as it should be, but you are defining dxdt(2), which should also be a scalar, as a 2x2 matrix. Hence the error. Since dxdt(2) = x'', you should be able to define it as a scalar function of t, x(1), and x(2), similarly to what is in the ode45 documentation for the van der Pol equation. The only difference is that your function explicitly involves t, whereas the vdP equation does not.

[Piotr Haciuk]

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;

[David Goodmanson]

It was not so clear before that you were looking for 2 degrees of freedom, in which case the resulting matrix is going to be nx4. For 2 DOF and a second order DE, there is a four dimensional column vector. Arbitrarily calling the two dimensional coordinates x1 and x2, then in the usual way of doing things the input and output of ode45 are the column vectors

x = [x1; x2; x1'; x2'] and xdot = [[x1'; x2'; x1''; x2'']

The defining function for ode45 puts

xdot(1) = x(3) and xdot(2) = x(4) which is to say xdot(1:2) = x(3:4)

(Since x and xdot are definded as column vectors, x(3:4) is a column vector etc.) Now you need to define the second derivatives, i.e. column vector xdot(3:4) in terms of x, dx/dt and t, that is in terms of x(1:2), x(3:4) and t. It's quite common at this point to use 2x2 matrices for that task.

After that there will be four initial conditions to specify.

[Piotr Haciuk]

That explains a lot, Thank you!

So I rewrote the script as follows:

Function:

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

Main:

clc
clear all
close all
x10=0;
x1dot0=0;
x20=0;
x2dot0=0;
IC=[x10; x1dot0; x20; x2dot0];
t=0:0.1:10;
[t,x]=ode45( @MDOF1, t, IC );

figure(1)
plot(t,x(:,1))
xlabel('x'); ylabel('t');
hold on
plot(t,x(:,2))
hold on
plot(t,x(:,3))
hold on
plot(t,x(:,4))
legend('x1','x2','x3','x4')
hold off

Now my question is if

xdt(3:4)=x(3:4);
dxdt(1:4)=inv(M)*F-(inv(M)*C)*x(3:4)-(inv(M)*K)*x(1:2);
dxdt=[dxdt(1); dxdt(2);dxdt(3); dxdt(4)];

is correct as well as:

IC=[x10; x1dot0; x20; x2dot0];

[Piotr Haciuk]

And the result is in the jpg file above and I'm not sure if this is correct :-(

[David Goodmanson]

after a quick look, here are some observations about the MDOF1 function.

If you predefine dxdt as a column vector which you can do at the beginning with the statement dxdt = zeros(4,1), then you won't need the statement at the end,

dxdt=[dxdt(1); dxdt(2);dxdt(3); dxdt(4)];

The first line after F = ... should be dxdt(1:2)=x(3:4); as I mentioned

and the next line needs to be an equation with a 2x1 column vector on each side:

dxdt(3:4) = something

I haven't seen the differential equation so I can't comment further except that

-(inv(M)*C)*x(3:4)-(inv(M)*K)*x(1:2)

is a column vector but inv(M)*F is a 2x2 matrix. So you get an error.

[Piotr Haciuk]

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]

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]

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.

[Piotr Haciuk]

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]

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.