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?
1 view (last 30 days)
Show older comments
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 Comment
James Tursa
on 9 Aug 2022
Edited: James Tursa
on 9 Aug 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?
Answers (1)
Chrissy Kimemwenze
on 9 Aug 2022
2 Comments
James Tursa
on 10 Aug 2022
I don't see the equation, and I don't see your code that tries to use ode45( ).
See Also
Categories
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!