MATLAB: First Order Differential Equation

MATLABodeode45plotting

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?

Best Answer

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,:))
Related Question