%Using MATLAB for math applications
%linear algebra
A=[5 -3 2; -3 8 4; 2 4 -9]; %initialize a matrix A
b=[10 20 9]'; %initialize column vector b
x=A\b %backslash solves the system
C=[A b] %concatenate to get augmented matrix
Cr=rref(C) %rref will automatically put matrix in reduced row echelon form
%can also write your own Gaussian elimination function, numerical analysis
%assignment once upon a time
%we can also use rref to find the inverse of a matrix
C=[A eye(size(A))] %concatenate A with identity matrix of the same size
Cr=rref(C)
A_inv = Cr(:,4:6) %keeps the last three columns as A_inv
A*A_inv == eye(3) %check to verify A_inv is correct
%machine epsilon is the cause of the discrepancy
A_inv2=inv(A) %built-in MATLAB inverse function
%----------------------------------------------
%eigenvalues and eigenvectors
eig(A) %outputs eigenvalues of A
[V,D]=eig(A) %V will contain eigenvectors, D will will be diagonal matrix with eigenvectors
%This satisfies A=V*D*inv(V) with machine epsilon error
%we can double check to make sure it is true
v2=V(:,2); %select second column
lambda=D(2,2); %select corresponding eigenvalue
A*v2-lambda*v2 %check the difference
%generalized eigenvalue problem: $Av=lambda Bv
%often studied in physical instances, such as coupled oscillators
%can be written as a matrix equation Mx_dd=-Kx
%note: this is equivalent to the standard eigenvalue problem if B is
%invertible, so let's do an example in which B is not invertible
B=[5 -3 1;2 9 -4;7 6 -3];
[V,D]=eig(A,B); %this satisfies A*V=B*V*D until an 'improper' eigenvalue is reached
%for large systems, you may only want the largest eigenvalues
eigs(A,2) %will return the two largest eigenvalues of A
%-----------------------------------------------
%matrix factorizations
%LU decomp
[L,U]=lu(A); %U upper triangular, L (almost) lower triangular such that LU=A
[L,U,P]=lu(A); %same as before, but P is the permutation matrix such that LU=PA
%QR decomp
[Q,R]=qr(A); %Q orthogonal, R upper triangular such that QR=A
%singular value decomp
[U,D,V]=svd(A); %U,V orthogonal matrices, D diagonal such that UDV=A
%---------------------------------------------------
%nonlinear algebraic equations
my_fun=@(x) cos(x)-x; %define a transcendental function for which you want to find zero
fzero(my_fun,0) %takes in function and initial guess for a zero
%the downside to fzero is that it can only find one zero (usually nearest
%initial guess)  for polynomial functions use roots instead
c=[1 0 -3 1 0 -9]; %declare vecotr of coeffiecients in descending order
roots(c) %MATLAB will also find complex zeros of a polynomial
%----------------------------------------------------------
%----------------------------------------------------------
%data analysis and regression
%basic fitting with MATLAB's built-in feature
x=[5 10 20 50 100]; y=[15 33 53 140 301]; %initialize data that is roughly linear
plot(x,y,'o')
%tools > basic fitting > show equations
%be mindful of "overfitting" the data with too high degree, show degree 4
%nonlinear fit
x=0:pi/50: pi;
y=sin(x)+rand(size(x))/10; %sin x with some noise
plot(x,y,'o')
%DYI polynomial curve fitting
%given n+1 data points, goal is to fit nth degree polynomial
a=polyfit(x,y,3) %last argument is degree of precision, so this will fit a cubic
%result is a vector containing coefficients
y2=polyval(a,x) %on the other hand, polyval will fit output data with the given coeffients
%exponential fitting
%if data appears to be exponential, convert to log scale and do a linear
%fit, then convert back
%let's construct y=P_0*exp(-x/Tau) fit
%so ln(y)=ln(P_0)*-1/Tau*x
x=[0 0.5 1 5 10 20]; y=[760 625 528 85 14 0.16]; %model exponential decay
ybar=log(y); %computes log of the output data
plot(x,ybar); %plot on log scale, looks linear
a=polyfit(x,ybar,1); %coefficients for a linear fit
p0=exp(a(2)); Tau=-1/a(1); %compute coefficient
xnew=linspace(0,20,20); %refines data
ynew=p0*exp(-xnew/Tau); %computes converted data points
figure(2), plot(x,y,'o',xnew,ynew)
%interpolation
x=[0 0.785 1.570 2.356 3.141 3.927 4.712 5.497 6.283];
y=[0 0.707 1.000 0.707 0.000 -.707 -1.00 -.707 -0.00]; %data that approximates sinx
xi=linspace(0,2*pi,50); %generate inputs to interpolate
yn=interp1(x,y,xi,'nearest'); %interpolates using nearest
subplot(2,2,1), plot(x,y,'o',xi,yn), title('nearest interp')
yi=interp1(x,y,xi,'linear'); %interpolates using linear (default if not stated)
subplot(2,2,2), plot(x,y,'o',xi,yi), title('linear interp')
yc=interp1(x,y,xi,'pchip'); %interpolates using cubics
subplot(2,2,3), plot(x,y,'o',xi,yc), title('cubic interp')
ys=spline(x,y,xi); %interpolates using cubic splines
subplot(2,2,4), plot(x,y,'o',xi,ys), title('cubic splines')
%CAUTION: for all interpolating functions, input data x needs to be
%monotonic
%can also do 2 and 3 variable interpolation with interp2,interp3

