%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% matlab_notes.txt
% written by Amy Cassidy 09-08-2008
% Introduction to MATLAB
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Basics of MATLAB operations, matricies, eigenvectors and eigenvalues,
% complex numbers, symbolic math, plotting functions, plotting data, 
% solving ordinary differential equations symbolically defining 
% functions and basic programming structures.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Online Help:
%http://www.mathworks.com/access/helpdesk/help/techdoc/matlab.shtml

%A Tutorial:
%http://www.math.ufl.edu/help/matlab-tutorial/matlab-tutorial.html

help        %to get help
help solve  %to get help on solve
exit        %to exit

x=5         %define variable 'x' 
y=5;        %the semi-colon at the end of the expression supresses output
who         %see what variables are active:
clear x     %remove variable 'x'
Inf         %infinity

%another hint:  if you want to modify previous commands, use the up (and down) 
%arrows to access previous commands

%~~~~~~~~~~~~
%Operators
%~~~~~~~~~~~~
a+b   %add 
a-b   %subtract
a*b   %multiply
a^b   %raise to power

%other mathematical operations:
sqrt(x)  %square root of x
cos(x)   %cosine of x in radian
acos(x)  %arccosine of x
cosd(x)  %cosine of x in degrees
exp(x)   %exponential
log(x)   %natural logarithm of x
log10(x) %log base 10

%~~~~~~~~~~~~
%Matrices:
%~~~~~~~~~~~~
%Generating Matrices  
N=5
M=7
zeros(N,M) %NxM matrix of zeros (N rows and M columns):
ones(N)    %square matrix of ones
eye(N)     %NxN identity matrix
rand(N)    %generates a NxN matrix of random numbers on unit interval

%building a matrix by inputing data

A=[1 2 3; 4 5 6; 7 8 9] %generate a 3x3 matrix 
B=[1 2 3; 4 5 6; 7 8 9; 10 11 12] %generate a 4x3 matrix 
whos A   	%finding out information about 'A':
length(B) 	%length of longest dimension of B
ndims(B)   	%number of dimensions
size(B)    	%length of each dimension of B
numel(B)   	%total number of elements in B

%indexing matrices (in matlab, matrix indices run from 1 to N)
A(1,1)   	%returns the first row and first column of matrix A
A(n,m)    	%returns the 'n'th row and 'm'th column
A(:,m)    	%returns the 'm'th column
A(1:2,1:2) 	%returns a block containing the columns 1-2 or rows 1-2

det(A)  	%determinant of square matrix, A
inv(A)   	%inverse of matrix A
trace(A)  	%trace of 'A' (sum over diagonals)

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Matrix and Array operations:
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%array operations:  for all of the following, A and B must have the same
%size, unless one is a scalar.  Array operations are element by element.
A+B   	%add matrices (A and B must have the same size, unless one is a scalar)
A-B   	%subtract matrices (A and B must have the same size, unless one is a scalar)
C=A.*B  %array multiplication:  C(i,j)=A(i,j)*B(i,j) 
C=A./B  %array right division: C(i,j)=A(i,j)/B(i,j) 
C=A.\B  %array left division: C(i,j)=B(i,j)/A(i,j) 
C=A.^B  %array power: C(i,j)=A(i,j)*B(i,j).  Also C=power(A,B)

%matrix operations
A*B   	%matrix multiplication. Columns of A must equal rows of B, unless one is a scalar
C=A^B   %matrix A raised to the power B
A'      %complex conjugate transpose. Also ctranspose(A)
A.'     %transpose (non-complex conjugate).  Also transpose(A)

%two others that I'm not going to explain, but you can look up if you'd like
A/B  
A\B

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Eigenvalues and Eigenfunctions:
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
E=eig(A)  	%returns a vector containing the eigenvalues of square matrix, A
[V,D]=eig(A)  	%D is a diagonal matrix with the eigenvalues on the diagonal
              	%V is a matrix which contains the corresponding eigenvectors
           
