function [ t, x, y, floq ] = periodic ( mu, y0 )
% Periodic
%
% This m-file is part of the material of the lectures on
%   Dynamical Systems and Nonlinear Analysis
% at the Math Department of Halle University, Germany.
%
% Author :      Martin Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Nov 30, 2009
%
% Periodic solutions of a parameter dependent dynamic system
%
% Parameters:
%   mu     (input)  : parameter value
%   y0     (input)  : initial guess y0 [optional]
%   t      (output) : independent variable of periodic solution
%   x      (output) : periodic solution, component x
%   y      (output) : periodic solution, component y
%   floq   (output) : vector of Floquet multipliers
%
% Example:
%   [ t, x, y, floq ] = periodic ( 1.5+0.01 );
%   plot ( x, y )
%
% Note:
%
% There is an unexpected behaviour for parameter values
%
% \mu = 1.538
% x0 = 0.0,   y0 = 2.293331357299218e-001
% x0 = 0.0,   y0 = 2.295615551878201e-001

% -> define initial values, rhs etc.
mu0 = 1.5;

x0 = 0.0;
if nargin<2
    y0 = 10 * ( mu - mu0 );
end

tolode = 1.0e-12;
tolnew = 1.0e-6;
options = odeset ( 'AbsTol', tolode, 'RelTol', tolode );

% -> initial time integration
tmax = 100;

rhs = @(t,u) ( perrhs ( t, u, mu ) );
[ tcur, ucur ] = ode45 ( rhs, [ 0 tmax ], [ x0 y0 ]', options );

sigvec = ucur(:,1)*ucur(2,1)<0;
[ sigaux ineg ] = max ( sigvec );
[ sigaux ipos ] = min ( sigvec(ineg:end) );

Tcur = mean ( tcur(ineg+ipos+(-2:-1)) );

% -> shooting method to get the periodic solution
itnew  = 1;
errnew = 0;

while ( ( errnew > tolnew + tolnew*norm([y0 Tcur]) ) || ( itnew==1 ) )

    y0n   = y0;
    Tcurn = Tcur;
    
    rhs = @(t,xvec) ( tausens ( t, xvec, mu, Tcurn ) );
    [ tcur, xvec ] = ode45 ( rhs, [ 0 1 ], [ x0 y0n 0 1 0 0]', options );

    f0  = ( xvec(end,1:2) - [ x0 y0 ] )';
    jac = [ xvec(end,3:4)'-[0;1] xvec(end,5:6)' ];

    diff = - jac \ f0;
    errnew = norm ( diff );
    
    gamma = 1.0;
    conv = 0;

    while ( conv==0 )

        y0n   = y0   + gamma * diff(1);
        Tcurn = Tcur + gamma * diff(2);

        rhs = @(t,u) ( taurhs ( t, u, mu, Tcurn ) );
        [ tcur, ucur ] = ode45 ( rhs, [ 0 1 ], [ x0 y0n ]', options );
        fn  = ( ucur(end,:) - [ x0 y0n ] )';

        if ( norm(jac\fn)<=(1-gamma/2)*errnew )
            conv = 1;
        else
            gamma = gamma / 2;
        end

    end
    
    y0    = y0n;
    Tcur  = Tcurn;
    itnew = itnew + 1;
    
end

rhs = @(t,xvec) ( persens ( t, xvec, mu ) );
[ t, xvec ] = ode45 ( rhs, [ 0 Tcur ], [ x0 y0 1 0 0 1 ]', options );
x = xvec(:,1);
y = xvec(:,2);

floq = eig ( [ xvec(end,3:4)' xvec(end,5:6)' ] );

end

% ----------------------------------------------------------------- %

function uprime = perrhs ( t, u, mu )
% PerRHS
%
% Right Hand Side of parameter dependent dynamic system with
% Periodic solution
%
% Parameters:
%   t      (input)  : current time (not used)
%   u      (input)  : current state u=(x,xp)
%   mu     (input)  : system parameter
%   uprime (output) : function value of right hand side

% -> read input data
x = u(1);
y = u(2);

% -> initialize output
uprime = zeros ( size(u) );

% -> evaluate right hand side
uprime(1) = mu * ( x - y ) - sin ( x*x + y*y );
uprime(2) = (2+mu)*x + (mu-3)*y - y*y*y;

end

% ----------------------------------------------------------------- %

function uprime = taurhs ( t, u, mu, T )
% TauRHS
%
% Right Hand Side of parameter dependent dynamic system with
% Periodic solution, Modified time scale
%
% Parameters:
%   t      (input)  : current time (not used)
%   u      (input)  : current state u=(x,xp)
%   mu     (input)  : system parameter
%   T      (input)  : scaling parameter of time scale \tau
%   uprime (output) : function value of right hand side

% -> output
uprime = T * perrhs ( t, u, mu );

end

% ----------------------------------------------------------------- %

function xvecp = persens ( t, xvec, mu )
% TauSens
%
% Right Hand Side of sensitivity equations
%
% Parameters:
%   t      (input)  : current time (not used)
%   xvec   (input)  : vector of unknowns [ u, du/dy0 du/dT ]
%   mu     (input)  : system parameter
%   xvecp  (output) : function value of right hand side

% -> read input data
u     = xvec(1:2);
dudx0 = xvec(3:4);
dudy0 = xvec(5:6);

x = u(1);
y = u(2);

% -> initialize output
xvecp = zeros ( size(xvec) );

% -> Jacobian of scaled right hand side
dfdu = [ mu-2*cos(x*x+y*y)*x  -mu-2*cos(x*x+y*y)*y
              2+mu                 mu-3-3*y*y      ];

% -> output
xvecp(1:2) = perrhs ( t, u, mu );
xvecp(3:4) = dfdu * dudx0;
xvecp(5:6) = dfdu * dudy0;

end

% ----------------------------------------------------------------- %

function xvecp = tausens ( t, xvec, mu, T )
% TauSens
%
% Right Hand Side of sensitivity equations, Modified time scale
%
% Parameters:
%   t      (input)  : current time (not used)
%   xvec   (input)  : vector of unknowns [ u, du/dy0 du/dT ]
%   mu     (input)  : system parameter
%   T      (input)  : scaling parameter of time scale \tau
%   xvecp  (output) : function value of right hand side

% -> read input data
u     = xvec(1:2);
dudy0 = xvec(3:4);
dudT  = xvec(5:6);

x = u(1);
y = u(2);

% -> initialize output
xvecp = zeros ( size(xvec) );

% -> Jacobian of scaled right hand side
dfdu = T * [ mu-2*cos(x*x+y*y)*x  -mu-2*cos(x*x+y*y)*y
                  2+mu                 mu-3-3*y*y      ];

% -> output
xvecp(1:2) = taurhs ( t, u, mu, T );
xvecp(3:4) = dfdu * dudy0;
xvecp(5:6) = dfdu * dudT + perrhs ( t, u, mu );

end

