Plotting Damped Nonlinear Oscillator
Show older comments
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?
, 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;
4 Comments
Athanasios Paraskevopoulos
on 7 May 2024
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
% 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;
You are right for 
the code works. Thank you!
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;
Accepted Answer
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.
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)
%% 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;
More Answers (0)
Categories
Find more on Time and Frequency Domain Analysis 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!
