First Order Differential Equation

I need to solve the below first ODE with a set time step say (Delt = 0.01s) and from 0 - 20s. Initial condiitons Y0 = [ 0 ; 0 ; 0 ]

dof = 3; %Three Story 
I = 1.5796*10^(-2); %m^4
E = 0.209*10^9; %kN/m^2
m = 500000; %kg
k = 849116255; %N/m


m1 = 2*m;
m2 = 2*m;
m3 = m;

k1 = 3*k;
k2 = 3*k;
k3 = 2*k;

M = [ m1 0 0 ; 0 m2 0 ; 0 0 m3 ]; %Mass Matrix
K = [(k1+k2) -k 0 ; -k2 (k2+k3) -k3 ; 0 -k3 k3 ]; %Stiffness Matrix
C = [ 7 -3 0 ; -3 3.2 -1.4 ; 0 -1.4 1.4 ]*10^5; %Damping Matrix
o = [ 0 0 0 ; 0 0 0 ; 0 0 0 ]; %Zero Matrix

A = [ M o ; o eye(dof) ] 
B = [ o K ; -eye(dof) o ]
P = 2*sin(4*(4*pi*t))
ft = [P ; 0 ; 0]

Yt+1 = Yt+ delt(inv(A)*(ft-B*Yt))

Once this has been computed, how does one store the results for each itteration of time and subsequently plot it?

[Alan Stevens]

Your matrices don't seem to be consistent. e.g. B is 6x6, so it has 6 columns, but you multiply it by Yt, which, if the initial value is to be believed, has just 3 rows.

What is/are the actual ODE's you are trying to solve?

[Christopher Wijnberg]

Hi Alan

My apologies -

ft = [ P ; 0 ; 0 ; 0 ; 0 ; 0] 

I no longer have the original ODE's as they have been fransferred from second order set of n differential equations to a 2n system of first order differential equations by making use of the State Space Form.

Answer 1

[Alan Stevens]

Well, the following shows how to progress through time. However, the explicit form you've adopted is unstable. I've included an implicit form as an alternative. Also, I notice that you haven't included the damping in the equations yet.

dof = 3; %Three Story 
I = 1.5796*10^(-2); %m^4
E = 0.209*10^9; %kN/m^2
m = 500000; %kg
k = 849116255; %N/m

m1 = 2*m;
m2 = 2*m;
m3 = m;

k1 = 3*k;
k2 = 3*k;
k3 = 2*k;

M = [ m1 0 0 ; 0 m2 0 ; 0 0 m3 ]; %Mass Matrix
K = [(k1+k2) -k 0 ; -k2 (k2+k3) -k3 ; 0 -k3 k3 ]; %Stiffness Matrix
C = [ 7 -3 0 ; -3 3.2 -1.4 ; 0 -1.4 1.4 ]*10^5; %Damping Matrix
o = [ 0 0 0 ; 0 0 0 ; 0 0 0 ]; %Zero Matrix

A = [ M o ; o eye(dof) ]; 
B = [ o K ; -eye(dof) o ];
u = ones(6,1);

delt = 0.01;
Tend = 2;
t = 0:delt:Tend;
Y = zeros(6,length(t));

% Explicit - unstable
% for i = 1:length(t)-1
% 
%         P = 2*sin(4*(4*pi*t(i)));
%         ft = [P ; 0 ; 0; 0; 0; 0];
% 
%         Y(:,i+1) = Y(:,i)+ delt*A\(ft-B*Y(:,i));
% 
% end

% Implicit - stable  
for i = 1:length(t)-1

        P = 2*sin(4*(4*pi*t(i)));
        ft = [P ; 0 ; 0; 0; 0; 0];      
        Y(:,i+1) = (Y(i) + delt*A\ft)./(u + delt*A\(B*u)); 
end

plot(t,Y(1,:))

[Alan Stevens]

Note: I've also replaced inv(A)*... by A\..., which Matlab tells me is faster and more accurate!

[Christopher Wijnberg]

Thanks Alan

