%%
%% 7.3 Duffing oscillator (Computational example)
% The Duffing oscillator is a one-dimensional driven oscillator with friction and a nonlinear potential. It is described by the
% Equation of motion: \ddot x+r\dot x=K+f \sin(\Omega t) where K=-{\partial \over \partial x}V and V(r)=-{1\over 2}\omega^2x^2+{1\over 4}\beta x^4
% The equation of motion is solved numerically and the solutions are examined for attractive or chaotic behavior.
% Considering:
% \circ Phase space trajectories (x(t),v(t))
% \circ Divergence and convergence of adjacent trajectories
% \circ Poincare cuts
% Depending on the drive amplitude f, the oscillator has different attractors (1-cycle, 2-cycle, 4-cycle, ...) and chaotic behavior.
%

clear all
syms omega Omega beta r x f x_0 v_0
Par=[r==0.25 omega==1 beta==1 Omega==1.4 x_0==0.001,v_0==0]
F=[0.34 0.38 0.388 0.4]; %Drive amplitude
T=600;
N=4000;
Tn=linspace(0,T,N);
opt=odeset('RelTol',3e-14);

%
%% 1 Potential of the Duffing oscillator
% V=-{1\over 2}\omega^2x^2+{1\over 4}\beta x^4


V=-omega^2/2*x^2+beta*x^4/4
sube(V,Par)
fplot(ans,[-1.5 1.5])
xlabel('Position x')
title('Potential {\sl V}')

%
%% 2 Equations of motion
% \ddot x+r\dot x=K+f \sin(\Omega t) mit $K =  - {\partial  \over {\partial x}}V \to \left\{ {\matrix{
%    {\dot x = v} \quad  \qquad  \qquad  \qquad  \qquad  \qquad  \cr
%    {\dot v =  - rv + {\omega ^2}x - \beta {x^3} + f\sin \left( {\Omega t} \right)} \cr
%
%  } } \right.$

syms x y(t) v(t)
K=-sube(diff(V,x),x==y(t));
syms x(t)
K=sube(K,y(t)==x(t))
diff(x(t),t,2)+r*diff(x(t),t)==K+f*sin(Omega*t)
isolate(ans,diff(x(t),t,2))
DGL=[v(t),sube(rhs(ans),diff(x(t),t)==v(t))]

%
%% 3 Numerical solution and phase space representation

sube(DGL,Par)
Fun=odeFunction(ans,[x(t) v(t)],f);
xv1=double(sube([x_0 v_0],Par));
xv2=double(sube([x_0+1e-7 v_0],Par));
l=0.8;
figure
for n=1:4
    [t_n,x1]=ode45(@(t,x)Fun(t,x,F(n)),Tn,xv1,opt);
    %subplot(2,2,n)
    subplot(4,4,n)
    plot([0 2*l],[0 0],':k')
    hold on
    plot([0 0],[-0.7 0.7],':k')
    plot(x1(1000:N,1),x1(1000:N,2),'b')
    text(0.8,l+.1,sprintf('Drive f=%3.3f',F(n)),...
        'HorizontalAlignment','center','Color','k')
    axis([0 1.5 -.7 .7])
    axis equal
    xlabel('Displacement x(t)')
    ylabel('Velocity v(t)')
    hold off
end

% For the four drive amplitudes considered, the results are 1-cycle, 2-cycle, 4-cycle, and chaotic behavior with a strange attractor.
%
%% 4 Convergence and divergence of the phase space trajectories
% Two trajectories of the Duffing oscillator are considered, which initially have a small spatial distance of x_2(0)-x_1(0)=10^{-7} and the velocities v_1(0)=v_2(0)=0.
% The distance in phase space is $dx\left( t \right) = \left| {\left( {\matrix{
%    {{x_2}\left( r \right) - {x_1}\left( r \right)}  \cr
%    {{v_2}\left( r \right) - {v_1}\left( r \right)}  \cr
%
%  } } \right)} \right|$
% The Lyapunov exponent $\lambda $ is defined by dx(t)=dx(0)e^{\lambda t}

n1=[1 1 1 1];
n2=[0.5*N 0.3*N 0.6*N N/2];
for n=1:4
    [t_n,x1]=ode45(@(t,x)Fun(t,x,F(n)),Tn,xv1,opt);
    [t_n,x2]=ode45(@(t,x)Fun(t,x,F(n)),Tn,xv2,opt);
    subplot(2,2,n)
    dx=log10((x1(:,1)-x2(:,1)).^2+(x1(:,2)-x2(:,2)).^2)/2;
    plot(t_n,dx)
    k=n1(n):n2(n);
    p=polyfit(t_n(k),dx(k),1);
    hold on
    plot([0 T],polyval(p,[0 T]),'r')
    text(300,-6,sprintf('Drive f=%3.3f',F(n)),...
        'HorizontalAlignment','center','Color','b')
    text(300,-8,['\lambda = ',sprintf('%3.3f',p(1)*log(10))],...
        'HorizontalAlignment','center','Color','r')
    axis([0 T -15 0])
    xlabel('Time t')
    ylabel('log(dx)')
    hold off
end

%
%% 5 Poincare cuts
% The displacement x(t_n) of the oscillator at times t_n={2\pi n \over \Omega} (n: cycles) is depicted.
% This leads to a stroboscopic effect

N4=100;
dt=2*pi/Omega
dt=double(sube(dt,Par));
Tn4=0:dt:T;
for n=1:4
    [t_n,x1]=ode45(@(t,x)Fun(t,x,F(n)),Tn4,xv1,opt);
    subplot(2,2,n)
    plot(t_n,x1(:,1),'o','MarkerSize',2)
    xlabel('Time t')
    ylabel('x(t_n)')
    title(sprintf('Drive f=%3.3f',F(n)))
end