%% Nonlinear Dynamics coursework analysis tools (YFX1520)
%
% The coursework consists of two parts. The first part requires you to analyse a 2-D system and the second part requires the analysis of a 3-D system. Completed coursework may be handed over in two parts at any time during the semester.
%% Part 1: Analysis of a 2-D system
%  Perform linear analysis:
%  Find a fixed point or points of your system.
%  Linearise your system about the fixed point or points.
% Plot the linearized phase portrait/s (strong suggestion).
% Determine the type of the linearized fixed point or points.
% Determine if/how changes in the control parameter values influence the dynamics (type of fixed point/s) of your system.
% Perform nonlinear analysis of the full homogeneous system using a computer:
% Compare the type of the nonlinear fixed point or points with the corresponding linearized fixed point or points.
%  Plot the nonlinear phase portrait/s. Compare it/them with the corresponding linearized one/s.
%  Explain any discrepancies between the nonlinear and linearized systems if any occur.
%  Perform nonlinear analysis of the non-autonomous system (if applicable) using a computer:
%  Analyse the influence of explicitly time dependant part of the system on the system dynamics.
%  Study the nonlinear non-autonomous phase trajectories. Compare them against the corresponding homogeneous phase portrait.
%  Explain the obtained and presented results.
%  Draw overall conclusions and comment on the presented results
%% Note: In Matlab you can write functions a few different ways, however, here system functions are written locally (you can find them at the bottom of the code). When studying your systems you need to change the equations there. More about creating functions in Matlab can be found: function
% Also, a good reading tip is the ode solvers and how to select the best for your system: Choose ode solver
%% Linear analysis tool
% Use syms-command if you want to study using symbolic variables, if not write the needed constants. You can also define them straight at the local function below.

clear all; close all; clc

% Define symbolic variables or constans
syms a b real
% Constants
% a = 1;
% b = 1;

fun = @(x) func2D(0,x,a,b); % Calling the function using t = 0
syms x y
% Solving the f.p., jacobian, eigenvectors, eigenvalues, determinant and trace for the system
fixedpoints = solve(fun([x,y])==[0,0],[x,y]) 
jcb_mat = jacobian(fun([x,y]),[x,y])
% Jacobian with f.p.
x = fixedpoints.x;
y = fixedpoints.y;
% replaces all the variables in the symbolic expression
jcb_mat = subs(jcb_mat)
[eigenvector, eigenvalues] = eig(jcb_mat)
Determinant = det(jcb_mat)

%

Trace = trace(jcb_mat)

%
%% 2D Phase portrait plotting tool

% Constants 
a = 1;
b = 1;

% Initial conditions
t0 = [0.1 0.1];
% Time span 
tspan = [0 5];

% Solving the ode function
[t,x] = ode45(@(t,x) func2D(t,x,a, b), tspan, t0);
figure;
plot(t,x,'-') % Plotting the system
ylabel('Solutions')
xlabel('Time')
legend('x(t)','y(t)')
title('Phase plane for the system')
hold off
%% Creating the vector field for the system and plotting with phase potrait using 
% few different i.c.

% Generates two matrices one for all the x-values on the grid, and one for
% all the y-values on the grid. Note that x and y are matrices of the same
% size and shape, in this case 20 rows and 20 columns
y1 = linspace(min(x(:,1))-2,max(x(:,1))+2,20); 
y2 = linspace(min(x(:,2))-2,max(x(:,2))+2,20);
% linspace(X1, X2, N) generates N points between X1 and X2.
% Change the verctor fields boundaries (X1 and X2) if nessesary

% Computing the vector field
[xdot,ydot] = meshgrid(y1,y2);

% Generate two matrices, u and v
u = zeros(size(xdot));
v = zeros(size(xdot));

% Using a single loop over each element to compute the derivatives at
% each point (y1, y2)
t=0; % derivatives at each point at t=0, i.e. the starting time
for i = 1:numel(xdot)
    Yprime = func2D(t,[xdot(i); ydot(i)],a, b);
    u(i) = Yprime(1);
    v(i) = Yprime(2);
