First Order Differential Equation

Question

Christopher Wijnberg 2020-8-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/575008-first-order-differential-equation

评论： Christopher Wijnberg 2020-8-5

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 个评论
显示无隐藏无

Alan Stevens 2020-8-4

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 2020-8-4

在 MATLAB Online 中打开

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 2020-8-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/575008-first-order-differential-equation#answer_475300

在 MATLAB Online 中打开

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 个评论
显示 7更早的评论隐藏 7更早的评论

Christopher Wijnberg 2020-8-5

在 MATLAB Online 中打开

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 2020-8-5

在 MATLAB Online 中打开

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 2020-8-5

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..