Problem with damping constant affecting displacement in a mass system displacement system

3 views (last 30 days)
When no damping constant is present in my code it works as intended however when a damping constant is added the displacement increases which I believe shouldn't happen for my scenario. I cannot spot where the error is in my code.
% Define parameters
m = 10; % Mass (kg)
k = 160; % Spring constant (N/m)
b = 2; % Damping coefficient (N s/m)
F0 = 200; % Amplitude of forcing term (N)
v0 = [2, 3, 4, 5]; % Frequency of forcing term (rad/s)
x0 = -3; % Initial displacement (m)
v0_dot = 0; % Initial velocity (m/s)
% Define time span
tspan = linspace(0, 60, 1000); % from 0 to 20 seconds
% Define the function for ODE solver
for i = 1:length(v0)
omega = v0(i);
f = @(t, Y) [Y(2); (F0 * sin(omega * t) - k * Y(1) - b * Y(2)) / m];
% Solve using ODE45
[t_ode45{i}, Y_ode45{i}] = ode45(f, tspan, [x0, v0_dot]);
% Solve symbolically
syms x(t)
ode = diff(x, t, 2) == (F0 * sin(omega * t) - k * x - b * diff(x, t)) / m;
cond = [x(0) == x0, subs(diff(x), 0) == v0_dot]; % Corrected indexing for symbolic differentiation
x_sol = dsolve(ode, cond);
x_symbolic{i} = matlabFunction(x_sol);
% Plot results
subplot(length(v0), 1, i);
plot(t_ode45{i}, Y_ode45{i}(:, 1), 'b-', 'LineWidth', 1.5);
hold on;
t_symbolic = linspace(0, 60, 1000);
plot(t_symbolic, x_symbolic{i}(t_symbolic), 'r--', 'LineWidth', 1.5);
title(['Frequency = ', num2str(omega), ' rad/s']);
xlabel('Time (s)');
ylabel('Displacement (m)');
legend('ODE45', 'Symbolic');
grid on;
hold off;

Accepted Answer

David Goodmanson
David Goodmanson on 9 Feb 2024
Hi Daneil,
I don't see any signs that your calculation is incorrect, although in some cases you may have to run the calculation out further in time so that the solution can settle down. If you define
gamma = b/m
w0 = sqrt(k/m)
then for the underdamped case that you have, the requirement is gamma < 2*w0 and the expected amplitude for large times is
A = (F0/m)/sqrt((w^2-w0^2)^2 + gamma^2*w^2)
which agrees pretty well. With underdamping, the particular solutions fall off like exp((-gamma/2)*t) and a suitable wait time for the solution to down to a constant amplitude (so you get a correct result) is
waitime = 8/gamma
Here k/m = w0 = 4 and one of your frequencies is w=4, so in that case you would expect the solution to be quite large for small b, which is what happens. For fixed w, A decreases as b increases and I didn't see any indications otherwise.

More Answers (0)


Find more on Symbolic Math Toolbox 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!