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)
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
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?

Sign in to comment.

Answers (1)

Chrissy Kimemwenze
Chrissy Kimemwenze on 9 Aug 2022
Here is the equation I am trying to solve. I tried using ode45 but I can get it to work.
  2 Comments
Chrissy Kimemwenze
Chrissy Kimemwenze on 11 Aug 2022
% HERE IS MY ODE45 CODE how can I input the earthquake load since the load is changing at each time step
tspan=[0:0.02:53.76];
z0=[0;0;0.01;0;0.05;0];
[t,z]=ode45(@forced_vibration,tspan,z0);
plot(t,z(:,1),'color','b','LineWidth',2);
grid on
xlabel('Time (t) (sec)')
ylabel('Displacement x1 (m)')
title('Displacement floor 1 Vs Time')
plot(t,z(:,2),'color','b','LineWidth',2)
grid on
xlabel('Time (t) (sec)')
ylabel('Velocity V1 (m/s)')
title('Velocity Vs Time')
plot(t,z(:,3),'color','b','LineWidth',2)
grid on
xlabel('Time (t) (sec)')
ylabel('Displacement x2 (m)');
title('Displacement floor 2 Vs Time')
plot(t,z(:,4),'color','b','LineWidth',2)
grid on
xlabel('Time (t) (sec)')
ylabel('Velocity V2 (m/s)')
title('Velocity Vs Time')
plot(t,z(:,5),'color','b','LineWidth',2)
grid on
xlabel('Time (t) (sec)')
ylabel('Displacement x3 (m)');
title('Displacement floor 3 Vs Time')
plot(t,z(:,6),'color','b','LineWidth',2)
grid on
xlabel('Time (t) (sec)')
ylabel('Velocity V3 (m/s)')
title('Velocity Vs Time')
function f = forced_vibration(t,z)
% Mass(kg)
m = (500*1000);
M = [2*m 0 0; 0 2*m 0; 0 0 m];
[nd,nd]= size(M);
% Damping coefficient (Ns/m)
C = 10^5 * [7 -3 0; -3 3.2 -1.4; 0 -1.4 1.4];
% Stiffness coefficient (N/m)
E = 0.209*10^9; I =1.5796*10^-2;
h = 3.6;
K1 = 36*E*I/h^3;
K2 = 36*E*I/h^3;
K3 = 24*E*I/h^3;
K =[K1+K2 -K2 0; -K2 K2+K3 -K3; 0 -K3 C(3,3)];
m1 = M(1,1);
m2 = M(2,2);
m3 = M(3,3);
c1 = C(1,1);
c2 = -C(2,1);
c3 = C(2,2);
c4 = C(2,3);
k1 = K(1,1);
k2 = -K(2,1);
k3 = K(2,2);
k4 = -K(2,3);
force = [ 0 , 0, 0];
f = [z(2); 1/m1 * (force(1,1) - c1 * z(2) + c2 * z(4) - k1 * z(1) + k2 * z(3)); z(4);1/m2 * (force(1,2) + c2 * z(2) - c3 * z(4) + c4 * z(6) + k2 * z(1) - k3 * z(3) + k4 * z(5)) ; z(6); 1/m3 * (force(1,3) + c4 * z(4) - c4 * z(6) + k4 * z(3) - k4 * z(5))] ;
end

Sign in to comment.

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!