how to solve 2DOF linear time variant differential equations in matlab with non-diagonal mass, stiffness and damping matrixes

Hello everyone!

i hope u already have started an amazing summer.

it's been a while i am struggling to solve 2DOF linear time variant differential equations in matlab with non-diagonal mass, stiffness and damping matrixes using newmark beta method. The associated differential equation is uploaded to the chat related to the equivalent bean considering the first mode of vibration and a moving spring mass damper model that applies dynamic force of G(t). the picture of the physical and the mathematical equation of the system is attached to the chat. according to the uploaded matlab code, my bridge response (or generally let's say system response) depends on the time decrement of the numerical modeling. i dont know why for each value of time decrement "dt" i get ddifferent results. can someone kindly take a look on the equation and my code and give me some comments. much appreciate it!

clc;
clear;
close all;
% ---------------------- PARAMETERS ----------------------
m_p = 70; ksi_p = 0.03; f_p = 2.33; omega_p = 2*pi*f_p;
k_p = m_p * omega_p^2;
c_p = 2 * ksi_p * m_p * omega_p;
m_b = 867.5; kisi_b = 0.008; f_b = 4.58; omega_b = 2*pi*f_b;
k_b = m_b * omega_b^2;
c_b = 2 * kisi_b * m_b * omega_b;
% Load properties
L = 22.9; n_mode = 1; fw = 140/60;
V = 0.714 * fw + 0.055;
% Time parameters
dt = 0.001; T = 13.28; t = 0:dt:T; n = length(t);
% GRF input
DLF = [0.4, 0.1, 0.042, 0.041, 0.01]; phi_phase = [0, 0, 0, 0, 0];
W = m_p * 9.81;
phi2 = sin(pi * V * t / L);
GRF = (W + W * (DLF(1) .* cos(1*2*pi*fw*t + phi_phase(1)) + ...
                DLF(2) .* cos(2*2*pi*fw*t + phi_phase(2)) + ...
                DLF(3) .* cos(3*2*pi*fw*t + phi_phase(3)) + ...
                DLF(4) .* cos(4*2*pi*fw*t + phi_phase(4)) + ...
                DLF(5) .* cos(5*2*pi*fw*t + phi_phase(5)))) .* phi2;
% Initialize state
x = zeros(2, n); xdot = zeros(2, n); xddot = zeros(2, n);
% Initial matrices at t=0
Phi_t = sin(n_mode * pi * V * t(1) / L);
M = [1, m_p*Phi_t; 0, m_p];
C = [c_b, 0; -c_p*Phi_t, c_p];
K = [k_b, 0; -k_p*Phi_t, k_p];
F0 = [GRF(1)/m_b; 0];
xddot(:,1) = M \ (F0 - C*xdot(:,1) - K*x(:,1));
% Newmark-Beta constants
beta = 1/4; gamma = 0.5;
a0 = 1 / (beta * dt^2);
a1 = gamma / (beta * dt);
a2 = 1 / (beta * dt);
a3 = gamma / beta;
a4 = 1 / (2 * beta);
a5 = dt / 2 * (gamma / beta - 2);
a6 = dt * (1 - gamma / (2 * beta));
maxIter = 20; tol = 1e-6;
for k = 1:n-1
    tk = t(k); Phi_t = sin(n_mode * pi * V * tk / L);
    M = [1, m_p*Phi_t; 0, m_p];
    C = [c_b, 0; -c_p*Phi_t, c_p];
    K = [k_b, 0; -k_p*Phi_t, k_p];
    Fk = [GRF(k)/m_b; 0];
    Keff = a0 * M + a1 * C + K;
    A = a2 * M + a3 * C;
    B = a4 * M + a5 * C;
    dFeff = Fk + A * xdot(:,k) + B * xddot(:,k);
    delta_u = Keff \ dFeff;
    % Iteration variables
    u = x(:,k); du = xdot(:,k); ddu = xddot(:,k);
    for i = 1:maxIter
        u_trial = u + delta_u;
        du_trial = du + (a1 * delta_u - a3 * du + a6 * ddu);
        ddu_trial = ddu + (a0 * delta_u - a2 * du - a4 * ddu);
        R = Fk - (M * ddu_trial + C * du_trial + K * u_trial);
        if norm(R) < tol
            break;
        end
        delta_u_corr = Keff \ R;
        delta_u = delta_u + delta_u_corr;
    end
    % Update states
    x(:,k+1) = u + delta_u;
    xdot(:,k+1) = du + (a1 * delta_u - a3 * du + a6 * ddu);
    xddot(:,k+1) = ddu + (a0 * delta_u - a2 * du - a4 * ddu);
