# sandbox/easystab/david/PhasePortrait_NonLinear.m

function [] = PhasePortrait_NonLinear()	

# Phase portrait of a NON-LINEAR dynamical system of dimension 2.

This program draws a phase portrait for a dynamical system written under the form

\displaystyle \frac{d}{dt} X(t) = F(X(t)).

With \displaystyle X = \left[ \begin{array}{c} x_1 \\ x_2 \end{array} \right]

and F is a function defined at the end of the program.

The drawing of the phase portrait is done using the same representation as for the linear case treated by the program PhasePortrait_Linear.m.

close all;clear all;
% parametres :
global typeproblem xmin xmax ymin ymax;

typeproblem = 'Pendulum';
% Select the problem to consider among the following cases :
% 'Pitchfork', 'Brusselator', 'VanDerPol', 'LotkaVolterra', 'BuffaloWolf', 'Trefethen',
% 'Custom' [ , ... ]

%plotEnergy=0; % use 1 to plot the energy in figure 2

% Dimensions of the figure. We give default values here. These can also be adjusted in the definition of the function.
xmin = -2; xmax = 2;
ymin = -2; ymax = 2;
Z = F_PortraitdePhase([0.;0.]); % (trick) call to the function to get the values defined there if not using default..
eps = 1e-6;
% Pour le trace de trajectoires en cliquant sur la figure : nombre de trajectoires, Tmax
Ntraj = 10;
Tmax = 50; % max time for integratation of trajectory
Tmaxbak = 5; % max time for integratation of trajectory (backwrards)
Tmaxplot = 5; % max time for plot in figure 2.

% Fin des parametres


## Illustration of the flow with vectors

This is done again with the quiver command.

    xP = linspace(xmin,xmax,21);
yP = linspace(ymin,ymax,21);

[xG, yG] = meshgrid(xP,yP);

for i=1:length(xP)
for j= 1:length(yP)
F = F_PortraitdePhase([xP(i),yP(j)]);
ux(j,i) = F(1);
uy(j,i) = F(2);
end
end

figure(1);
%   subplot(4,2,1:6)
hold on
quiver(xG, yG, ux, uy, 'Color', 'b');
title({['Phase portrait for ',typeproblem ],'Left-click to draw a trajectory ; right-click to stop'});
xlabel('x_1');ylabel('x_2');
axis([xmin xmax ymin ymax]);
hold on

%    if(plotEnergy==1)
%    figure(1);subplot(4,2,7:8);hold on;
%title('Amplitude vs. time');
%    xlabel('t');ylabel('(x_1^2+x_2^2)^{1/2}');
%    end

pause(1);

## Drawing of a few trajectories selected by clicking on the figure :

   figure(1);
%  subplot(4,2,1:6)
[xp,yp,button] = ginput(1)
while(button==1);
xinit = [xp ; yp];
dFdx1 = (F_PortraitdePhase([xp+eps,yp])-F_PortraitdePhase([xp-eps,yp]))/(2*eps);
dFdx2 = (F_PortraitdePhase([xp,yp+eps])-F_PortraitdePhase([xp,yp-eps]))/(2*eps);
disp([' Drawing trajectory starting from [x1,x2] = [ ',num2str(xp), ' , ' num2str(yp), ' ] ']);
disp([' Gradient matrix at starting point :   [ [ ' ,num2str(dFdx1(1)), ' , ',num2str(dFdx2(1)), ' ] ;']);
disp(['                                     [ ' ,num2str(dFdx1(2)), ' , ',num2str(dFdx2(2)), ' ] ]']);

[t,xtraj] = ode45(@(t,X)F_PortraitdePhase(X),linspace(0, Tmax,500),xinit);
[tbak,xtrajback] = ode45(@(t,X)F_PortraitdePhase(X),linspace(0,-Tmaxbak,500),xinit);
plot(xtraj(:,1),xtraj(:,2),'r',xtraj(1,1),xtraj(1,2),'ro');
hold on;
%    subplot(4,2,1:6);
plot(xtrajback(:,1),xtrajback(:,2),'m:');
axis([xmin xmax ymin ymax]);
xp=xtraj(end,1); yp=xtraj(end,2);
dFdx1 = (F_PortraitdePhase([xp+eps,yp])-F_PortraitdePhase([xp-eps,yp]))/(2*eps);
dFdx2 = (F_PortraitdePhase([xp,yp+eps])-F_PortraitdePhase([xp,yp-eps]))/(2*eps);
disp([' This trajectory ends at point [x1,x2] = [ ',num2str(xp), ' , ' num2str(yp), ' ] ']);
disp([' Gradient matrix at this point :         [ [ ' ,num2str(dFdx1(1)), ' , ',num2str(dFdx2(1)), ' ] ;']);
disp(['                                         [ ' ,num2str(dFdx1(2)), ' , ',num2str(dFdx2(2)), ' ] ]']);
disp(' ');
%    if(plotEnergy==1)
%figure(2);
%      subplot(4,2,7:8);
%      plot(t,sqrt(xtraj(:,1).^2+xtraj(:,2).^2));
%figure(1);
%     end
[xp,yp,button] = ginput(1);
end

