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);

[nick]

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);

Answer 1

[Swastik Sarkar]

Hi @Athanasios Paraskevopoulos,

Upon reviewing the paper linked in the question (https://www.colorado.edu/amath/sites/default/files/attached-files/good_sample_project_0.pdf), it specifies a time range of [0, 3000] and focuses on the initial and final few points for its plots, whereas the current solution plots the range [0, 200] with 3000 points, as indicated usinglinspace(0, 200, 3000) upon the entire duration. So, changing the code to the below one will give you desired plots.

% 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 

tspan = [0 3000];  % Increase the number of points 

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

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

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

% Define the head (e.g., first 10% of the time interval) 
head_start = 1;  % Starting index for the tail 
head_end = floor(0.1 * length(t));  % Ending index for the tail 
% Plot the phase portrait 
figure 
plot(y(head_start:head_end, 1), y(head_start:head_end, 2), '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; 

% Plot the tail of the solution 
figure 
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), '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; 

% Plot the extreme tail of the solution 
figure 
plot(y(extreme_tail_start:tail_end, 1), y(extreme_tail_start:tail_end, 2), '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);