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 :
    % 'Pendulum', 'RotatingPendulum', 'InvertedPendulum','SaddleNode', 'Transcritical', 
    % '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');
    
    Figure 1 : Phase portrait of a 2x2 nonsystem. This example corresponds to the Van der Pol oscillator

    Figure 1 : Phase portrait of a 2x2 nonsystem. This example corresponds to the Van der Pol oscillator

    ## 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];
            % adjust the figure axes 
            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];
            % adjust the figure axes 
            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];
            % adjust the figure axes 
            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];
            % adjust the figure axes 
            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] ;
            
            
            % adjust the figure axes 
            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];
            % adjust the figure axes 
            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)) ];
    
        % adjust the figure axes 
     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];
     % adjust the figure axes 
     xmin = -1 ; xmax = 1; ymin = -1; ymax = 1;   
    

    Case ‘SaddleNode’ :

    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];
            % adjust the figure axes 
            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];
             % adjust the figure axes 
             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];
             % adjust the figure axes 
             xmin = -1 ; xmax = 1; ymin = -.5; ymax = .5;         
    %}
    
    
        case('Custom')
            % add your own case here !
        
        otherwise
            error(['Error in PhasePortrait_NonLinear : unknown type ', typeproblem]);
     
    end%swith
    
    end%function