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?

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

[James Tursa]

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?

Answer 1

[Chrissy Kimemwenze]

Here is the equation I am trying to solve. I tried using ode45 but I can get it to work.

[James Tursa]

I don't see the equation, and I don't see your code that tries to use ode45( ).

[Chrissy Kimemwenze]

% 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