end
% Output
acc1 = xddot(1, :); acc2 = xddot(2, :);
fprintf('Max bridge acc: %.6f m/s²\n', max(abs(acc1)));
fprintf('Max pedestrian acc: %.6f m/s²\n', max(abs(acc2)));
% Plot
figure; plot(t, acc1); title('Bridge Acceleration'); grid on;
% FFT
Fs = 1/dt;
N = length(acc1);
Y = fft(acc1); P = 2 * abs(Y / N); P = P(1:N/2 + 1);
f_vector = Fs * (0:(N/2)) / N;
figure; plot(f_vector, P, 'LineWidth', 2); xlim([0 15]);
title('FFT of Bridge Acceleration');

[William Rose]

@Elyar Ghaffarian Dallali,

The quantity in the figure, and in the equation that is part of the figure, appears to correspond the the variable phi2=sin(pi*V*t/L) in the code. Am I understanding correctly?

[William Rose]

@Elyar Ghaffarian Dallali, the figure and the equation in the figure refer to and to . The code does not include GRS or G_RS, but the dcode does have GRF=... . Am I correct that GRF defined in the code corresponds to in the equation in the figure?

[Elyar Ghaffarian Dallali]

Yes exactly

[William Rose]

@Elyar Ghaffarian Dallali,

I think there is an inconsistency between the equation in the figure and the implementation in your code. The equation in the figure suggests that element C(1,1) in your code should equal . But your code defines c_b = 2*kisi_b*m_b*omega_b and C(1,1)=c_b. In other words, your C(1,1) includes a factor m_b which is not present in the equation in the figure.

A similar problem appears with array K. The equation in the figure suggests that element K(1,1) in your code should equal . But your code defines k_b=m_b*omega_b^2 and K(1,1)=k_b. In other words, your K(1,1) includes a factor m_b which is not present in the equation in the figure.

Analysis of the presumed units also suggests that the code is incorrect:

The forcing term in the code is Fk = [GRF(k)/m_b; 0] (which is consistent with the equation in the figure). If GRF(k) has units of force, then Fk(1) has units of (force/mass)=acceleration. In that case, you do not want a mass term in K(1,1) or C(1,1).

You said to @Torsten: "however, the acceleration response of the bridge DOF is irrational according to the provided code". Do you mean the acceleration has an imaginary component? If so, could it be related to the values of C(1,1) and K(1,1)?

[William Rose]

@Elyar Ghaffarian Dallali,

What is the source of the figure showing the bridge and the coupled equations? I am asking because the equation in the figure looks wrong to me. It appears to have inconsistent units.

I assume and y(t) are displacements (units of length), and is dimensionless. The first part of the equation in the figure is

.

When the top row of the matrix multiplies the column vector, the result is . This represents adding acceleration () to (mass times acceleration) (). But one may not add quantities with different units. Therefore the equation in the figure is wrong, or my assumptions about the units for the variables are wrong.

Perhaps the equation in the figure should be

where and . The units are consistent in this equation. You could use @Torsten's approach or mine to solve this equation.

Answer 1

[Torsten]

x0 = [0;0;0;0];
tspan = [0 13.28];
[T,X] = ode45(@fun,tspan,x0);
for i = 1:numel(T)
  dxdt = fun(T(i),X(i,:));
  q1ddot(i) = dxdt(2);
  yddot(i) = dxdt(4);
end
figure(1)
plot(T,X(:,1))   %q1

figure(2)
plot(T,X(:,3))   %y

figure(3)
plot(T,X(:,2))    %q1dot

●

