%% Van der Pol Oscillator
%  File: vanderpol.m
%  Directory: 4_gyujtemegy/11_CCS/_1_ccs/ccs_2020/vanderpol
%  Author: Peter Polcz (ppolcz@gmail.com) 
%  
%  Created on 2020. October 08. (2020b)
%
%% Simulate the van der Pol model

fig = figure(1);
delete(fig.Children)

syms t x1 x2 real
x = [x1;x2];

%%% Van der Pol oscillator
% The nondimensional input/output equation of the system is:
% 
% $\ddot{y} - \mu(1 - y^2) \dot y + y = 0$.
% 
% Setting $x_1 = y$ and $x_2 = \dot y$, we obtain the following state-space
% model:
% 
% $\begin{aligned} &\dot{x}_1 = x_2 \\ &\dot{x}_2 = -x_1 + \mu(1 - x_1)^2 x_2 \end{aligned}$
% 
% Original van der Pol system, which has globally attractive a limit cycle.
% Every trajectory will tend to this limit cycle,
vdp_oscillator_sym = [
    x2
    (1 - x1^2)*x2 - x1
    ];
vdp_oscillator_fh = matlabFunction(vdp_oscillator_sym,'vars', {t x});

%%% Time inverted model
% If we compute the time inverted model, we obtain a locally asymptotically
% stable system, for which the domain of attraction is bounded by the
% original system's limit cycle.
% 
% Let $x(t) \leftarrow x(-t)$, then:
% 
% $\begin{aligned} &\dot{x}_1 = - x_2 \\ &\dot{x}_2 = x_1 - \mu(1 - x_1)^2 x_2 \end{aligned}$
% 
vdp_inverted_sym = [
    -x2
    -(1 - x1^2)*x2 + x1
    ];
vdp_inverted_fh = matlabFunction(vdp_inverted_sym,'vars', {t x});

%%%
% Simulate a trajectory of the original van der Pol oscillator, to obtain
% the limit cycle
simulate(vdp_oscillator_fh,[0 9],[2.01;0],'LineWidth',2,'Color',[0 0.5 0])

%%%
% Simulate some convergent trajectories of the time-inverted (locally
% asymptotically stable) van der Pol system
simulate(vdp_inverted_fh,[0 9],[
    0 1.5 1 -1.5 -1  1.2
    1 1   2  0   -2 -0.5 ])

%%%
% Simulate some diverging trajectories of the time-inverted system
simulate(vdp_inverted_fh,[0 2],[-2;-1],'b')

%%%
% Visualization
xlabel('x1(t)');
ylabel('x2(t)');

TxT = text(0.9,2.8,'Limit cycle');
TxT.HorizontalAlignment = 'center';
TxT.FontSize = 12;
TxT.FontWeight = 'bold';
TxT.Color = [0 0.5 0];

function simulate(fh,Tspan,ics,varargin)

    figure(1);
    hold on

    for x0 = ics
    
        [~,x] = ode45(fh,Tspan,x0);

        Pl = plot(x(:,1),x(:,2),varargin{:});
        plot(x(1,1),x(1,2),'.','MarkerSize',20,'Color',Pl.Color)

    end
        
    grid on
    box on
    axis equal
end