I am working through this now - Just wondering what is the matrix u and why did you introduce it?

u = ones(6,1);

Further more - how would I account for initial displacement conditions? say -

Y0 = [0 ; 0 ; 0 ; 0 ; 0 ; 0.005]

Where the above matrix describes an initial displacement of 0.005m on the third floor of a three story (3dof) building?

[Christopher Wijnberg]

Note - Y0 is comprised of a 3x1 velocity vector Xdot and a 3x1 displacement vector X as described below.

Xdot = [0 ; 0 ; 0]
X = [ 0 ; 0 ; 0.005]

Y0 = [Xdot ; X]

[Alan Stevens]

I introduced u in order to make the denominator of the implicit equation work (A\B on its own produces a 6x6 matrix, but you want a 6x1 vector).

For your initial displacement vector:

Y = zeros(6,length(t));
Y0(:,1) = [0 ; 0 ; 0 ; 0 ; 0 ; 0.005];

[Christopher Wijnberg]

When I intruduce

Y = zeros(6,length(t));
Y0(:,1) = [0 ; 0 ; 0 ; 0 ; 0 ; 0.005];

and then set P = 0 (Dynamic machine is turned off and I am trying to only analyse the response of the structure to the initial displacement) - I get a plot of zero's for the various times.

dof = 3; %Three Story 
I = 1.5796*10^(-2); %m^4
E = 0.209*10^9; %kN/m^2
m = 500000; %kg
k = 849116255; %N/m


m1 = 2*m;
m2 = 2*m;
m3 = m;

k1 = 3*k;
k2 = 3*k;
k3 = 2*k;

M = [ m1 0 0 ; 0 m2 0 ; 0 0 m3 ]; %Mass Matrix
K = [(k1+k2) -k 0 ; -k2 (k2+k3) -k3 ; 0 -k3 k3 ]; %Stiffness Matrix
C = [ 7 -3 0 ; -3 3.2 -1.4 ; 0 -1.4 1.4 ]*10^5; %Damping Matrix
o = [ 0 0 0 ; 0 0 0 ; 0 0 0 ]; %Zero Matrix
A = [ M zeros(dof , dof) ; zeros(dof , dof) eye(dof) ] 
B = [ zeros(dof , dof) K ; -eye(dof) zeros(dof , dof) ]
u = ones(6,1);
% X   = [0 ; 0 ; 0 ]
% dX  = [0 ; 0 ; 0 ]
% ddX = [0 ; 0 ; 0 ]

delt = 0.01;
Tend = 2;
t = 0:delt:Tend;
Y = zeros(6,length(t));
Y0(:,1) = [0 ; 0 ; 0 ; 0 ; 0 ; 0.005];
  
for i = 1:length(t)-1
        P = 0
        %P = 2*sin(4*(4*pi*t(i)));
        ft = [P ; 0 ; 0; 0; 0; 0];      
        Y(:,i+1) = (Y(i) + delt*A\ft)./(u + delt*A\(B*u)); 
end
plot(t,Y(1,:))

Am I doing something wrong?

[Alan Stevens]

What does plot(t, Y(6,,:) look like? I'm busy for the moment I'll take a look myself later.

[Christopher Wijnberg]

Its a flat line with all values of Y being zero.

[Alan Stevens]

There is no Y0. Try replacing

Y0(:,1) = [0 ; 0 ; 0 ; 0 ; 0 ; 0.005];

by either

Y(:,1) = [0 ; 0 ; 0 ; 0 ; 0 ; 0.005];

or

Y(6,1) = 0.005;

[Christopher Wijnberg]

Sorry Alan

I am just not getting the right outputs, while I should be getting an output for Y(t) at each floor/ dof I am getting a single plot... I am trying not be be useless here or waste your time but please forgive me this is my first time coding in any language let alone in Matlab..

This was a graph taken from a Youtube video for a similar example coded in Python.

Answer 2

[Alan Stevens]

Try plot(t, Y(3:6,:))

[Christopher Wijnberg]

Thank you Alan, that works a treat!

And thanks for ALL your help prior!

Regards