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 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaire
James Tursa
le 9 Août 2022
Modifié(e) : James Tursa
le 9 Août 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?
Réponses (1)
Chrissy Kimemwenze
le 9 Août 2022
2 commentaires
James Tursa
le 10 Août 2022
I don't see the equation, and I don't see your code that tries to use ode45( ).
Voir également
Catégories
En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!