% MATLAB seminar
% November 11, 1997
% R.W. Lessard
%
warning off
echo on
% Starting MATLAB
%    london> matlab
%    >>
%   
% MATLAB is a technical computing enviroment for high-performance numeric
% computation and visualization. MATLAB integrates numerical analysis,
% matrix computation, signal processing, and graphics.
% To see what MATLAB can do enter
%
%    >> demo
%
% at the command line. To learn the MATLAB basics enter
%
%    >> tour
%
% at the command line. To get help on MATLAB enter
%
%    >> helpdesk
%
% at the command line.
%
pause
clc
% MATRIX OPERATIONS
% The power of MATLAB is really rooted in its array and matrix processing.
% In fact MATLAB was originally written to provide easy access to matrix
% software developed by the LINPACK and EISPACK projects, which together
% represent the state of the art in software for matrix computation
% A MATLAB matrix can be entered in the form:
%
A = [1 2 3; 4 5 6; 7 8 9]
B = A'
pause
clc
%
% Matrix operations include addition, subtraction, multiplication and 
% division. High order operations such as powers (for a square matrix),
%
A^3/B^2
pause
clc
%
% and Trancendental and Elementary functions are easily applied
%
expm(A)
pause
clc
%
% ARRAY OPERATIONS
% Array operations refers to element-by-element arithmetic operations.
% Here .* ./ .^ operations refer to the element-by-element operation.
%
x = [1 2 3]
y = [4 5 6]
z = x.*y
pause
clc
%
% Mathematical operations can also be applied to arrays on an 
% element-by-element basis.
%
sin(pi*A/180)
erf(A)
pause
clc
%
% MATLAB derives much of its power from its matrix functions. Some
% are built into the MATLAB processor itself while others are available
% through a library of external script files. These classes of functions
% include:
%    triangular factorization
%    orthogonal factorization
%    singular value decomposition
%    eigenvalue decomposition
%
[L,U] = lu(A)
pause
clc
X = inv(A)
X = inv(U)*inv(L)
pause
clc
% GRAPHICS
% MATLAB has a large collection of plotting routines.
%
% Simple 1D Plots
%
t = 0:pi/100:2*pi;
x = sin(t);
y1 = sin(t+.25);
y2 = sin(t+.50);
plot(x,y1,'r-',x,y2,'g--')
grid;
title('Phase Shift');
xlabel('x=sin(t)');
ylabel('y=sin(t+)');
pause
plot(x);
hold on;
plot(y1,'--');
plot(y2,'-.');
xlabel('Number of Points');
title('Adding Plots');
grid;
hold off
pause
clc
%
% Matrix Plotting
%
plot(peaks);
pause
clc
%
% 3D Plotting
%
surf(peaks);
colormap(hot)
title('Mesh Drawing of Peaks');
pause
surfl(peaks);
pause
clc
%
% Images
% 
% NONLINEAR EQUATIONS AND OPTIMIZATION FUNCTIONS
% Solutions to nonlinear equations and optimization can be obtained
% from MATLAB using: fmin, fmins and fzero
%
% Curve Fitting a nonlinear equation to a set of data.
% FMINS performs a simplex method of minimizing a nonlinear
% function of several variables.
% Given some data:
%
global Data Plothandle
Data = [ 0.0000,   5.8955; 0.1000,  3.5639; 0.2000,   2.5173; ...
         0.3000,   1.9790; 0.4000,  1.8990; 0.5000,   1.3938; ...
         0.6000,   1.1359; 0.7000,  1.0096; 0.8000,   1.0343; ...
         0.9000,   0.8435; 1.0000,  0.6856; 1.1000,   0.6100; ...
         1.2000,   0.5392; 1.3000,  0.3946; 1.4000,   0.3903; ...
         1.5000,   0.5474; 1.6000,  0.3459; 1.7000,   0.1370; ...
         1.8000,   0.2211; 1.9000,  0.1704; 2.0000,   0.2636];
