%% The Duffing Equation
% The Duffing equation is an externally forced and damped oscillator equation. It is given by\ddot{x} +\delta \;\dot{x} +\alpha \;x+\beta \;x^3 =\gamma \;\cos \;\omega \;t.
%
% We will investigate the equation by turning on a few constants at a time. We begin with simple harmonic motion: \ddot{x} +\alpha \;x=0.

clear
alpha = 1; delta = 0; beta = 0; % SHM
duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3];
tspan = [0 30];
y0 = [0 0.01]';
[t,y] = ode45(duffing,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Simple Harmonic Motion')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('SHM Trajectory')
axis equal

%
% Now add damping, \delta >0\ldotp

clear
alpha = 1; delta = 0.1; beta = 0;  % damped SHM
duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3];
tspan = [0 30];  % tspan = 0:.1:30;  % Makes solution smoother
y0 = [0 0.1]';
[t,y] = ode45(duffing,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Damped Simple Harmonic Motion')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('Damped SHM Trajectory')
axis equal

%
% Now, consider the undamped Duffing equation \ddot{x} +\alpha \;x+\beta \;x^3 =0\ldotp As a first order system, we have \dot{x}=y\\
% \dot{y}=-\alpha x-\beta x^3. Note the equilibrium points are \left(0,0\right),\left(\pm \sqrt{\frac{-\alpha }{\beta }},0\right).

clear
alpha = -1; delta = 0; beta = 1; % Undamped Duffing Equation
duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3];
tspan = [0 10];  
y0 = [1 1]';
[t,y] = ode45(duffing,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Undamped, Unforced Duffing Equation')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('Duffing Equation Trajectory')
axis equal

%
% Its useful to look at the phase plot using the quiver function. N is used to normalize the arrows in the direction field. The first plot shows one trajectory.

[X, Y] = meshgrid(-1.5:0.25:1.5, -1.5:0.25:1.5);
U = Y; V=-delta*Y-alpha*X-beta*X.^3;
N=sqrt(U.^2+V.^2);
quiver(X, Y, U./N, V./N)

hold on
plot(y(:,1),y(:,2),'r')
hold off
xlabel('x')
ylabel('y')
title('Duffing Equation Quiver Plot')

%
% The following shows how to add several trajectories.

hold on
for j=1:5
    y0 = [1 j/5]';
    [t,y] = ode45(duffing,tspan,y0);
    plot(y(:,1),y(:,2),'r')
end
hold off

%
% Now, we add back some damping. Changing the initial condition, y0 = [1 1]'; will reveal different types of behavior.

clear
alpha = -1; delta = 0.1; beta = 1;  % Damped Duffing Equation
duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3];
tspan = [0 30];  
y0 = [1 1]';
[t,y] = ode45(duffing,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Damped Duffing Equation')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('Duffing Equation Trajectory')
axis equal

%
% Here we add the forcing.

clear
alpha = -1; delta = 0.1; beta = 1; gamma = .1; omega=1.4; % Forced DE
duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3+gamma*cos(omega*t)];
tspan = 0:0.1:200;  
y0 = [0,0]';
[t,y] = ode45(duffing,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Damped, forced Duffing Equation')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('Duffing Equation Trajectory')
axis equal

%
% We can examine the tail of the trajectory and notice a periodic behavior with period T=\frac{2\pi }{\omega }=4\ldotp 488\ldotp

% Plot tail of trajectory
figure
L=length(y(:,1));
Npts=45;         % Npts+1 points of the tail
plot(y(L-Npts:L,1),y(L-Npts:L,2))
%plot(y(L-44:L,1),y(L-44:L,2))  % roughly 45*.1=period
xlabel('x')
ylabel('y')
title(['Tail of Trajectory \gamma = ',num2str(gamma)])
axis equal
T=2*pi/omega

%
% If we change the amplitude to \gamma =0\ldotp 318, then we notice period doubling.

clear   % Change gamma to 0.318 - Period Doubling to 2 cycles
alpha = -1; delta = 0.1; beta = 1; gamma = .318; omega=1.4; % Forced DE
duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3+gamma*cos(omega*t)];
NT = 800;  
tspan = 0:0.1:NT;  
y0 = [0,0]';
[t,y] = ode45(duffing,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Undamped, Unforced Duffing Equation')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('Duffing Equation Trajectory')
axis equal

% Plot tail of trajectory
figure
L=length(y(:,1));
Npts=180;         % Npts+1 points of the tail
%plot(y(L-Npts:L,1),y(L-Npts:L,2))
plot(y(L-90:L,1),y(L-90:L,2))  % 90*.1=period - double before
xlabel('x')
ylabel('y')
title(['Tail of Trajectory \gamma = ',num2str(gamma)])
axis equal

% Try gamma = 0.338 but 0.35 is chaotic region

%
% For \gamma =0\ldotp 338, we get four cycles, or period 4.

clear   % Change gamma to 0.338 - Period Doubling to 4 cycles
alpha = -1; delta = 0.1; beta = 1; gamma = .338; omega=1.4; % Forced DE
duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3+gamma*cos(omega*t)];
NT = 2000; 
tspan = 0:0.1:NT;  
y0 = [0,0]';
[t,y] = ode45(duffing,tspan,y0);

