plotting the bifurcation diagram of Duffing oscillations in the (x1, omega) plane?

function dx = duffing(t,x);

global alpha delta beta F omega

dx=[x(2);

-delta*x(2)-alpha*x(1)-beta*x(1)^3+F*cos(x(3));

omega];

% Define parameters

global alpha delta beta F omega

alpha=0.1;

delta=0.2;

beta=0.3;

F=0.5;

omega_range = 0.8:0.1:1; % Range of omega values

x0 = [0.1, 0.1, 0.1]; % Initial conditions

% Set simulation step

dt = 0.001;

% Presave a matrix of NaNs to save points of intersection

M = NaN * zeros(1000, length(omega_range));

pos = 0; % Indexing dummy variable

% Set ode options

options = odeset('RelTol', 1e-5, 'AbsTol', 1e-5);

for omega = omega_range

omega % Print omega, just to track progress

pos = pos + 1; % Increase index

% Simulate the system

[t, x] = ode45(@duffing, 0:dt:1000, x0, options);

% Discard transient

index = t > 400;

X = x(index, :); % Save states without transient in a new vector

l = length(X);

% Save the final x(1) value

M(1:l, pos) = X(:, 1);

end

% Plot the bifurcation diagram

figure;

plot(omega_range, M, '.b', 'MarkerSize', 2);

xlabel('\omega');

ylabel('x_1');

title('Bifurcation Diagram of Duffing System (x_1 vs. \omega)');

grid on;

Can anyone help troubleshoot my code for plotting the bifurcation diagram of Duffing oscillations in the (x1, omega) plane? It's not functioning as expected. Thanks!

Answer 1

[Athanasios Paraskevopoulos]

function dx = duffing(t, x)
    global alpha delta beta F omega
    dx = [x(2);
          -delta * x(2) - alpha * x(1) - beta * x(1)^3 + F * cos(omega * t);
          omega];
end

% Define parameters
global alpha delta beta F omega
alpha = 0.1;
delta = 0.2;
beta = 0.3;
F = 0.5;
omega_range = 0.8:0.01:1.0; % Finer range of omega values
x0 = [0.1, 0.1, 0.1]; % Initial conditions

% Set simulation step and time
dt = 0.01;
simulation_time = 1000; % Total simulation time
transient_time = 400; % Time to discard transient

% Preallocate matrix for results
num_points = floor((simulation_time - transient_time) / dt) + 1;
M = NaN(num_points, length(omega_range));

% Set ODE options
options = odeset('RelTol', 1e-5, 'AbsTol', 1e-5);

% Loop over omega values
for pos = 1:length(omega_range)
    omega = omega_range(pos);
    fprintf('Processing omega = %.2f\n', omega); % Track progress

    % Simulate the system
    [t, x] = ode45(@duffing, 0:dt:simulation_time, x0, options);
    
    % Remove transients
    index = t > transient_time;
    X = x(index, :);
    
    % Save the x(1) values after transient
    M(1:length(X), pos) = X(:, 1);
end

Processing omega = 0.80
Processing omega = 0.81
Processing omega = 0.82
Processing omega = 0.83
Processing omega = 0.84
Processing omega = 0.85
Processing omega = 0.86
Processing omega = 0.87
Processing omega = 0.88
Processing omega = 0.89
Processing omega = 0.90
Processing omega = 0.91
Processing omega = 0.92
Processing omega = 0.93
Processing omega = 0.94
Processing omega = 0.95
Processing omega = 0.96
Processing omega = 0.97
Processing omega = 0.98
Processing omega = 0.99
Processing omega = 1.00

% Plot the bifurcation diagram
figure;
hold on;
for pos = 1:length(omega_range)
    omega_vals = omega_range(pos) * ones(sum(~isnan(M(:, pos))), 1);
    plot(omega_vals, M(~isnan(M(:, pos)), pos), '.b', 'MarkerSize', 2);
end
xlabel('\omega');
ylabel('x_1');
title('Bifurcation Diagram of Duffing System (x_1 vs. \omega)');
grid on;
hold off;