Plotting Damped Nonlinear Oscillator

The given system is a damped nonlinear oscillator with a nonlinear resistance, described by the differential equation:

where .

I tried to simulate the above problem for different values of .

As you can see below I have results for , but the code it gives me nothing for different values of . What should I change?

• 

% Define common parameters
v = 1; % Damping coefficient
a = 1; % Linear stiffness coefficient
T = 100; % Total simulation time
dt = 0.01; % Time step
t = 0:dt:T; % Time vector
b = (1/20) * a; % Nonlinear stiffness coefficient
x0 = 0.1; % Initial displacement
y0 = 0; % Initial velocity
m_values = [0.1, 1, 10]; % Different mass values

% Initialize a figure
figure(1);
hold on

% Loop through each mass value
for m = m_values
    N = length(t); % Number of time steps

    % Initialize arrays for displacement and velocity
    x = zeros(1, N);
    y = zeros(1, N);

    % Set initial conditions
    x(1) = x0;
    y(1) = y0;

    % Simulate the oscillator's response
    for i = 2:N
        % Calculate the acceleration and velocity changes using equations of motion
        f = -(v/m) * y(i-1) - (a/m) * x(i-1) + (b/m) * x(i-1)^3;
        g = y(i-1);

        % Update velocity and displacement using Euler's method
        y(i) = y(i-1) + dt * f;
        x(i) = x(i-1) + dt * g;
    end

    % Plot the displacement and velocity for the current mass value
    plot(t, x, 'DisplayName', ['x, m = ', num2str(m)]);
    plot(t, y, '--', 'DisplayName', ['y, m = ', num2str(m)]);
end

% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Amplitude');
title('Damped Nonlinear Oscillator Response for Different Mass Values');
legend('show');
grid on;
hold off;

• 

% Define common parameters
v = 1; % Damping coefficient
m = 1; % Mass
T = 50; % Total simulation time
dt = 0.01; % Time step
t = 0:dt:T; % Time vector
x0 = 0.1; % Initial displacement
y0 = 0; % Initial velocity
a_values = [0.1, 1, 10]; % Different stiffness values

% Initialize a figure
figure(2);
hold on

% Loop through each stiffness value
for a = a_values
    b = (1/20) * a; % Nonlinear stiffness coefficient dependent on a
    N = length(t); % Number of time steps

    % Initialize arrays for displacement and velocity
    x = zeros(1, N);
    y = zeros(1, N);

    % Set initial conditions
    x(1) = x0;
    y(1) = y0;

    % Simulate the oscillator's response
    for i = 2:N
        % Calculate the acceleration and velocity changes using equations of motion
        f = -(v/m) * y(i-1) - (a/m) * x(i-1) + (b/m) * x(i-1)^3;
        g = y(i-1);

        % Update velocity and displacement using Euler's method
        y(i) = y(i-1) + dt * f;
        x(i) = x(i-1) + dt * g;
    end

    % Plot the displacement and velocity for the current stiffness value
    plot(t, x, 'DisplayName', ['x, a = ', num2str(a)]);
    plot(t, y, '--', 'DisplayName', ['y, a = ', num2str(a)]);
end

% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Amplitude');
title('Damped Nonlinear Oscillator Response for Different Stiffness Values');
legend('show');
grid on;
hold off;

• 

  % Define common parameters
m = 1; % Mass, kept constant
a = 1; % Linear stiffness coefficient
T = 100; % Total simulation time
dt = 0.01; % Time step
t = 0:dt:T; % Time vector
b = (1/20) * a; % Nonlinear stiffness coefficient
x0 = 0.1; % Initial displacement
y0 = 0; % Initial velocity
v_values = [0.1, 1, 10]; % Different damping coefficients

% Initialize a figure
figure(3);
hold on

% Loop through each damping coefficient value
for v = v_values
    N = length(t); % Number of time steps

    % Initialize arrays for displacement and velocity
    x = zeros(1, N);
    y = zeros(1, N);

    % Set initial conditions
    x(1) = x0;
    y(1) = y0;

    % Simulate the oscillator's response
    for i = 2:N
        % Calculate the acceleration and velocity changes using equations of motion
        f = -(v/m) * y(i-1) - (a/m) * x(i-1) + (b/m) * x(i-1)^3;
        g = y(i-1);

        % Update velocity and displacement using Euler's method
        y(i) = y(i-1) + dt * f;
        x(i) = x(i-1) + dt * g;
    end

    % Plot the displacement and velocity for the current damping value
    plot(t, x, 'DisplayName', ['Displacement, v = ', num2str(v)]);
    plot(t, y, '--', 'DisplayName', ['Velocity, v = ', num2str(v)]);
end

% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Amplitude');
title('Damped Nonlinear Oscillator Response for Different Damping Values');
legend('show');
grid on;
hold off;

• 

% Define common parameters
v = 1; % Damping coefficient
a = 1; % Linear stiffness coefficient
T = 10; % Total simulation time
dt = 0.01; % Time step
t = 0:dt:T; % Time vector
b = (1/20) * a; % Nonlinear stiffness coefficient
m = 1; % Mass (constant in this case)
x0_values = [0.1, 1, 10]; % Different initial displacements
y0 = 0; % Initial velocity

% Initialize a figure
figure(4);
hold on