%~~~~~~~~~~~~~~~~
%Complex Numbers
%~~~~~~~~~~~~~~~~~
%both 'i' and 'j' are equal to the sqrt(-1) and con be used for inputting
%complex numbers:
z=5+6j
%is the same as 
z=5+6i
%is the same as
z=5+6*sqrt(-1)

%for z=x+i*y = |z|exp(i*theta)
real(z)  %real part of 'z' = x
imag(z)  %imaginary part of 'z' = y
abs(z)   %absolute value of 'z' = sqrt(x^2+y^2)
angle(z) %phase angle of 'z' = theta
conj(z)  %complex conjugate of 'z;'

%~~~~~~~~~~~~~~~~
%Symbolic Math :
%~~~~~~~~~~~~~~~~

%creating symbolic variables:
u=sqrt(sym(2))
x=sym('x')     	%create symbolic variable x
syms x         	%equivalent to the previous statement
s=sym('sigma')
syms a b c d    %create symbolic variables, a,b,c,d

fx=expand((x+3)^2*(x-1))   %symbolic expansion: carry out multiplication of a symbolic 
			   %expression. Works on trig functions and exponentials as well.
factor(fx)   		   %factor the expresion fx

%symbolic expression:
f=sin(2*x)
g=a*x^2+b*x + 10   	%for a,b,x already defined as symbolic variables

%substitution
subs(f,pi/4)      	%substitutes pi/4 for x
subs(g,a,5)       	%substitutes 5 for a
subs(subs(g,a,5),b,-3)  %substitute 5 for 'a', then -3 for 'b'

simplify('cos(2*x)^2+sin(x)')   %simplify an expression
solve('y=x^2+3')    		%solves for 'y' in terms of 'x' (which has been defined
				%as a symbolic variable).

syms y
[x,y]=solve('x+y=3', 'x^2=y')  %solve system of two equations for x,y

%differentiation
diff(f)  	%differentiate function 'f'
diff(f,2)  	%2nd derivative of function 'f'
diff(g,x)  	%differentiate 'g' with respect to 'x'
diff(g,2,x)  	%differentiate 'g' twice with respect to 'x'

%integration
help sym/int  	%help for symbolic integration
int(f)         	%indefinite integral of 'f'
int(f,a,b)     	%indefinite integral of 'f' from a to b
int(g,x,c,d)    %indefinite integral of 'f' from c to d
int(g,x,0,1)    %indefinite integral of 'f' from 0 to 1


%~~~~~~~~~~~~~~~~~~~~~
%Plotting
%~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~
%Plotting Functions
%~~~~~~~~~~~~~~~~~~~~~
ezplot('f(x)') 	%plot function written as a string.  Default domain is -2pi, 2pi
ezplot('f(x)', [xmin, xmax])  %plot function with specified domain
fplot('f(x)',  [xmin, xmax])  %plot function with specified domain. Domain required.
ezsurf('f(x,y)', [xmin, xmax, ymin,ymax])	%3d surface plot
ezcontour('f(x,y)')				%contourplot

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Examples - plotting functions
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The following are equivalent, except for the ranges:
ezplot('sin(2*x)')
fplot('sin(2*x)', [0, pi])
y='sin(2*x)'
ezplot(y,[0, pi])
fplot(y,[0, pi])
fplot(@(x)sin(2*x), [0, pi])

%the following are equivalent
ezplot('x^2-2xy-y^2=1')
ezplot('x^2-2xy-y^2-1')

%note that ezcontour plots variables in alphabetic order (first goes on x-axis, second on y-axis).  Unless you specify.   See last expression.
ezcontour('x*sin(y)+x^2')	%contour plot
ezcontourf('x*sin(y)+x^2')	%filled contour plot
ezcontour(@(y,x)x*sin(y)+x^2)	%changes axes

ezsurf('x*sin(y)+x^2')		%surface plot
ezsurfc('x*sin(y)+x^2')		%surface + contour plot

