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

Question

Elyar Ghaffarian Dallali 2025-5-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2177097-how-to-solve-2dof-linear-time-variant-differential-equations-in-matlab-with-non-diagonal-mass-stiff

编辑： William Rose 2025-5-15

SMD6.m

在 MATLAB Online 中打开

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');

5 个评论
显示 3更早的评论隐藏 3更早的评论

William Rose 2025-5-12

@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 2025-5-15

@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 2025-5-12

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2177097-how-to-solve-2dof-linear-time-variant-differential-equations-in-matlab-with-non-diagonal-mass-stiff#answer_1565082

编辑：Torsten 2025-5-12

在 MATLAB Online 中打开

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Elyar Ghaffarian Dallali 2025-5-12

编辑：Torsten 2025-5-12

在 MATLAB Online 中打开

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 2025-5-12

编辑：Torsten 2025-5-12

在 MATLAB Online 中打开

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 2025-5-14

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 2025-5-14

编辑：Torsten 2025-5-14

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 2025-5-12

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2177097-how-to-solve-2dof-linear-time-variant-differential-equations-in-matlab-with-non-diagonal-mass-stiff#answer_1565087

@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

.

5 个评论
显示 3更早的评论隐藏 3更早的评论

William Rose 2025-5-14

beamDynamicLoad.m

@Elyar Ghaffarian Dallali,

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

William Rose 2025-5-15

编辑：William Rose 2025-5-15

beamDynamicLoad2.m

@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.

请先登录，再进行评论。

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

回答（2 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

5 个评论 显示 3更早的评论隐藏 3更早的评论

回答（2 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

5 个评论 显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论