end
% quiver command to plot the vector field
quiver(xdot,ydot,u,v,'b'); figure(gcf)
title('Phase potrait with vector field')
xlabel('x(t)')
ylabel('y(t)')
axis tight equal;
hold on
% Plotting the f.p. 
fixedpoints= subs(fixedpoints);
plot(fixedpoints.x,fixedpoints.y, '.r','MarkerSize',20)
% Change the initial conditions to see how the trajactory changes
% Here both i.c. coordinates have the same value
for initial = [0 0.5 1 1.5]
    [ts,ys] = ode45(@(t,x) func2D(t,x,a, b),tspan,[initial;initial]); 
    plot(ys(:,1),ys(:,2))
    plot(ys(1,1),ys(1,2),'k*','MarkerSize',10) % starting point
    plot(ys(end,1),ys(end,2),'kx','MarkerSize',10) % ending point
end
hold off

%
%
%% Part 2: Analysis of a 3D system
% Compare the known properties of strange attractors agains your system
% Draw conclusion based on your analysis
%% 3D Phase trajectory plotting tool
%


syms x y z 
[x0,y0,z0] = solve(func3D(0,[x,y,z])==[0,0,0]) % Solve the f.p. for the system
% The time we want to calculate
tspan = [0 20];
% initial conditions for the first calculus
xin =-1;
yin =-1;
zin =-1.0001;
% initial conditions for the second calculus
xin1 =-1;
yin1 =-1;
zin1 =-1;

initial = [xin, yin, zin];          %column vector containing initial values
initial1 = [xin1, yin1, zin1];          %column vector containing initial values
% Calculating the ode functions
[t,xyz] = ode23(@func3D, tspan, initial);
[t1,xyz1] = ode23(@func3D, tspan, initial1); 

% Plotting the both solutions to the same figure to comparement, were the
% red dot is the fixed point of the system
figure
plot3(xyz(:,1),xyz(:,2), xyz(:,3));
hold on
plot3(xyz1(:,1),xyz1(:,2), xyz1(:,3));
plot3(x0,y0, z0, '.r','MarkerSize',20);
grid on
legend('1st', '2nd', 'f.p.')
xlabel('x(t)')
ylabel('y(t)')
zlabel('z(t)')
figure
subplot(1,3,1);
plot(xyz(:,1),xyz(:,3))
hold on
plot(xyz1(:,1),xyz1(:,3))
plot(x0,z0,'.r','MarkerSize',10)
xlabel('x(t)')
ylabel('z(t)')
subplot(1,3,2);
plot(xyz(:,1),xyz(:,2))
hold on
plot(xyz1(:,1),xyz1(:,2))
plot(x0,y0,'.r','MarkerSize',10)
xlabel('x(t)')
ylabel('y(t)')
subplot(1,3,3);
plot(xyz(:,2),xyz(:,3))
hold on
plot(xyz1(:,2),xyz1(:,3))
plot(y0,z0,'.r','MarkerSize',10)
legend('1st', '2nd', 'f.p.')
xlabel('y(t)')
ylabel('z(t)')
hold off

figure
subplot(3,1,1);
plot(t,xyz(:,1))
hold on
plot(t1,xyz1(:,1))
ylabel('x(t)')
subplot(3,1,2);
plot(t,xyz(:,2))
hold on
plot(t1,xyz1(:,2))
ylabel('y(t)')
legend('1st', '2nd')
subplot(3,1,3);
plot(t,xyz(:,3))
hold on
plot(t1,xyz1(:,3))
xlabel('time')
ylabel('z(t)')

%
%% Local fuctions (Documentation)
% Write here the systems you want to analyze.
% In example code the system for 2D is:
%
%     \dot x = -x + ay+x^2y \\
%     \dot y = -b-ay-x^2y
% and for the 3D:
%
%     \dot x = 10(-x+y)\\
%     \dot y = 28x-y-xz\\
% \dot z =xy-8/3z
%

function dxdt = func2D(t,x,a,b) % 2D system
  dxdt1 = -x(1)+a.*x(2)+x(1).^2.*x(2);
  dxdt2 = b-a.*x(2)-x(1).^2.*x(2);
  dxdt = [dxdt1,dxdt2]';
end

function dxdydz = func3D(t, xyz) % 3D system
    dxdydz1 = 10.*(-xyz(1)+xyz(2)); 
    dxdydz2 = 28.*xyz(1)-xyz(2)-xyz(1).*xyz(3);
    dxdydz3 = xyz(1).*xyz(2)-8/3.*xyz(3);
    dxdydz = [dxdydz1,dxdydz2,dxdydz3]';
end

% This code is provided by
% Maria Vuin <maria.vuin@taltech.ee> 
% Date: 05/06/2023 % 1st Matlab version