Duffing equation:Transition to Chaos

The Original Equation is the following:

Let . This implies that

Then we rewrite it as a System of First-Order Equations

Using the substitution for , the second-order equation can be transformed into the following system of first-order equations:

Exploring the Effect of γ.

• 

% Define parameters
gamma = 0.338;
alpha = -1;   
beta = 1; 
delta = 0.1;
omega = 1.4;   

% Define the system of equations
odeSystem = @(t, y) [y(2); 
                     -delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];

% Initial conditions
y0 = [0; 0];  % x(0) = 0, v(0) = 0

% Time span
tspan = [0 2000];

% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);

% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;

% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase-Plane  $$\ddot{x}+\delta \dot{x}+\alpha x+\beta x^3=0$$ for $$\gamma=0.318$$, $$\alpha=-1$$,$$\beta=1$$,$$\delta=0.1$$,$$\omega=1.4$$','Interpreter', 'latex');
grid on;

% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t));  % Starting index for the tail
tail_end = length(t);  % Ending index for the tail

% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase-Plane  $$\ddot{x}+\delta \dot{x}+\alpha x+\beta x^3=0$$ for $$\gamma=0.318$$, $$\alpha=-1$$,$$\beta=1$$,$$\delta=0.1$$,$$\omega=1.4$$','Interpreter', 'latex');
grid on;

ax = gca;
chart = ax.Children(1);
datatip(chart,0.5581,-0.1126);

Then I used but the results were not that I expected for

. The paper I study is Duffing Equation

My code gives me the following. Any suggestion?

% Define parameters
gamma = 0.35;
alpha = -1;   
beta = 1; 
delta = 0.1;
omega = 1.4;   

% Define the system of equations
odeSystem = @(t, y) [y(2); 
                     -delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];

% Initial conditions
y0 = [0; 0];  % x(0) = 0, v(0) = 0

% Time span
tspan = [0 3000];

% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);




% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase-Plane  $$\ddot{x}+\delta \dot{x}+\alpha x+\beta x^3=0$$ for $$\gamma=0.318$$, $$\alpha=-1$$,$$\beta=1$$,$$\delta=0.1$$,$$\omega=1.4$$','Interpreter', 'latex');
grid on;

% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t));  % Starting index for the tail
tail_end = length(t);  % Ending index for the tail

% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase-Plane  $$\ddot{x}+\delta \dot{x}+\alpha x+\beta x^3=0$$ for $$\gamma=0.35$$, $$\alpha=-1$$,$$\beta=1$$,$$\delta=0.1$$,$$\omega=1.4$$','Interpreter', 'latex');
grid on;

ax = gca;
chart = ax.Children(1);
datatip(chart,0.5581,-0.1126);

[Neelanshu]

Hello Athanasios,

Could you please share the expected output for the plot with changes in gamma? This would help in investigating the reason for the discrepancy.

[Athanasios Paraskevopoulos]

Hello @Neelanshu. The plots that I need are the following

[Athanasios Paraskevopoulos]

Hello @Neelanshu

I tried

% Time span with more points for better resolution
tspan = linspace(0, 200,3000);  % Increase the number of points

I feel that I get closer but still I can not create the other two plots yet

% Define parameters
gamma = 0.35;
alpha = -1;   
beta = 1; 
delta = 0.1;
omega = 1.4;   

% Define the system of equations
odeSystem = @(t, y) [y(2); 
                     -delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];

% Initial conditions
y0 = [0; 0];  % x(0) = 0, v(0) = 0

% Time span with more points for better resolution
tspan = linspace(0, 200,3000);  % Increase the number of points

% Solve the system with increased output refinement
options = odeset('Refine', 4);  % Additional refinement of the output
[t, y] = ode45(odeSystem, tspan, y0, options);

% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;

% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase-Plane  $$\ddot{x}+\delta \dot{x}+\alpha x+\beta x^3=0$$ for $$\gamma=0.318$$, $$\alpha=-1$$,$$\beta=1$$,$$\delta=0.1$$,$$\omega=1.4$$','Interpreter', 'latex');
grid on;

% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t));  % Starting index for the tail
tail_end = length(t);  % Ending index for the tail

% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase-Plane  $$\ddot{x}+\delta \dot{x}+\alpha x+\beta x^3=0$$ for $$\gamma=0.318$$, $$\alpha=-1$$,$$\beta=1$$,$$\delta=0.1$$,$$\omega=1.4$$','Interpreter', 'latex');
grid on;

% Adding a datatip for visualization (optional)
ax = gca;
chart = ax.Children(1);
datatip(chart,0.5581,-0.1126);