% This file allows to plot some nonlinear systems: 
% Lorenz, Duffing, Chua circuit, Colpitts and Rayleigh
% oscillators, Van der Pol oscillator (for µ=0.2, µ=1, µ=10), 
% and Henon map.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [1] Plot Lorenz dynamic by using ode45 command:
%
%    Lorenz's equations:
%           
%       dx/dt = sigma*(y-x)
%       dy/dt = x*(rho-z)-y
%       dz/dt = x*y-beta*z
%
%    with sigma=10, rho=28., beta=8./3,
%    Initial conditions x0=[10 10 10].

disp('[1] Lorenz ')
disp(' ')
disp(' ... Please Wait ...')
disp(' ')
disp('Then Press any Key to Continue ...')
disp(' ')
disp(' ')

global sigma r b
sigma=10.;  
r=28.;    
b=8./3.;
x0=[10 10 10];
t0=0; 
tf=100; 
[t, x] = ode45(@lorenzeqs, [t0, tf], x0);     
figure;
plot3(x(:,1),x(:,2),x(:,3)); 
title('Lorenz');
xlabel('x');
ylabel('y');   
zlabel('z');
axis([-10 10 -30 30 0 70]);
box;

pause
clear all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [2] Plot Duffing dynamic in the Phase space by using ode45 command
%
%    Duffing's equations (general form):
%           
%       du/dt = v
%       dv/dt = -(omega_0)^2*u-beta*u^3-delta*v+gamma*cos(omega*t+phi)
%
%    where u is the displacement, v is the velocity, and the 
%    parameters are constants.


disp('[2] Duffing')
disp(' ')
disp(' ... Please Wait ...')
disp(' ')
disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)')
disp(' ')
disp('Press any Key to Continue, ....')
disp(' ')


[t,z]=ode45(@duffing,[0,300],[3.01,4.01]); 
plot(z(:,1),z(:,2));
title('Duffing, phase space');

pause
clear all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [3] Plot Chua circuit dynamic by using ode45 command
%
%    Chua's autonomous circuit
%           
%       dx/dt = alpha*(-x+y-H(x))
%       dy/dt = x-y+z
%       dz/dt = -beta*y
%           
%    where H(x) = bx+0.5*(a-b)*(|x+1|-|x-1|),
%    with alpha arbitrary, alpha is selected here equal to 10, 
%    beta = 15, a = -1.3, b = -0.7


disp('[3] Chua circuit')
disp(' ')
disp(' ... Please Wait ...')
disp(' ')
disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)')
disp(' ')
disp('Press any Key to Continue, ....')
disp(' ')


y0=[1 0.2 0];
[t,z]=ode45(@Chua,[0 300],y0); 
plot(z(:,1),z(:,2));
title('Chua Circuit');


pause
clear all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [4] Plot Colpitts Oscillator dynamic by using ode45 command
%
%    Colpitts oscillator system:
%           
%       dx/dt = z-beta*f(y)
%       dy/dt = alpha*(-r*(y+gamma)-z-f(y))
%       dz/dt = delta*(-x+y+epsilon)-rho*z
%                                  
%    where f(y)=0 if y<=1, and y-1 otherwise.                   
%    with alpha=1, beta = 200, gamma = -20/3, delta = 5.5,
%    with r arbitrary, rho = 2,epsilon = 20/3.  


disp('[4] Colpitts')
disp(' ')
disp(' ... Please Wait ...')
disp(' ')
disp('Then Press any Key to Continue, ....')
disp(' ')


y0=[0.1 0 0];
s=30;Tm=[0:1/s:200];
[t,y]=ode45(@Colpitts,Tm,y0); 
plot3(y(:,1),y(:,2),y(:,3));
xlabel('Y_1');
ylabel('Y_2');
zlabel('Y_3');
title('Colpitts Oscillator');
axis([0 8 -2.5 2 -2 8])
view(70,25);
box;


pause
clear all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [5] Rayleigh oscillator with physical parameters
%          
%    c*dv/dt = i+id(-v)
%    b*di/dt = -v-u1-u2*sin(omega*t)
%
% where u1=1, u2=878, omega=60.3e3, a1=1.542, a2=0.5,
% b=10, c=2.377 and id(v)=-a1*v+a2*(-v)^3.

