I am trying find the response of a three story building using 4th order runge kutta but my response is diverging I am not sure why can someone help?
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
m = (500*1000)
M = [2*m 0 0; 0 2*m 0; 0 0 m]
[nd,nd]= size(M)
disp('mass matrix')
M
%insert elcentro data
t = e002
accn = e003
dt = t(2) - t(1)
plot(t,accn,'black'); set(gca,'Fontsize',13)
xlabel('TIME (SEC)')
ylabel('ACCELERATION (m/s2)')
%force vector calculations
for i = 1:nd
f(:,i)=-accn*M(i,i)
end
disp ('Force vector')
f
plot(t,f(:,1),'black'); set(gca,'Fontsize',13)
xlabel('TIME (SEC)')
ylabel('FORCE FLOOR 1 (N)')
plot(t,f(:,2),'black'); set(gca,'Fontsize',13)
xlabel('TIME (SEC)')
ylabel('FORCE FLOOR 2 (N)')
plot(t,f(:,3),'black'); set(gca,'Fontsize',13)
xlabel('TIME (SEC)')
ylabel('FORCE FLOOR 3 (N)')
%Damping Matrix input
C = 10^5 * [7 -3 0; -3 3.2 -1.4; 0 -1.4 1.4]
disp ('Damping matrix')
C
%Stiffness Matrix input
E = 0.209*10^9; I =1.5796*10^-2
K1 = 2*E*I
K2 = 3*E*I
K3 = 3*E*I
K =[K1+K2 -K2 0; -K2 K2+K3 -K3; 0 -K3 C(3,3)]
disp ('Stiffness matrix')
K
u1 = [0;0;0]
v1 = [0;0;0]
tt = 53.74
n= 150
n1 = n + 1
F = [0;0;0]
an1 = inv(M) * (F - C * v1 - K * u1)
t=0.0
for i = 2 : n1
F = [f(i,1);f(i,2);f(i,3)]
% for k1
ui = u1
vi = v1
ai = an1
d1 = vi
q1 = ai
% for k2
l = 0.5
x = ui + l * dt * d1
d2 = vi + l * dt * q1
q2 = inv(M) * (F - C * d2 - K * x)
% for k3
l = 0.5
x = ui + l * dt * d2
d3 = vi + l * dt * q2
q3 = inv(M) * (F - C * d3 - K * x)
% for k4
l = 1
x = ui + l * dt * d3
d4 = vi + l * dt * q3
q4 = inv(M) * (F - C * d4 - K * x)
x1 = u1 + dt * (d1 + 2 * d2 + 2 * d3 + d4)/6
ve1 = v1 + dt * (q1 + 2 * q2 + 2 * q3 + q4)/6
anc1 = inv(M) * (F - C * ve1 - K * x1)
u1 = x1
v1 = ve1
an1 = anc1
end
1 comentario
James Tursa
el 9 de Ag. de 2022
Editada: James Tursa
el 9 de Ag. de 2022
This code is hard to read. I suggest you completely rewrite this using a matrix-vector formulation along the lines of the examples you can find in the ode45( ) doc. That is, create a derivative function with a (t,y) signature, where y is the state vector. Then write your RK4 code using this function. That way you can compare your results directly with ode45( ) results using the exact same derivative function.
Also, can you post an image of the differential equation you are solving?
Respuestas (1)
Chrissy Kimemwenze
el 9 de Ag. de 2022
2 comentarios
James Tursa
el 10 de Ag. de 2022
I don't see the equation, and I don't see your code that tries to use ode45( ).
Ver también
Categorías
Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!