% Phase Space model describing motion of a simple driven pendulum, driven by 
% a sinusoidal driving force. 
% This code illustrates the phase space of this system and illustrates the 
% concept of 'basin of attraction'
% ========================================================================
% Nonlinear Dynamical Systems
% date : January 23,2006, Monday
% IES: http://www.cavs.msstate.edu/hse/ies/projects/nsf_nonlinear
% Created by: Saurabh Prasad and Sundararajan Srinivasan
% ========================================================================

clear;
clc;
% Setting System Parameters
wD = 2/3;
q = 2;
g = 1.15;
N = 10000;
% Setting Initial Conditions
x1 = 0.5;
y1 = 0.5;
z1 = 0.5;
dt = 0.01;
x1_cat = [];
y1_cat = [];
z1_cat = [];

% The following FOR loop uses Runge-Kutta Numerical Integration to generate the observables of a simple pendulum
  for i = 1 : 1 : N  

    dx1 = ( -x1/q - sin(y1) + g*cos(z1) ) * dt;
    dy1 = ( x1 ) * dt;
    dz1 = ( wD ) * dt;

    dx2 = ( -(x1+dx1/2)/q - sin(y1+dy1/2) + g*cos(z1+dz1/2) ) * dt;
    dy2 = ( x1+dx1/2 ) * dt;
    dz2 = ( wD ) * dt;

    dx3 = ( -(x1+dx2/3)/q - sin(y1+dy2/3) + g*cos(z1+dz2/3) ) * dt;
    dy3 = ( x1+dx2/3 ) * dt;
    dz3 = ( wD ) * dt;

    dx4 = ( -(x1+dx3)/q - sin(y1+dy3) + g*cos(z1+dz3) ) * dt;
    dy4 = ( x1+dx3 ) * dt;
    dz4 = ( wD ) * dt;

    x1 = x1 + dx1/ 6 + dx2/ 3 +dx3/ 3 +dx4/ 6;    
    y1 = y1 + dy1/ 6 + dy2/ 3 +dy3/ 3 +dy4/ 6;    
    z1 = z1 + dz1/ 6 + dz2/ 3 +dz3/ 3 +dz4/ 6;
    
    y1 = mod(y1, 2*pi);
    if y1 > pi
        y1 = - (2*pi - y1);
    end;
    
    x1_cat = [x1_cat x1];
    y1_cat = [y1_cat y1];
    z1_cat = [z1_cat z1];
  
  end;
  
  
plot(y1_cat(5000:end),x1_cat(5000:end), '.');  
  
x1 = .80;
y1 = .80;
z1 = .80;
dt = 0.01;
x1_cat = [];
y1_cat = [];
z1_cat = [];

% The following FOR loop uses Runge-Kutta Numerical Integration to generate the observables of a simple pendulum 
  for i = 1 : 1 : N  

    dx1 = ( -x1/q - sin(y1) + g*cos(z1) ) * dt;
    dy1 = ( x1 ) * dt;
    dz1 = ( wD ) * dt;

    dx2 = ( -(x1+dx1/2)/q - sin(y1+dy1/2) + g*cos(z1+dz1/2) ) * dt;
    dy2 = ( x1+dx1/2 ) * dt;
    dz2 = ( wD ) * dt;

    dx3 = ( -(x1+dx2/3)/q - sin(y1+dy2/3) + g*cos(z1+dz2/3) ) * dt;
    dy3 = ( x1+dx2/3 ) * dt;
    dz3 = ( wD ) * dt;

    dx4 = ( -(x1+dx3)/q - sin(y1+dy3) + g*cos(z1+dz3) ) * dt;
    dy4 = ( x1+dx3 ) * dt;
    dz4 = ( wD ) * dt;

    x1 = x1 + dx1/ 6 + dx2/ 3 +dx3/ 3 +dx4/ 6;    
    y1 = y1 + dy1/ 6 + dy2/ 3 +dy3/ 3 +dy4/ 6;    
    z1 = z1 + dz1/ 6 + dz2/ 3 +dz3/ 3 +dz4/ 6;
    
    y1 = mod(y1, 2*pi);
    if y1 > pi
        y1 = - (2*pi - y1);
    end;
    
    x1_cat = [x1_cat x1];
    y1_cat = [y1_cat y1];
    z1_cat = [z1_cat z1];
  
  end;
  
  hold on;
  plot(y1_cat(5000:end),x1_cat(5000:end),'r.');
  title('The dynamical system started from two different initial conditions converges to the basin of attraction in the phase space.')
  xlabel('theta')
  ylabel('omega')