%~~~~~~~~~~~~~~~~~~
% Plotting Data   
%~~~~~~~~~~~~~~~~~~
plot(x,f(x))		%plot data arrays x, f(x)
plot(x,f(x), y,g(y))    %plot data arrays x,f(x) and y,f(y)
semilogx(x,f(x))	%log10 scale for x-axis
semilogy(x,f(x))	%log10 scale for y-axis
loglog(x,f(x))		%plot f(x) vs. x on log-log axes
scatter(x,f(x))		%scatter plot
plot3(x,y,z)		%plots lines from columns of x,y,z in 3D

contour(z)		%contour plot of matrix 'z'
contour(x,y,z)		%countour plot of matrix 'z'
surface(x,y,z)		%surface plot of z=z(x,y) vs x,y

%some options:
xlabel('x label text')
ylabel('ylabel text')
legend('data set 1', 'data set 2')
axis([xmin xmax ymin ymax])	%control plotting range of current plot
ylim([ymin ymax]) 		%set y limits.  also xlim(), zlim()

%~~~~~~~~~~~~~~~~~~~~~~~~~
%Examples - plotting data
%~~~~~~~~~~~~~~~~~~~~~~~~~
xmin=-2*pi
xmax=2*pi
nSteps=100
x=linspace(xmin, xmax, nSteps)      %defines an array from xmin to xmax with nSteps

figure(1)		% where to plot the data
plot(x, sin(x))		% data to plot
ylabel('sin(x)')	% label for the y-axis
xlabel('x')		% label for the x-axis
title(['Plot of sin(x). nSteps=',num2str(nSteps)]); % title
text(-1, 1, "Text label at -1,1")   %generate a label 

%plotting multiple plots on a single graph
subplot(2,2,1)
plot(x,sin(x))

subplot(2,2,2)
ezplot('sin(x)')
title('Using ezplot')

subplot(2,2,3)
fplot('sin(x)', [-5, 5])
title('Using fplot')

subplot(2,2,4)
plot(x,sin(x), x,cos(x))
legend('sin(x)','cos(x)')
title('Two plots')

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Solving ODE's - Symbolically
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%First Order Differential Equation: solve dy/dt=f(y,t)
dsolve('Dy=f(y,t)','y(t0)=y0')		
%Second Order Differential Equation: solve d^2y/dt^2=f(Dy,y,t) 
dsolve('D2y=f(Dy,y,t)','y(t0)=y0', 'Dy(t0)=ydot0')
%System of two first order differential equations
dsolve('Dy1=f(y1,y2,t)', 'Dy2=f(y1,y2,t)', 'y1(t0)=y10', 'y2(t0)=y20')

%initial conditions:
%   y(t_0)=y0 specifies the value of y at t0
%   Dy(t_0)=ydot0 specifies the value of dy/dt at t0
%Note that 't' is the default independent variable