% Loop through each initial displacement value
for x0 = x0_values
    N = length(t); % Number of time steps

    % Initialize arrays for displacement and velocity
    x = zeros(1, N);
    y = zeros(1, N);

    % Set initial conditions
    x(1) = x0;
    y(1) = y0;

    % Simulate the oscillator's response
    for i = 2:N
        % Calculate the acceleration and velocity changes using equations of motion
        f = -(v/m) * y(i-1) - (a/m) * x(i-1) + (b/m) * x(i-1)^3;
        g = y(i-1);

        % Update velocity and displacement using Euler's method
        y(i) = y(i-1) + dt * f;
        x(i) = x(i-1) + dt * g;
    end

    % Plot the displacement and velocity for the current initial displacement value
    plot(t, x, 'DisplayName', ['Displacement, x_0 = ', num2str(x0)]);
    plot(t, y, '--', 'DisplayName', ['Velocity, x_0 = ', num2str(x0)]);
end

% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Amplitude');
%title('Damped Nonlinear Oscillator Response for Different Initial Displacements');
legend('show');
grid on;
hold off;

[Torsten]

I removed the title of the last graph. Do you see the y-scale value of 1e252 ?

[Athanasios Paraskevopoulos]

yes! Should I take smaller values of x_0?

[Torsten]

Seems x0 = 10 is too large.

% Define common parameters
v = 1; % Damping coefficient
a = 1; % Linear stiffness coefficient
T = 10; % Total simulation time
dt = 0.1; % Time step
t = 0:dt:T; % Time vector
b = (1/20) * a; % Nonlinear stiffness coefficient
m = 1; % Mass (constant in this case)
x0_values = [0.1,1,10]; % Different initial displacements
y0 = 0; % Initial velocity
% Initialize a figure
figure(4);
hold on
% Loop through each initial displacement value
for x0 = x0_values
    [t,z] = ode45(@(t,z)[z(2);-(v/m) * z(2) - (a/m) * z(1) + (b/m) * z(1)^3],t,[x0;y0]);
    x = z(:,1);
    y = z(:,2);
    % Plot the displacement and velocity for the current initial displacement value
    plot(t, x, 'DisplayName', ['Displacement, x_0 = ', num2str(x0)]);
    plot(t, y, '--', 'DisplayName', ['Velocity, x_0 = ', num2str(x0)]);
end

Warning: Failure at t=9.360545e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t.

% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Amplitude');
title('Damped Nonlinear Oscillator Response for Different Initial Displacements');
legend('show');
grid on;
hold off;

[Athanasios Paraskevopoulos]

You are right for the code works. Thank you!

% Define common parameters
v = 1; % Damping coefficient
a = 1; % Linear stiffness coefficient
T = 10; % Total simulation time
dt = 0.1; % Time step
t = 0:dt:T; % Time vector
b = (1/20) * a; % Nonlinear stiffness coefficient
m = 1; % Mass (constant in this case)
x0_values = [0.1,0.5,1]; % Different initial displacements
y0 = 0; % Initial velocity
% Initialize a figure
figure(4);
hold on
% Loop through each initial displacement value
for x0 = x0_values
    [t,z] = ode45(@(t,z)[z(2);-(v/m) * z(2) - (a/m) * z(1) + (b/m) * z(1)^3],t,[x0;y0]);
    x = z(:,1);
    y = z(:,2);
    % Plot the displacement and velocity for the current initial displacement value
    plot(t, x, 'DisplayName', ['Displacement, x_0 = ', num2str(x0)]);
    plot(t, y, '--', 'DisplayName', ['Velocity, x_0 = ', num2str(x0)]);
end

% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Amplitude');
title('Damped Nonlinear Oscillator Response for Different Initial Displacements');
legend('show');
grid on;
hold off;

Answer 1

[Sam Chak]

Hi @Athanasios Paraskevopoulos

You could have solved the system's equations to determine the critical points of the Nonlinear Oscillator and then selected appropriate initial values. Previously, the selected initial displacement was outside the stability region, which caused the trajectory to diverge. Here is a demo for reference. I have also corrected the information displayed in the legends.

%% Determine the critical points
syms x y
m       = 1;        % mass
v       = 1;        % damping coefficient
a       = 1;        % linear stiffness
b       = a/20;     % nonlinear stiffness
eq1     = y == 0;
eq2     = (- v*y - a*x + b*x^3)/m == 0;
eqn     = [eq1, eq2];
sol     = solve(eqn);
xcrit   = double(sol.x)

xcrit = 3x1

0 -4.4721 4.4721

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% Damped Nonlinear Oscillator
% Define common parameters
T       = 10;                       % Total simulation time
dt      = 0.1;                      % Time step
t       = 0:dt:T;                   % Time vector
x0      = [1, xcrit(3), 4.4723];    % Different initial displacements
y0      = 0; % Initial velocity

% Initialize a figure
figure(4);
hold on
% Loop through each initial displacement value
for j = 1:length(x0)
    [t, z]  = ode45(@(t, z) [z(2); ...
                           -(v/m)*z(2) - (a/m)*z(1) + (b/m)*z(1)^3], t, [x0(j); y0]);
    x       = z(:,1);
    y       = z(:,2);
    % Plot the displacement and velocity for the current initial displacement value
    plot(t, x,       'DisplayName', ['init disp, x_0 = ', num2str(x0(j))]);
    plot(t, y, '--', 'DisplayName', ['init vel,  y_0 = ', num2str(y0)]);
end
% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Amplitude');
title('Responses of Nonlinear Oscillator for different initial values');
legend('show', 'location', 'best');
grid on;
hold off;