figure(4)
plot(T,X(:,4))    %ydot

figure(5)
plot(T,q1ddot)  %q1ddot

figure(6)
plot(T,yddot)   % yddot

function dxdt = fun(t,x)
  q1 = x(1);
  q1dot = x(2);
  y = x(3);
  ydot = x(4);
  % ---------------------- PARAMETERS ----------------------
  m_p = 70; ksi_p = 0.03; f_p = 2.33; omega_p = 2*pi*f_p;
  k_p = m_p * omega_p^2;
  c_p = 2 * ksi_p * m_p * omega_p;
  m_b = 867.5; kisi_b = 0.008; f_b = 4.58; omega_b = 2*pi*f_b;
  k_b = m_b * omega_b^2;
  c_b = 2 * kisi_b * m_b * omega_b;
  % Load properties
  L = 22.9; n_mode = 1; fw = 140/60;
  V = 0.714 * fw + 0.055;
  % GRF input
  DLF = [0.4, 0.1, 0.042, 0.041, 0.01]; phi_phase = [0, 0, 0, 0, 0];
  W = m_p * 9.81;
  phi2 = sin(pi * V * t / L);
  GRF = (W + W * (DLF(1) .* cos(1*2*pi*fw*t + phi_phase(1)) + ...
                  DLF(2) .* cos(2*2*pi*fw*t + phi_phase(2)) + ...
                  DLF(3) .* cos(3*2*pi*fw*t + phi_phase(3)) + ...
                  DLF(4) .* cos(4*2*pi*fw*t + phi_phase(4)) + ...
                  DLF(5) .* cos(5*2*pi*fw*t + phi_phase(5)))) .* phi2;
  Phi_t = sin(n_mode * pi * V * t / L);
  M = [1, m_p*Phi_t; 0, m_p];
  C = [c_b, 0; -c_p*Phi_t, c_p];
  K = [k_b, 0; -k_p*Phi_t, k_p];
  Fk = [GRF/m_b; 0];
  acc = M\(Fk - C*[q1dot;ydot] - K*[q1;y]);
  dxdt = [q1dot;acc(1);ydot;acc(2)];
end

[Elyar Ghaffarian Dallali]

Thank u so much @Torsten for the provided code however, the acceleration response of the bridge DOF is irrational according to the provided code. i got another code in which the equation is solved using state space representation of the coupled system and then RK4 numerical modeling is applied to solve it. which gives way better and more reasonable results. what do u think about this code? thanks? %% RK4 Simulation of Time-Varying State-Space System with GRF Input

clc;
clear;
close all;
TV = 13.26;
% === Time Vector ===
T = 0:0.01:TV;
dt = T(2) - T(1);
nSteps = length(T);
% === Pedestrian Properties ===
fp = 140/60;
omegap = fp * 2 * pi;
mp = 70;
kisip = 0.38 * log(mp) - 1.42;
kp = omegap^2 * mp;
cp = 2 * kisip * omegap * mp;
% === Bridge Properties ===
omegab = 4.58 * 2 * pi;
kisib = 0.008;
mb = 1735 / 2;
kb = omegab^2 * mb;
cb = 2 * kisib * omegab * mb;
L = 22.9;
% === Walking Parameters ===
fw = fp;
V = 0.714 * fw + 0.055;
phi = (pi * V / L);
% === GRF Harmonic Parameters ===
DLF = [0.4, 0.1, 0.042, 0.041, 0.01];
phi_phase = [0, 0, 0, 0, 0];
W = mp * 9.81;
% === Initialize State Vector ===
x0 = zeros(4, 1);           % [q1; q2; q1_dot; q2_dot]
X = zeros(4, nSteps);       % Store full response
X(:, 1) = x0;
% === Precompute Time-Varying System Matrices ===
Bmat = cell(1, nSteps);
for t = 1:nSteps
    z = T(t);
    Mmat = [1, mp*sin(phi*z)/mb;
        0, mp];
    Kmat = [omegab^2, 0;
        -kp*sin(phi*z)/mb, kp];
    Cmat = [2*kisib*omegab, 0;
        -cp*sin(phi*z)/mb, cp];
    mk = -Mmat \ Kmat;
    mc = -Mmat \ Cmat;
    Bmat{t} = [zeros(2), eye(2); mk, mc];