figure(1);
xlabel('x_1');  ylabel('x_2');title('Phase portrait of a nonlinear 2x2 system');
set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','PhasePortraitNonLinear.png');


## Definition of the function F :

end

function F = F_PortraitdePhase(X)

% Fonction pour tracer le portrait de phase d'un systeme 2-2 de la forme
% d X /d T = F(X).
%
global typeproblem xmin xmax ymin ymax;
x1 = X(1); x2 = X(2);

switch (typeproblem)


### Case “Pendulum” :

Equation for the motion of a damped pendulum.

\displaystyle m L^2 \ddot \theta = - \mu_f \dot \theta - m g L \sin \theta

Considering this problem as a dynamical system of order 2 for X = [x_1 ; x_2] = [\theta, \dot \theta] leads to :

\displaystyle \frac{ d}{dt} [x_1 ; x_2] = [ x_2 ; - \mu x_2 - \omega_0^2 \sin(x_1) ]

with \omega_0 = \sqrt{\frac{g}{L}} ; \mu = \frac{\mu_f}{M L^2}.

Exercice : Observe the structure of the phase portrait for various values of the damping parameter mu. What do we observe for mu=0 ?

    case('Pendulum')

mu = .5;
omega0 = 1; % omega_0 = sqrt(g/L)
F = [x2;-omega0^2*sin(x1)-mu*x2];
xmin = -3*pi ; xmax = 3*pi; ymin = -2*pi; ymax = 2*pi;


### Case “RotatingPendulum” :

Equation for the motion of a pendulum with imposed precession

\displaystyle m L^2 \ddot \theta = - \mu_f \dot \theta - m g L \sin \theta + m L^2 \Omega^2 \sin \theta \cos \theta

Considering this problem as a dynamical system of order 2 for X = [x_1 ; x_2] = [\theta, \dot \theta] leads to :

\displaystyle \frac{ d}{dt} [x_1 ; x_2] = [ x_2 ; - \mu x_2 - \omega_0^2 \sin x_1 (1 - R \cos x_1 ) ]

with \omega_0 = \sqrt{\frac{g}{L}} ; \mu = \frac{\mu_f}{M L^2} ; R = \frac{\Omega^2}{\omega^2}.

Exercice : Observe the structure of the phase portrait for various values of the rotation parameter R. Draw the corresponding bifurcation diagram.

    case('RotatingPendulum')

omega0 = 1; mu = .5;
R = 1.0;
F = [x2;-omega0^2*sin(x1)*(1-R*cos(x1))-mu*x2];
xmin = -pi/2 ; xmax = pi/2; ymin = -.5; ymax =.5;


### Case “InvertedPendulum” :

Equation for the motion of an inverted pendulum with a spring (offset angle alpha)

\displaystyle m L^2 \ddot \theta = - \mu_f \dot \theta + m g L \sin \theta - K (\theta-\alpha)

Considering this problem as a dynamical system of order 2 for X = [x_1 ; x_2] = [\theta, \dot \theta] leads to :

\displaystyle \frac{ d}{dt} [x_1 ; x_2] = [ x_2 ; - \mu x_2 + \omega_0^2 \sin x_1 - k ( x_1- \alpha) ]

with \omega_0 = \sqrt{\frac{g}{L}} ; \mu = \frac{\mu_f}{m L^2} ; k = \frac{K}{m L^2}.

Exercice : Observe the structure of the phase portrait for various values of the offset angle alpha. Draw the corresponding bifurcation diagram.

    case('InvertedPendulum')

omega0 = 1; mu = .5;k = .9;
alpha = 0.03;

F = [x2;+omega0^2*sin(x1)-k*(x1-alpha)-mu*x2];
xmin = -pi/2 ; xmax = pi/2; ymin = -.5; ymax =.5;


### Case ‘VanDerPol’ :

The VanDerPol oscillator is a well-known model of self-sustained oscillator. It is defined as follows :

\displaystyle m \ddot x +\omega_0^2 x = (r - \delta x^2) \dot x

Exercice : Observe the structure of the phase portrait for various choices of the parameters. Reconstruct the corresponding bifurcation diagram in terms of the bifurcation parameter r.

    case('VanDerPol');
omega0 = 1;
r = 0.1;
deltahttp://basilisk.fr/_edit/sandbox/easystab/david/PhasePortrait_NonLinear.m = 1;
F = [x2 ; -omega0^2*x1+(r-delta*x1^2)*x2];
xmin = -3 ; xmax = 3; ymin = -3; ymax = 3;


### Case ‘Brusselator’ :

The ‘Brusselator’ (more precisely the homogeneous brusselator) is a simple model for a chemical reaction displaying unsteadiness. It is defined as follows :

