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

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;
end

Answer 1

[David Goodmanson]

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.