end
% === RK4 Integration Loop with GRF Input ===
for k = 1:nSteps-1
    t = T(k);
    B = Bmat{k};
    z = T(k);
    % --- Compute GRF (force applied to bridge DOF only) ---
    phi2 = sin(pi * V * z / L);  % Mode shape at pedestrian location
    harmonic_sum = 0;
    for n = 1:5
        harmonic_sum = harmonic_sum + DLF(n) * cos(n * 2*pi*fw*t + phi_phase(n));
    end
    GRF = (W + W * harmonic_sum) * phi2;
    % --- External force vector: applied to bridge DOF acceleration (3rd row) ---
    F = [0; 0; GRF./mb; 0];
    % --- RK4 integration step ---
    f = @(x) B * x + F;
    k1 = f(X(:, k));
    k2 = f(X(:, k) + 0.5 * dt * k1);
    k3 = f(X(:, k) + 0.5 * dt * k2);
    k4 = f(X(:, k) + dt * k3);
    X(:, k+1) = X(:, k) + (dt/6) * (k1 + 2*k2 + 2*k3 + k4);
end
% === Optional: Plot Bridge Acceleration ===
bridge_acc = diff(X(3, :)) / dt;
pedest_acc = diff(X(4, :)) / dt;
figure;
plot(T(1:end-1), bridge_acc, 'k', 'LineWidth', 1.5);
xlabel('Time (s)', 'FontSize', 18);
ylabel('Bridge Acceleration (m/s^2)', 'FontSize', 18);
title('Bridge Acceleration Response', 'FontSize', 20);
grid on;
set(gca, 'FontSize', 14);

Signal = bridge_acc;
Fs = 100;
% Perform FFT
N = length(Signal);
Y = fft(Signal);
P2 = abs(Y / N);
P1 = 2 * P2(1:N/2 + 1);
f_vector = Fs * (0:(N/2)) / N;
figure;
plot(f_vector, P1);

disp(max(bridge_acc));

    2.3560

[Torsten]

Thank u so much @Torsten for the provided code however, the acceleration response of the bridge DOF is irrational according to the provided code.

Your settings have changed - thus you get different results. That's what I get with your actual parameters and equations using ode45.

Note that in your code from above, you forgot to multiply F by inv(Mmat).

Further, your equations depend explicitly on t - thus you must make your function handle f depend on t and x for application of the Runge-Kutta method.

x0 = [0;0;0;0];
tspan = [0 13.28];
[T,X] = ode45(@fun,tspan,x0);
for i = 1:numel(T)
  dxdt = fun(T(i),X(i,:));
  q1ddot(i) = dxdt(2);
  yddot(i) = dxdt(4);
end
figure(1)
plot(T,X(:,1))   %q1

figure(2)
plot(T,X(:,3))   %y

figure(3)
plot(T,X(:,2))    %q1dot

figure(4)
plot(T,X(:,4))    %ydot

figure(5)
plot(T,q1ddot)  %q1ddot

figure(6)
plot(T,yddot)   % yddot

Signal = q1ddot;
Fs = 100;
% Perform FFT
N = length(Signal);
Y = fft(Signal);
P2 = abs(Y / N);
P1 = 2 * P2(1:N/2 + 1);

Warning: Integer operands are required for colon operator when used as index.

f_vector = Fs * (0:(N/2)) / N;
figure;
plot(f_vector, P1);

disp(max(q1ddot));

    2.3904