\displaystyle \dot x_1 = -(\beta + 1) x_1 + x_1^2 x_2 + \alpha ; \displaystyle \dot x_2 = \beta x_1 - x_1^2 x_2 ;

Exercice : Observe the structure of the phase portrait for various choices of the parameters. Construct the bifurcation diagram as function of the bifurcation parameter beta.

   case('Brusselator')

alpha = 1;
beta = 2.1;
F = [-(beta+1)*x1+x1^2*x2+alpha ; beta*x1-x1^2*x2] ;

xmin = 0 ; xmax = 6; ymin = 0; ymax = 6;


### Case ‘LotkaVolterra’ :

The Lotka-Volterra is a simple system representing the competition of two animal species (for instance hares and lynxes). It is defined as follows :

\displaystyle \dot x_1 = x_1( \alpha - \beta x_2) ; \displaystyle \dot x_2 = x_2( -\gamma + \delta x_1) ;

Exercice : Observe the structure of the phase portrait for various choices of the parameters. Is this problem conservative or not ?

   case('LotkaVolterra')
% we take the same values as wikipedia...
alpha = 2/3;
beta = 4/3;
gamma = 1;
delta = 1;
F = [(alpha-beta*x2)*x1 ; (delta*x1-gamma)*x2];
xmin = 0 ; xmax = 2.5; ymin = 0; ymax = 2.5;


### Case ‘BuffaloWolf’ :

This problem is a

\displaystyle \dot x_1 = r x_1 - A x_1 (x_2+ E x_2^2) - B x_1^2; \displaystyle \dot x_2 = -C x_2 + D x_1 (x_2+ E x_2^2) ) ;

Exercice : Observe the structure of the phase portrait for various choices of the parameters. Reconstruct the corresponding bifurcation diagram as the parameter r is varied

    case('BuffaloWolf')
r = .4;
A = .3;
B = 0.1;
C = 1;
D = .2;
E = 0.15;
F = [(r*x1 -A*x1*(x2+E*x2^2) - B*x1^2) ; (-C*x2 + D*x1*(x2+E*x2^2)) ];

xmin = 0 ; xmax = 10; ymin = 0; ymax = 4;


### Case ‘Trefethen’ :

This is a simple model which is asymptotically stable but may lead to sucritical transition with very small initial conditions thanks to transient growth.

This problem is taken from Schmid & Henningson (Exercice 5.1, p. 525) and is actually a slight modification of the model of Trefethen et al. (science, 1993).

\displaystyle \dot x_1 = -(1/R) x_1 -\sqrt{x_1^2+x_2^2} x_2; \displaystyle \dot x_2 = -(2/R) x_2 + x_1 + \sqrt{x_1^2+x_2^2} x_1;

Exercice : observe the behaviour of the system for very small values of the initial condition (try with various values of R?)

    case('Trefethen')
R = 5;
NL = 0;% select 0 or 1
normX = sqrt(x1^2+x2^2)  ;
F1 = -(1/R)*x1 + NL*normX*x2;
F2 = -(2/R)*x2+x1 - NL*normX*x1;
F = [F1;F2];
xmin = -1 ; xmax = 1; ymin = -1; ymax = 1;


We investigate a saddle-node bifurcation in dimension 2 by considering the following model:

\displaystyle m \ddot x + \mu \dot x = r - x^2

If the friction is dominant over the inertia this equation reduces to the classical one-dimensional saddle-node bifurcation.

Exercice : Observe the structure of the phase portrait for various choices of the parameters. Reconstruct the corresponding bifurcation diagram.

   case('SaddleNode')
mass = .5; friction = 1; gravity = 1;
r= 1.5;
F = [x2 ; -friction*x2/mass + (r-x1^2)*gravity];
xmin = -3 ; xmax = 3; ymin = -1; ymax = 1;

### Case ‘Pitchfork’ :

We investigate a transcritical bifurcation in dimension 2 by considering the following model:

\displaystyle m \ddot x + \mu \dot x = r x - \delta x^3

If the friction is dominant over the inertia this equation reduces to the classical one-dimensional pitchfork bifurcation.

Exercice : Observe the structure of the phase portrait for various choices of the parameters. Reconstruct the corresponding bifurcation diagram.

   case('Pitchfork')
mass = .5; friction = 2;
r= .5; delta = 1;
F = [x2 ; -friction*x2/mass+(r*x1-delta*x1^3)/mass];
xmin = -1 ; xmax = 1; ymin = -1; ymax = 1;
%}


### Case ‘TransCritical’ :

We investigate a transcritical bifurcation in dimension 2 by considering the following model:

\displaystyle m \ddot x + \mu \dot x = r x - x^2

If the friction is dominant over the inertia this equation reduces to the classical one-dimensional transcritical bifurcation.

Exercice : Observe the structure of the phase portrait for various choices of the parameters. Reconstruct the corresponding bifurcation diagram.

    case('TransCritical')
mass = .5; friction = 1;
r= .5;
F = [x2 ; -friction*x2/mass+(r*x1-x1^2)/mass];
end%function