First Order Differential Equation

7 views (last 30 days)
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?
  2 Comments
Alan Stevens
Alan Stevens on 4 Aug 2020
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
Christopher Wijnberg on 4 Aug 2020
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.

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 5 Aug 2020
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,:))
  9 Comments
Alan Stevens
Alan Stevens on 5 Aug 2020
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
Christopher Wijnberg on 5 Aug 2020
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.

Sign in to comment.

More Answers (1)

Alan Stevens
Alan Stevens on 5 Aug 2020
Try plot(t, Y(3:6,:))
  1 Comment
Christopher Wijnberg
Christopher Wijnberg on 5 Aug 2020
Thank you Alan, that works a treat!
And thanks for ALL your help prior!
Regards

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!