function dxdt = fun(t,x)
  q1 = x(1);
  q1dot = x(2);
  y = x(3);
  ydot = x(4);
  % === Pedestrian Properties ===
  fp = 140/60;
  omegap = fp * 2 * pi;
  mp = 70;
  kisip = 0.38 * log(mp) - 1.42;
  kp = omegap^2 * mp;
  cp = 2 * kisip * omegap * mp;
  % === Bridge Properties ===
  omegab = 4.58 * 2 * pi;
  kisib = 0.008;
  mb = 1735 / 2;
  kb = omegab^2 * mb;
  cb = 2 * kisib * omegab * mb;
  L = 22.9;
  % === Walking Parameters ===
  fw = fp;
  V = 0.714 * fw + 0.055;
  phi = (pi * V / L);
  % === GRF Harmonic Parameters ===
  DLF = [0.4, 0.1, 0.042, 0.041, 0.01];
  phi_phase = [0, 0, 0, 0, 0];
  W = mp * 9.81;
  Mmat = [1, mp*sin(phi*t)/mb;
          0, mp];
  Kmat = [omegab^2, 0;
         -kp*sin(phi*t)/mb, kp];
  Cmat = [2*kisib*omegab, 0;
          -cp*sin(phi*t)/mb, cp];
  % --- Compute GRF (force applied to bridge DOF only) ---
  phi2 = sin(pi * V * t / L);  % Mode shape at pedestrian location
  harmonic_sum = 0;
  for n = 1:5
     harmonic_sum = harmonic_sum + DLF(n) * cos(n * 2*pi*fw*t + phi_phase(n));
  end
  GRF = (W + W * harmonic_sum) * phi2;
  % --- External force vector: applied to bridge DOF acceleration (3rd row) ---
  F = [GRF./mb; 0];
  acc = Mmat\(F - Cmat*[q1dot;ydot] - Kmat*[q1;y]);
  dxdt = [q1dot;acc(1);ydot;acc(2)];
end

[Elyar Ghaffarian Dallali]

Thank you so much for providing the code @Torsten, however, as u may have noticed in the main differential equation for the mass, stiffness, and damping we have this form:

but in the code that u have provided:

Mmat = [1, mp*sin(phi*t)/mb;

0, mp];

Kmat = [omegab^2, 0;

-kp*sin(phi*t)/mb, kp];

Cmat = [2*kisib*omegab, 0;

-cp*sin(phi*t)/mb, cp]; which means that M(1,2), C(1,2), and K(1,2), which is the main source of our time variancy in the problem, has been normalized by mb (mass of the bridge) that we can not see the same thing in the main governing differential equation. meanwhile, when we remove this component, the response of the system becomes odd and nonrealistic. what do u think about it?

thank you so much and have a great day!

[Torsten]

I just took the matrices as you defined them in your code and trusted in your expertise. But it shouldn't be a problem to modify them according to your needs, is it ?

Answer 2

[William Rose]

@Elyar Ghaffarian Dallali,

@Torsten has given a great answer with plots. Here is another way of looking at the problem - which shoudl lead to the exact same results as @Torsten's approach.

Your original equation is

(eq.1)

where (eq.2).

Define (eq.3). Then we can write

(eq.4)

Combine eq.3 and eq.4:

(eq.5)

Equation 5 has the form needed for ode45(), as long as you use eq.2 when evaluating the right hand side of eq.5, at every time step. You could define

and then eq.5 could be written as

And that is what @Torsten has done and solved. His vector dxdt is like my , except the elements are in a different order. He has , and I have .

[Elyar Ghaffarian Dallali]

Thank you so much for the mathematical background

[William Rose]

@Elyar Ghaffarian Dallali,

The attached script solves the system with ode45(), like @Torsten's script. My dxdt function specifies the 4x4 matrices I defined in my answer above. I did not use k_b or c_b. Instead, I used and . See my comment and questions above. My solutions look different from @Torsten's solutions.

I certainly could have made mistakes in my calculations or my code. Therefore I recommend that you check all my work carefully.

[Elyar Ghaffarian Dallali]

thank u so much @William Rose but can u please give me the access to the code u ve run?

have a great day

[William Rose]

@Elyar Ghaffarian Dallali,

Sorry I forgot to attach the code earlier. Here is its.

[William Rose]

@Elyar Ghaffarian Dallali,

[Edit: Add plot of acceleration.]

I solved the system using the alternate equation, which I proposed in a comment above.

The equation in the figure you posted indicates the addition of quantities with different units, which is an error. That is why I propose the alternate equation - see comment above for details.. When I solve with the alternate equation, I get the solution shown. Code attached.