t = Data(:,1);
y = Data(:,2);
plot(t,y,'ro');
hold on;
Plothandle = plot(t,y,'-');
pause
clc
%
% We would like to fit the following function with 2 linear
% parameters and 2 nonlinear parameters to the data:
%
%    y =  c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t)
% 
% The equation is stored in the MATLAB function fitfun.
%
lam = [1 0]';
trace = 0;
tol = .1;
pause
lambda = fmins('fitfun',lam,[trace tol])
hold off
pause
clc
%
% SIMPLE PENDULUM
% Numerical Solutions of Ordinary Differential Equations
% The equation of motion for the Simple Pendulum can be expressed as
%
% thetadotdot + (g/l) * sin(theta) = 0	(1)
%
% This second order differential equation must be written as a system of
% coupled first order differential equations. Given the substitutions
%
% y1 = theta ; y2 = thetadot
%
% equation (1) is written as
%
% y1dot = y2
% y2dot = -(g/l) sin(y1)
%
g = 1;
l = 1;
[t,y] = ode23('simpen',[0 10],[1 0]);
plot(t,y);
legend('Position \theta','Velocity \theta\prime');
xlabel('Time t(s)');
pause
plot(y(:,1),y(:,2));
xlabel('\theta');
ylabel('\theta\prime');
title('Phase Plane');
%
% 
pause
clc
%
% Solutions of Definite Integrals (Numerical Quadrature)
% The solution for the Period as a function of initial amplitude, thetanot,
% can be obtained from the equation of conservation of energy T + U = E.
% From this we obtain the equation
%
%    thetadot = 2* sqrt(g/l)*sqrt(sin^2(thetanot/2) - sin^2(theta/2))
%
% which can be integrated to solve for the period tau.
%
thetanot = 0.5;
tau = 2*sqrt(l/g)*quad8('simpenT',0,thetanot-thetanot*0.00001,[],[],thetanot)
%
% or the equation may be placed into a MATLAB function and plotted with fplot
% for the range of thetanot from 0 to 1 radian.
%
fplot('simpenP',[0 1]);
xlabel('Initial Amplitude \theta_o');
ylabel('Period \tau (s)');
pause
clc
%
% CHAOTIC MOTION in 3-D
% The nonlinear differential equation, known as the Lorentz strange attractor
% can be written as
%
%   ydot = A*y
%
% with a vector valued function y(t) and a matrix A, which depend upon y:
%
%          [ -8/3  0  y(2) ]
%   A(y) = [  0   -10  10  ]
%          [-y(2)  28  -1  ]
%
A = [-8/3 0 0;0 -10 10;0 28 -1];
y = [35 -10 -7]';
h = 0.01;
fig = figure;
colordef(fig,'black');
p = plot3(y(1),y(2),y(3),'.r','EraseMode','none','MarkerSize',2);
axis([0 50 -25 25 -25 25]);
hold on;
echo off
for i=1:10000
   A(1,3) = y(2);
   A(3,1) = -y(2);
   ydot = A*y;
   y = y + h*ydot;
   set(p,'XData',y(1),'YData',y(2),'ZData',y(3));
   drawnow;
end
hold off
echo on

****************************************************************************
* simpen.m                                                                 *
****************************************************************************
*cut here*
function ydot = simpen(t,y)
%
% Simple Pendulum
%
% adotdot + (g/l) * sin(a) = 0	(1)
%
% This second order differential equation must be written as a system of
% coupled first order differential equations. Given the substitutions
%
% y1 = a ; y2 = adot
%
% equation (1) is written as
%
% y1dot = y2
% y2dot = -(g/l) sin(y1)
%
g = 9.81;
l = 1;
ydot = zeros(2,1);
ydot(1) = y(2);
ydot(2) = -(g/l)*sin(y(1));

****************************************************************************
* simpenP.m                                                                *
****************************************************************************
*cut here*
function y = simpenP(x)
%
y = 2*quad8('simpenT',0,x-x*0.00001,[],[],x);

****************************************************************************
* simpenT.m                                                                *
****************************************************************************
*cut here*
function y = simpenT(x,xo)
%
y = 1./sqrt(sin(xo/2)^2 - sin(x/2).^2);

****************************************************************************
* fitfun.m                                                                 *
****************************************************************************
*cut here*
function err = fitfun(lambda)
%   FITFUN(lambda) returns the error between the data and the
%   values computed by the current function of lambda.
%   FITFUN assumes a function of the form
%
%     y =  c(1)*exp(-lambda(1)*t) + ... + c(n)*exp(-lambda(n)*t)
%
%   with n linear parameters and n nonlinear parameters.

global Data Plothandle

t = Data(:,1);
y = Data(:,2);
A = zeros(length(t),length(lambda));
for j = 1:length(lambda)
   A(:,j) = exp(-lambda(j)*t);
end
c = A\y;
z = A*c;
set(Plothandle,'ydata',z)
drawnow
err = norm(z-y);