figure
L=length(y(:,1));
Npts=180;         % Npts+1 points of the tail
plot(y(L-Npts:L,1),y(L-Npts:L,2))
xlabel('x')
ylabel('y')
title(['Tail of Trajectory \gamma = ',num2str(gamma)])
axis equal

%
% Can we find period 8? Maybe near \gamma =0\ldotp 3431\ldotp

clear   % - Period Doubling to more cycles?
alpha = -1; delta = 0.1; beta = 1; gamma = .3431; omega=1.4; % Forced DE
duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3+gamma*cos(omega*t)];
NT = 6400; 
tspan = 0:0.1:NT;  
y0 = [0,0]';
[t,y] = ode45(duffing,tspan,y0);

figure
L=length(y(:,1));
Npts=360;         % Npts+1 points of the tail
plot(y(L-Npts:L,1),y(L-Npts:L,2))
xlabel('x')
ylabel('y')
title(['Tail of Trajectory \gamma = ',num2str(gamma)])
axis equal

%
% Possibly for \gamma =0\ldotp 35 there is only chaos.

clear   % - Chaos?
alpha = -1; delta = 0.1; beta = 1; gamma = .35; omega=1.4; % Forced DE
duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3+gamma*cos(omega*t)];
NT = 10000; 
tspan = 0:0.1:NT;  
y0 = [0,0]';
[t,y] = ode45(duffing,tspan,y0);

figure
L=length(y(:,1));
Npts=NT;         % Npts+1 points of the tail
plot(y(L-Npts:L,1),y(L-Npts:L,2))
xlabel('x')
ylabel('y')
title(['Tail of Trajectory \gamma = ',num2str(gamma)])
axis equal

%
% Let's look at a Poincaré Surface of Section for \gamma =0\ldotp 35.

% Strange Attractor
clear   
alpha = -1; delta = 0.1; beta = 1; gamma = .35; omega=1.4; 
%alpha = -1; delta = 0.25; beta = 1; gamma = .40; omega=1; % Standard Shape
%alpha = -1; delta = 0.2; beta = 1; gamma = .30; omega=1; % Standard Shape

duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3+gamma*cos(omega*t)];
NT = 6000;  % or 2000 for gamma = 0.338
period=2*pi/omega;
dt=period/10;
tspan = 0:dt:20000;  
y0 = [0,0]';
[t,y] = ode45(duffing,tspan,y0);


P_x = y(mod(t, period) < 1e-6, 1);
P_y = y(mod(t, period) < 1e-6, 2);

figure
scatter(P_x,P_y,'.')
xlabel('Position');
ylabel('Velocity');
title('Poincaré Surface of Section');

%
% Here we add more initial conditions and find the same attractor. In fact, it is a strange attractor.

% Add more ICs
figure
scatter(P_x,P_y,'.')
xlabel('Position');
ylabel('Velocity');
title('Poincaré Surface of Section');
hold on
for j=1:50
    y0 = [j/50,0]';
    [t,y] = ode45(duffing,tspan,y0);
    P_x = y(mod(t, period) < 1e-6, 1);
    P_y = y(mod(t, period) < 1e-6, 2);
    L=length(P_y);
    Npts=180;
   % scatter(P_x,P_y,'.r')
    scatter(P_x(L-Npts:L),P_y(L-Npts:L),'.r')    
end
hold off

%
% A more standard strange attractor for the Duffing equation is found using a unit amplidute. Here are two possibles cases with varying parameters.

% Strange Attractor
clear   
alpha = -1; delta = 0.25; beta = 1; gamma = .40; omega=1; % Standard Shape
%alpha = -1; delta = 0.2; beta = 1; gamma = .30; omega=1; % Standard Shape

duffing = @(t,y) [y(2); -delta*y(2)-alpha*y(1)-beta*y(1).^3+gamma*cos(omega*t)];
NT = 6000;  % or 2000 for gamma = 0.338
period=2*pi/omega;
dt=period/10;
tspan = 0:dt:20000;  
y0 = [0,0]';
[t,y] = ode45(duffing,tspan,y0);


P_x = y(mod(t, period) < 1e-6, 1);
P_y = y(mod(t, period) < 1e-6, 2);

figure
scatter(P_x,P_y,'.')
xlabel('Position');
ylabel('Velocity');
title('Poincaré Surface of Section');