disp('[5] Rayleigh oscillator')
disp(' ')
disp(' ... Please Wait ...')
disp(' ')
disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)')
disp(' ')
disp('Press any Key to Continue, ....')
disp(' ')

y0=[0.1 0];
q=10e5;
vtm=[0:1/q:20e-3];
[t,y]=ode45(@Rayleigh,vtm,y0); 
plot(y(:,1),y(:,2));
xlabel('y_1');
ylabel('y_2');
title('Rayleigh Oscillator');


pause
clear all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [6] The van der Pol equation (y_1)"-µ*(1-(x_1)^2)*y'+(y_1)=0 is equivalent 
%    to a system of coupled first-order differential equations
%       (y_1)'=(y_2)
%       (y_2)'=µ*[1-(y_1)^2]*(y_2)-(y_1).
%   Evaluate the van der Pol Ode for mu = 0.2. 
%   The ode command is as follows:


disp('[6] van der Pol oscillator for µ=0.2')
disp(' ')
disp(' ... Please Wait ...')
disp(' ')
disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)')
disp(' ')
disp('Press any Key to Continue ...')
disp(' ')


x0=[0.1 0.1];
[t,y] = ode45(@vdpeq1,[0 100],x0); 
plot(t,y(:,1),'-',t,y(:,2),'-.');
xlabel('time t');
ylabel('solution y_1 and y_2');
title('van der Pol with \mu=0.2');

pause

disp(' ')
disp('[ . 6 bis] ')
disp(' ')
disp('Press any Key to Continue ...')
disp(' ')


plot(y(:,1),y(:,2));
xlabel('y_1');
ylabel('y_2');
title('van der Pol with \mu=0.2');


pause
clear all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [7] The van der Pol equation (y_1)"-µ*(1-(x_1)^2)*y'+(y_1)=0 is equivalent 
%    to a system of coupled first-order differential equations
%       (y_1)'=(y_2)
%       (y_2)'=µ*[1-(y_1)^2]*(y_2)-(y_1).
%   Evaluate the van der Pol Ode for mu = 1. 
%   The ode command is as follows:


disp('[7] van der Pol oscillator for µ=1')
disp(' ')
disp('Press any Key to Continue ....')
disp(' ')

x0=[2 0];
[t,y] = ode45(@vdpeq2,[0 20],x0); 
plot(t,y(:,1),'-',t,y(:,2),'-.');
xlabel('time t');
ylabel('solution y_1 and y_2');
title('van der Pol with \mu=1');

pause

disp(' ')
disp('[ . 7 bis] ')
disp(' ')
disp('Press any Key to Continue ....')
disp(' ')


plot(y(:,1),y(:,2));
xlabel('y_1');
ylabel('y_2');
title('van der Pol with \mu=1')


pause
clear all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [8] The van der Pol equation (y_1)"-µ*(1-(x_1)^2)*y'+(y_1)=0 is equivalent 
%    to a system of coupled first-order differential equations
%       (y_1)'=(y_2)
%       (y_2)'=µ*[1-(y_1)^2]*(y_2)-(y_1).
%   Evaluate the van der Pol Ode for mu = 10. 
%   The ode command is as follows:

disp(' ')
disp('[8] van der Pol oscillator for µ=10')
disp(' ')
disp('Press any Key to Continue, ....')
disp(' ')

x0=[2 0];
[t,y] = ode15s(@vdpeq3,[0 100],x0); 
plot(t,y(:,1),'-');
xlabel('time t');
ylabel('solution y_1');
title('van der Pol with \mu=10')

pause
disp(' ')
disp('[ . 8 bis] ')
disp(' ')
disp('Press Any Key to Continue ...')
disp(' ')


plot(y(:,1),y(:,2));
xlabel('y_1');
ylabel('y_2');
title('van der Pol with \mu=10')


pause
clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [9] Plot the Henon map by using the function henonmap(iterations);
%   (The number of iterations has to be selected between 100 and 2000).

disp('[9] Henon map')
disp(' ')
disp(' ... Please Wait Few Seconds ...')
disp(' ')
disp('(To Stop the Simulation, Press "Ctrl+c" in the Command Window.)')
disp(' ')


henonmap(1500);   
   
 
% "Complex and Chaotic Nonlinear Dynamics. 
% Advances in Economics and Finance, 
% Mathematics and Statistics" 
% T.Vialar, Springer 2009
% Copyright(c).