First Order Differential Equation
조회 수: 3 (최근 30일)
이전 댓글 표시
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?
채택된 답변
Alan Stevens
2020년 8월 5일
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
Alan Stevens
2020년 8월 5일
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;
추가 답변 (1개)
참고 항목
제품
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!