% M21-126 Intro to Math Software
% Spring 2000
% Dr. David Handron

% lesson10.m (Introduction to Matlab)

		% This next command causes all the commands below to be shown
echo on		% echo off will only print out the answers.

% To use this file, start up matlab and issue the command
% >> lesson10

% To quit Matlab, issue the command >> quit


% To continue after a pause statment, hit any key to continue, or
% hit (control) C to terminiate the execution of a ".m" file like this.
pause

% ARRAY OPERATIONS

% The default arithmetic for arrays reverts to MATRIX operations; however,
% Matlab will do element by element (pointwise) operations on arrays of 
% the same size.
% To distinguish between matrix and pointwise multiplication, for example,
% matlab uses .* for the pointwise product (similarly ./ .^ etc.)
pause

x = 1:5
y = ones(1,5)

x + y 		% Element by element operation coincides with matrix operation
pause

x .* y		% Pointwise multiply
x .^ 2		% Pointwise square
sin(x)		% Pointwise sine
pause

x - 2		% Scalar subtraction
2*x		% Scalar multiplication
pause

% MATRIX OPERATIONS

x		% Recall x
x'		% The tick mark computes the transpose
pause

size(x)
size(x')
pause

x * y'		% (1x5)*(5x1) = 1x1
A = x' * y	% (5x1)*(1*5) = 5*5
pause

A * x'		% Matrix * Vector multiply
x * A		% Row vector * Matrix multiply
pause

A = hilb(5)	% A 5x5 Hilbert matrix
y = A \ x'	% This computes A^(-1) x'
pause

A * y - x'	% Check that A y = x'
pause

y = x / A	% Computes x A^(-1)
y * A - x	% Check y A = x
pause

A2  = A^2	% Compute A A
Am2 = A^(-2)	% Compute A^(-1) A^(-1)
pause

A2 * Am2	% Check
pause

rank(A)		% Rank of A
det(A)		% Determinant of A
pause

eig(A)		% Eigenvalues of A
pause

[X,D] = eig(A)	% (Generalized) eigenvctors and Jordan form of A
pause
X^(-1) * A * X	% Diagonalizing A
pause

% BASIC PLOTTING

x  = (1:5) *(pi/5)	% make x
y1 = sin(x)	% sin(x)
y2 = cos(x)	% cos(x)
pause

plot(y1)	% x-axis is 1,2,3,4,5
pause
plot([y1',y2'])	% Plot both of them - note they must be made into columns
pause

plot(x,y1)	% x-axis correct
pause
plot(x,[y1',y2'])	% Plot both
pause
plot(x,y1,x,y2)	% This is used when y2 has different x's
pause

% Plots can be embellished with title()'s, axes labels etc.
% Generally you make the graph first and then embellish it.

title('Sin and Cos Curves')
xlabel('x-axis')
ylabel('trig curves')
grid		% toggles the grid on and off
pause

% To make a two dimensional plot use the mesh() command

x = -3:3
y = x
pause

[X,Y] = meshgrid(x,y)	% This makes matricies of X and Y values
pause

x = -8:0.5:8	% Now for the real example
[X,Y] = meshgrid(x,x);	% These are big matricies
pause

% We want to plot sin(r) / r
R = (X.^2 + Y.^2).^(1/2);
Z = sin(R) ./ (R + 0.0000001);	% Don't divide by zero!
mesh(Z)
pause

contour(Z)	% All the contours are circular
axis equal	% Make x and y axes have equal scaling
pause
axis normal	% reset back to default rectangular box

% CASE STUDY: Deflection of a streched string

% A string under tension subjected to a force per unit length of f(x)
% will have deflection u(x) which satisfies u'' = -f. We will always
% assume that the string is tied at the ends (x=0 and x=1) so that
% u(0) = 0 and u(1) = 1. 

% The differential equation may be approximated by

%	   (-u(x-h) + 2 u(x) - u(x+h)) / h^2 = -f(x)	

% If we have this equation satisfied on a grid of uniformly spaced points
% we get a system of linear equations.
pause

N = 5
h = 1/(N-1)
x = ((1:N) - 1) * h	% Equal spaced x values
f = -5*sin(x*pi)	% Sinusoidal force per uinit length
pause

A = (2*diag(ones(1,N)) - diag(ones(1,N-1),1) - diag(ones(1,N-1),-1)) / h^2
pause

% The first equation should be changed to force u(0) = 0

A(1,1) = 1; A(1,2) = 0, f(1) = 0
pause

% And the last equation should be changed to force u(1) = 1

A(N,N) = 1; A(N,N-1) = 0, f(N) = 1
pause

u = A \ f'	% Solve the equations ...
plot(x,u)	% ... and plot the deflection