%some examples:
dsolve('Dx=-a*x')
dsolve('Dx=-a*x', 'x(0)=1)
dsolve('D2x+5*Dx=-6*x','x(0)=1','Dx(0)=0') %ans=3*exp(-2*t)-2*exp(-3*t)
S=dsolve('Dy1=t','Dy2=y1+t^2','y1(0)=0','y1(0)=1')
%to see the solutions:
S.y1
S.y2

%~~~~~~~~~~~~~~~~~~~~~
%Programming Basics
%~~~~~~~~~~~~~~~~~~~~~
%if-statement
if(x==0)
   'x is zero'
end

%if-else statement
if(x>0)
   'x is positive'
else
   'x is not positive'
end

%multiple conditions: && = AND, || = OR
if( a > 0 && b >0)
   'a and b are positive'
elseif ( (a > 0 && b < 0) || (a<0 && b> 0))
   'either a or b are greater than zero'
elseif ( a<=0 && b<=0)
   'a and b less than or equal to zero'
else
   'you messed something up'
end

%for loop [logistic map]
clear x;
r=3.4;
x(1)=0.25;
nSteps=100;
for i=1:nSteps
    x(i+1)=r*x(i)*(1-x(i);
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Functions and Scripts
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~
%Functions:
%~~~~~~~~~~~~
% user can define new functions in matlab.  The function must be in a file with name:
% functionname.m
% The first line of file must define the name of the new function. 
% Variables within the body of a function are local variables.
% Subfunctions can also be created in the same function definition file and
% can be used by that function, but not outside it.
% functions can be terminated by an END statement, but this is mostly optional

%Example:  File logistic.m
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Defines the logistic map.
function x=logistic(x0,r,nSteps)
clear x;
%nSteps=100;
x(1)=x0;
for i=1:nSteps
    x(i+1)=r*x(i)*(1-x(i));
end

%~~~~~~~~~~~~~~
%MATLAB Scripts
%~~~~~~~~~~~~~~
%A MATLAB script is a file that contains a sequence of matlab commands.  
%Scripts have file extension of ".m". Aka "M-files".

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Solving ODE's - Numerically
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%MATLAB has a number of ode solvers.  The best one to use depends on the problem.
%the recommended solver to start with is ode45.  Type 'help ode45' for more information
%and a list of other solvers

[T,Y]=ode45(fun, tspan, y0)

% input: 
%   fun = function to integrate
%   tspan = interval of integration
%   y0=initial conditions
% output:
%   T = time array (independent variable)
%   Y = array of solutions, corresponding to time 'T' (dependent variable)

%Example
%2D equations: Van der Pol:  d^2x/dt^2+mu*(x^2-1)dx/dt+x=0
%define the following function in a separate file
function dy = vdp(t,y,flag,a)
	dy=zeros(2,1);	%column vector
	%define the set of equations:
	dy(1)=y(2);
	dy(2)=a*(1-y(1)^2)*y(2)-y(1);
%----end of function
a=10;
tspan=[0 20]; %time range
y0=[2; 0];    %column vector of initial conditions
[T,Y]=ode45('vdp', tspan,y0,[],a);
 %-or-
function dy = vdp(t,y,a)
	dy=zeros(2,1);	%column vector
	%define the set of equations:
	dy(1)=y(2);
	dy(2)=a*(1-y(1)^2)*y(2)-y(1);
%----end of function
[T,Y]=ode45(@vdp, tspan,y0,[],a);

plot(t,y)  		%plot y1(t), y2(t)
plot(t,y(:,1), t,y(:,2))%same as above
plot(y(:,1),y(:,2))  	%phase plot:  y1(t) vs. y2(t)

%~~~~~~~~~~~~~~~~
%Exercises:
%~~~~~~~~~~~~~~~~
%(a) van der Pol equation: Plot x(t) as well as the phase space plot, dx/dt(t) vs. x(t) 
%    for different values of a and different intial conditions.
%(b) logistic map:  plot x vs n for different values of 'r' and initial conditions 0<=x0<=1
%(c) generate a scatter plot of the stable points of the logistic map for r [3.5,4]

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%One possible solution for (b) and (c):
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function x=logistic(x0,r,nSteps)
% iterate the logistic map: x_{n+1}=r[x_n(1-x_n)]
%    x0 = starting value 
%    nSteps = number of iterations

clear x;
%nSteps=100;
x(1)=x0;
for i=1:nSteps
    x(i+1)=r*x(i)*(1-x(i));
end
 
plot(x)
title(['Logistic Map. r=',num2str(r)])
xlabel('n')
ylabel('x_n')
xlim([0 nSteps])

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function [r,lo]=logistic_orbit(rmin,rmax,nrSteps)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculates the fixed points of the logistic map for 
% rmin <= r <= rmax
% and generates a scatter plot.
% Begin tracking values of x after 'nxStart' steps and
% continue until 'nxSteps'

x0=0.25;	%starting value of 'x'
nxStart=1000;	
nxSteps=1100;    
delta_r=(rmax-rmin)/(nrSteps-1);
r=[rmin:delta_r:rmax];
lo=zeros(nrSteps,nxSteps-nxStart+1);
 
for r_i=1:nrSteps
    xr=logistic(x0,r(r_i),nxSteps);
    lo(r_i,:)=xr(nxStart:nxSteps);
end
 
figure(1)
for nx=1:(nxSteps-nxStart-1)
    scatter(r,lo(:,nx),3,'b')
    hold on
end
hold off
title(['Logistic Map.'])
xlabel('n')
ylabel('x_n')
xlim([rmin rmax])