function doall=pendulum_self_sustained
clear all; close all;
% Eli Tziperman: example for course in nonlinear dynamics and
% chaos. 2002.

% A modified pendulum equation that allows for self-sustained
% oscillations even in the absence of forcing. This allows us to see
% the full quasi-periodicity route to chaos, including the stage of
% the torus in phase space, unlike thenormal pendulum which cannot
% display the torus. Change the forcing amplitude f below to see all
% stages in the route to chaos. Solved numerically using the "modified
% Euler scheme" (a second order, constant step scheme), see Strogatz
% page 33.

% The equation:
% d^2\theta/dt^2=(g/L)*\sin(\theta) 
%        + gamma*([d\theta/dt]-[d\theta/dt]^3) + f*cos(omega_f*t)
% The +gamma*d\theta/dt leads to exponential growth rather than decay,
% and the growth is then stopped by the cubic term once the amplitude
% grows to make it non-negligible.


%  First, translate to a dynamical system of three first order ODEs:
%  dx/dt = y
%  dy/dt = -(g/L)*sin(x) + gamma*(y-y^3) + f*cos(omega_f*z)
%  dz/dt = 1

% define parameters:
% ------------------
% pendulum parameters: gravity and length:
g=9.8;
L=9.8;
% Linear frequency:
omega_L=sqrt(g/L);
% Strength of friction (gamma1=1.0; gamma3=0.1; for chaos):
gamma1=1.0;
gamma3=0.1;
% f is forcing amplitude (3.0 gives chaos; 2.0 mode locked; 1.0 torus; 0.0 limit cycle):
f=1.0;
% Forcing frequency (1.5*sqrt(g/L) for chaos):
omega_f=1.5*sqrt(g/L);
% time step: needs to be small enough for stability. Make it also a
% multiple of the forcing period, so that we can then use delay
% coordinates for phase spare reconstruction based on a delay that's
% closely related to teh forcing period:
delta_t=0.01*2*pi/omega_f;
% Print (and save) every nprint steps:
nprint=1;
% Integration time:
T=1500.0;

% Number of steps:
N=floor(T/delta_t);
NP=floor(N/nprint);

% Initial conditions:
x=zeros(3,NP);
n=1;
x(:,n)=[pi/8.0,0,0]';

% Integrate:
for n=1:NP
  xp=x(:,n);
  for j=1:1:nprint
    x_tilde=xp+delta_t*rhs(xp,g,L,gamma1,gamma3,omega_f,f);
    xp=xp+0.5*(rhs(xp,g,L,gamma1,gamma3,omega_f,f)  ...
               + rhs(x_tilde,g,L,gamma1,gamma3,omega_f,f))*delta_t;
  end
  x(:,n+1)=xp;
end

% Plot time series of solution:
figure(1); clf
plot(x(3,:),x(2,:))
xlabel('time')
ylabel('d theta/dt (rad/time)')


% Plot reconstructed phase space using delay coordinates:
figure(2); clf
% make the delay equal to the forcing period
delay=floor(0.2*2*pi/(delta_t*omega_f));
b=floor(NP/2);
e=NP-3*delay;
plot3(x(2,b:e),x(2,b+delay:e+delay),x(2,b+2*delay:e+2*delay))
title('click ''Rotate 3D'' to see the attractor from all angles')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ff=rhs(x,g,L,gamma1,gamma3,omega_f,f)
% return the rhs of the dynamical system:

ff=[x(2),-(g/L)*sin(x(1)) + gamma1*x(2)-gamma3*x(2)^3 + f*cos(omega_f*x(3)),1.0]';
