Need Bifurcation Diagram from code please.

% Define the system of ODEs

function f = harmonic_oscillator(t, z, omega, gamma, gamma_drive)

x = z(1);

x_dot = z(2);

f = [x_dot; -2 * gamma * x_dot - omega^2 * x + gamma_drive * cos(1.25 * t)];

end

% Set the parameters

omega = 1.25;

gamma_values = [0.2, 0.3, 0.315, 0.37, 0.5, 0.8];

t_span = [0, 400];

% Loop through different gamma values

for gamma = gamma_values

% Set the initial conditions

initial_conditions = [1; 0];

% Define the function for the ODE solver

ode_func = @(t, z) harmonic_oscillator(t, z, omega, gamma, gamma);

% Use the ODE solver to obtain the solution

[t, sol] = ode45(ode_func, t_span, initial_conditions);

% Extract the solution

x = sol(:, 1);

x_dot = sol(:, 2);

% Plot the numerical solution

figure('Position', [100, 100, 1200, 400]);

subplot(1, 2, 1);

plot(t, x, 'DisplayName', 'x(t)');

xlabel('Time');

ylabel('Displacement');

title(['Numerical Solution for Gamma = ', num2str(gamma)]);

legend;

% Plot the phase portrait

subplot(1, 2, 2);

plot(x, x_dot, 'DisplayName', 'Phase Portrait');

xlabel('x');

ylabel('x_dot');

title(['Phase Portrait for Gamma = ', num2str(gamma)]);

legend;

drawnow;

end

This is my code and I am not sure on how to do the bifurcation plots now.

Answer 1

[Shubham]

Hi Rony,

To create a bifurcation diagram showing how the steady-state amplitude changes with the damping coefficient gamma, follow these steps:

• Simulate the harmonic oscillator for a range of gamma values. This will allow you to observe how the system's behavior changes with different damping.
• Remove transient oscillations by focusing on the solution after a certain time (e.g., t > 300). This ensures you're analyzing the steady-state behavior.
• Identify the peaks of x(t) in the steady-state. These peaks represent the maximum displacement values.
• Compute the average of these peaks for each gamma value and plot them. This will form your bifurcation diagram.

Here's the modified code for achieve this:

% Define the system of ODEs
function f = harmonic_oscillator(t, z, omega, gamma, gamma_drive)
    x = z(1);
    x_dot = z(2);
    f = [x_dot; -2 * gamma * x_dot - omega^2 * x + gamma_drive * cos(1.25 * t)];
end

% Parameters
omega = 1.25;                          % Natural frequency
gamma_values = linspace(0.2, 0.8, 50); % Range of damping coefficients
t_span = [0, 400];                     % Time range for simulation
num_skip = 300;                        % Ignore first 300 seconds (transients)
amplitudes = [];                       % Array to store steady-state amplitudes

% Loop through each gamma
for gamma = gamma_values
    % Initial conditions
    initial_conditions = [1; 0];
    % Define ODE function
    ode_func = @(t, z) harmonic_oscillator(t, z, omega, gamma, gamma);
    % Solve the ODE
    [t, sol] = ode45(ode_func, t_span, initial_conditions);
    
    % Extract steady-state portion
    x = sol(t > num_skip, 1); % Displacement after removing transients
    
    % Find peaks in the steady-state solution
    [peaks, ~] = findpeaks(x);
    amplitudes = [amplitudes; mean(peaks)]; % Average of steady-state peaks
end

% Plot bifurcation diagram
figure;
plot(gamma_values, amplitudes, 'o-', 'LineWidth', 1.5, 'MarkerSize', 6);
xlabel('Gamma (Damping Coefficient)');
ylabel('Steady-State Amplitude');
title('Bifurcation Diagram');
grid on;

Hope this helps.