Matlab Notes

Overview of the MATLAB Environment:
MATLAB is a high-level technical computing language and interactive

environment for algorithm development, data visualization, data analysis, and
numeric computation. Using the MATLAB product, you can solve technical
computing problems faster than with traditional programming languages, such as
C, C++, and Fortran.
You can use MATLAB in a wide range of applications, including signal and image
processing, communications, control design, test and measurement, financial
modeling and analysis, and computational biology. Add-on toolboxes (collections
of special-purpose MATLAB functions, available separately) extend the MATLAB
environment to solve particular classes of problems in these application areas.
MATLAB provides a number of features for documenting and sharing your work.
You can integrate your MATLAB code with other languages and applications, and
distribute your MATLAB algorithms and applications. Features include:
High-level language for technical computing
Development environment for managing code, files, and data
Interactive tools for iterative exploration, design, and problem solving
Mathematical functions for linear algebra, statistics, Fourier analysis,
filtering, optimization, and numerical integration
2-D and 3-D graphics functions for visualizing data
Tools for building custom graphical user interfaces
Functions for integrating MATLAB based algorithms with external
applications and languages, such as C, C++, Fortran, Java, COM, and
Microsoft
Excel
The MATLAB System
The MATLAB system consists of these main parts:

Desktop Tools and Development Environment
This part of MATLAB is the set of tools and facilities that help you use and
become more productive with MATLAB functions and files. Many of these tools
are graphical user interfaces. It includes: the MATLAB desktop and Command
Window, an editor and debugger, a code analyzer, and browsers for viewing help,
the workspace, and folders.
Mathematical Function Library
This library is a vast collection of computational algorithms ranging from
elementary functions, like sum, sine, cosine, and complex arithmetic, to more
sophisticated functions like matrix inverse, matrix eigenvalues, Bessel functions,
and fast Fourier transforms.

The Language
The MATLAB language is a high-level matrix/array language with control flow
statements, functions, data structures, input/output, and object-oriented
programming features. It allows both "programming in the small" to rapidly create
quick programs you do not intend to reuse. You can also do "programming in the
large" to create complex application programs intended for reuse.
Graphics
MATLAB has extensive facilities for displaying vectors and matrices as graphs, as
well as annotating and printing these graphs. It includes high-level functions for
two-dimensional and three-dimensional data visualization, image processing,
animation, and presentation graphics. It also includes low-level functions that
allow you to fully customize the appearance of graphics as well as to build
complete graphical user interfaces on your MATLAB applications.
External Interfaces
The external interfaces library allows you to write C/C++ and Fortran programs
that interact with MATLAB. It includes facilities for calling routines from

MATLAB (dynamic linking), for calling MATLAB as a computational engine, and
for reading and writing MAT-files.
Starting and Quitting the MATLAB Program
Starting a MATLAB Session
On Microsoft Windows
platforms, start the MATLAB program by double-

clicking the MATLAB shortcut on your Windows desktop.
On Apple Macintosh
platforms, start MATLAB by double-clicking the MATLAB

icon in the Applications folder.
On UNIX
platforms, start MATLAB by typing matlab at the operating system

prompt.
When you start MATLAB, by default, MATLAB automatically loads all the
program files provided by MathWorks for MATLAB and other MathWorks
products. You do not have to start each product you want to use.
There are alternative ways to start MATLAB, and you can customize MATLAB
startup. For example, you can change the folder in which MATLAB starts or
automatically execute MATLAB statements upon startup.

The Desktop
When you start MATLAB, the desktop appears, containing tools (graphical user
interfaces) for managing files, variables, and applications associated with
MATLAB.
The following illustration shows the default desktop. You can customize the
arrangement of tools and documents to suit your needs.

Quitting the MATLAB Program
To end your MATLAB session, select File > Exit MATLAB in the desktop, or
type quit in the Command Window. You can run a script file named finish.m each
time MATLAB quits that, for example, executes functions to save the workspace.

Confirm Quitting
MATLAB can display a confirmation dialog box before quitting. To set this
option, select File > Preferences > General > Confirmation Dialogs, and select
the check box for Confirm before exiting MATLAB.

Basic Matrix Operations
This is a demonstration of some aspects of the MATLAB language.
First, let's create a simple vector with 9 elements called a.
a = [1 2 3 4 6 4 3 4 5]
a =

1 2 3 4 6 4 3 4 5

Now let's add 2 to each element of our vector, a, and store the result in a new
vector.
Notice how MATLAB requires no special handling of vector or matrix math.
b = a + 2
b =

3 4 5 6 8 6 5 6 7

Creating graphs in MATLAB is as easy as one command. Let's plot the result of
our vector addition with grid lines.
plot(b)
grid on

MATLAB can make other graph types as well, with axis labels.
bar(b)
xlabel('Sample #')
ylabel('Pounds')

MATLAB can use symbols in plots as well. Here is an example using stars to mark
the points. MATLAB offers a variety of other symbols and line types.
plot(b,'*')
axis([0 10 0 10])

One area in which MATLAB excels is matrix computation.
Creating a matrix is as easy as making a vector, using semicolons (;) to separate
the rows of a matrix.
A = [1 2 0; 2 5 -1; 4 10 -1]
A =

1 2 0
2 5 -1
4 10 -1

We can easily find the transpose of the matrix A.
B = A'
B =

1 2 4
2 5 10

0 -1 -1

Now let's multiply these two matrices together.
Note again that MATLAB doesn't require you to deal with matrices as a collection
of numbers. MATLAB knows when you are dealing with matrices and adjusts your
calculations accordingly.
C = A * B
C =

5 12 24
12 30 59
24 59 117

Instead of doing a matrix multiply, we can multiply the corresponding elements of
two matrices or vectors using the .* operator.
C = A .* B
C =

1 4 0
4 25 -10
0 -10 1

Let's find the inverse of a matrix ...
X = inv(A)
X =

5 2 -2
-2 -1 1
0 -2 1

... and then illustrate the fact that a matrix times its inverse is the identity matrix.
I = inv(A) * A

I =

1 0 0
0 1 0
0 0 1

MATLAB has functions for nearly every type of common matrix calculation.
There are functions to obtain eigenvalues ...
eig(A)
ans =

3.7321
0.2679
1.0000

... as well as the singular value decomposition.
svd(A)
ans =

12.3171
0.5149
0.1577

The "poly" function generates a vector containing the coefficients of the
characteristic polynomial.
The characteristic polynomial of a matrix A is

p = round(poly(A))
p =

1 -5 5 -1

We can easily find the roots of a polynomial using the roots function.
These are actually the eigenvalues of the original matrix.
roots(p)
ans =

3.7321
1.0000
0.2679

MATLAB has many applications beyond just matrix computation.
To convolve two vectors ...
q = conv(p,p)
q =

1 -10 35 -52 35 -10 1

... or convolve again and plot the result.
r = conv(p,q)
plot(r);
r =

1 -15 90 -278 480 -480 278 -90 15 -1

At any time, we can get a listing of the variables we have stored in memory using
the who or whos command.
whos
Name Size Bytes Class Attributes

A 3x3 72 double
B 3x3 72 double
C 3x3 72 double
I 3x3 72 double
X 3x3 72 double
a 1x9 72 double
ans 3x1 24 double
b 1x9 72 double
p 1x4 32 double
q 1x7 56 double
r 1x10 80 double

You can get the value of a particular variable by typing its name.

A
A =

1 2 0
2 5 -1
4 10 -1

You can have more than one statement on a single line by separating each
statement with commas or semicolons.
If you don't assign a variable to store the result of an operation, the result is stored
in a temporary variable called ans.
sqrt(-1)
ans =

0 + 1.0000i

Matrix Manipulation
This demo examines some basic matrix manipulations in MATLAB.
We start by creating a magic square and assigning it to the variable A.
A = magic(3)
A =

8 1 6
3 5 7
4 9 2

Here's how to add 2 to each element of A.
Note that MATLAB requires no special handling of matrix math.
A+2
ans =

10 3 8
5 7 9
6 11 4

The apostrophe symbol denotes the complex conjugate transpose of a matrix.
Here's how to take the transpose of A.
A'
ans =

8 3 4
1 5 9
6 7 2

The symbol * denotes multiplication of matrices.
Let's create a new matrix B and multiply A by B.
B = 2*ones(3)
A*B
B =

2 2 2
2 2 2
2 2 2

ans =

30 30 30
30 30 30
30 30 30

We can also multiply each element of A with its corresponding element of B by
using the .* operator.
A.*B

ans =

16 2 12
6 10 14
8 18 4

MATLAB has functions for nearly every type of common matrix calculation. For
example, we can find the eigenvalues of A using the "eig" command.
eig(A)
ans =

15.0000
4.8990
-4.8990

This concludes our brief tour of some MATLAB matrix handling capabilities.

Using FFT
This demonstration uses the FFT function to analyze the variations in sunspot
activity over the last 300 years.
Sunspot activity is cyclical, reaching a maximum about every 11 years. Let's
confirm that. Here is a plot of a quantity called the Zurich sunspot relative number,
which measures both number and size of sunspots. Astronomers have tabulated
this number for almost 300 years.
load sunspot.dat
year=sunspot(:,1);
relNums=sunspot(:,2);
plot(year,relNums)
title('Sunspot Data')

Here is a closer look at the first 50 years.
plot(year(1:50),relNums(1:50),'b.-');

The fundamental tool of signal processing is the FFT, or fast Finite Fourier
Transform. To take the FFT of the sunspot data type the following.
The first component of Y, Y(1), is simply the sum of the data, and can be removed.
Y = fft(relNums);
Y(1)=[];
A graph of the distribution of the Fourier coefficients (given by Y) in the complex
plane is pretty, but difficult to interpret. We need a more useful way of examining
the data in Y.
plot(Y,'ro')
title('Fourier Coefficients in the Complex Plane');
xlabel('Real Axis');
ylabel('Imaginary Axis');

The complex magnitude squared of Y is called the power, and a plot of power
versus frequency is a "periodogram".
n=length(Y);
power = abs(Y(1:floor(n/2))).^2;
nyquist = 1/2;
freq = (1:n/2)/(n/2)*nyquist;
plot(freq,power)
xlabel('cycles/year')
title('Periodogram')

The scale in cycles/year is somewhat inconvenient. We can plot in years/cycle and
estimate the length of one cycle.
plot(freq(1:40),power(1:40))
xlabel('cycles/year')

Now we plot power versus period for convenience (where period=1./freq). As
expected, there is a very prominent cycle with a length of about 11 years.
period=1./freq;
plot(period,power);
axis([0 40 0 2e+7]);
ylabel('Power');
xlabel('Period (Years/Cycle)');

Finally, we can fix the cycle length a little more precisely by picking out the
strongest frequency. The red dot locates this point.
hold on;
index=find(power==max(power));
mainPeriodStr=num2str(period(index));
plot(period(index),power(index),'r.', 'MarkerSize',25);
text(period(index)+2,power(index),['Period = ',mainPeriodStr]);
hold off;

FFT for Spectral Analysis
This example shows the use of the FFT function for spectral analysis. A common
use of FFT's is to find the frequency components of a signal buried in a noisy time
domain signal.
First create some data. Consider data sampled at 1000 Hz. Start by forming a time
axis for our data, running from t=0 until t=.25 in steps of 1 millisecond. Then form
a signal, x, containing sine waves at 50 Hz and 120 Hz.
t = 0:.001:.25;
x = sin(2*pi*50*t) + sin(2*pi*120*t);
Add some random noise with a standard deviation of 2 to produce a noisy signal y.
Take a look at this noisy signal y by plotting it.
y = x + 2*randn(size(t));
plot(y(1:50))
title('Noisy time domain signal')

Clearly, it is difficult to identify the frequency components from looking at this
signal; that's why spectral analysis is so popular.
Finding the discrete Fourier transform of the noisy signal y is easy; just take the
fast-Fourier transform (FFT).
Y = fft(y,251);
Compute the power spectral density, a measurement of the energy at various
frequencies, using the complex conjugate (CONJ). Form a frequency axis for the
first 127 points and use it to plot the result. (The remainder of the points are
symmetric.)
Pyy = Y.*conj(Y)/251;
f = 1000/251*(0:127);
plot(f,Pyy(1:128))
title('Power spectral density')
xlabel('Frequency (Hz)')

Zoom in and plot only up to 200 Hz. Notice the peaks at 50 Hz and 120 Hz. These
are the frequencies of the original signal.
plot(f(1:50),Pyy(1:50))
title('Power spectral density')
xlabel('Frequency (Hz)')

Optimal Fit of a Non-linear Function
This is a demonstration of the optimal fitting of a non-linear function to a set of
data. It uses FMINSEARCH, an implementation of the Nelder-Mead simplex
(direct search) algorithm, to minimize a nonlinear function of several variables.
First, create some sample data and plot it.
t = (0:.1:2)';
y = [5.8955 3.5639 2.5173 1.9790 1.8990 1.3938 1.1359 1.0096 1.0343 ...
0.8435 0.6856 0.6100 0.5392 0.3946 0.3903 0.5474 0.3459 0.1370 ...
0.2211 0.1704 0.2636]';
plot(t,y,'ro'); hold on; h = plot(t,y,'b'); hold off;
title('Input data'); ylim([0 6])

The goal is to fit the following function with two linear parameters and two
nonlinear parameters to the data:
y = C(1)*exp(-lambda(1)*t) + C(2)*exp(-lambda(2)*t)
To fit this function, we've create a function FITFUN. Given the nonlinear
parameter (lambda) and the data (t and y), FITFUN calculates the error in the fit
for this equation and updates the line (h).
type fitfun
function err = fitfun(lambda,t,y)
%FITFUN Used by FITDEMO.
% FITFUN(lambda,t,y) 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.

% Copyright 1984-2004 The MathWorks, Inc.
% $Revision: 5.8.4.1 $ $Date: 2004/11/29 23:30:50 $

A = zeros(length(t),length(lambda));
for j = 1:length(lambda)
A(:,j) = exp(-lambda(j)*t);
end
c = A\y;
z = A*c;
err = norm(z-y);

Make a guess for initial estimate of lambda (start) and invoke FMINSEARCH. It
minimizes the error returned from FITFUN by adjusting lambda. It returns the final
value of lambda. Use an output function to plot intermediate fits.
start = [1;0];
% We use an anonymous function to pass additional parameters t, y, h to the
% output function.
outputFcn = @(x,optimvalues,state) fitoutputfun(x,optimvalues,state,t,y,h);
options = optimset('OutputFcn',outputFcn,'TolX',0.1);
estimated_lambda = fminsearch(@(x)fitfun(x,t,y),start,options)
estimated_lambda =

3.5897
0.0030

Integer Arithmetic
This gives examples of performing arithmetic on signal and image integer data.
Load Integer Signal Data
Load measurement datasets comprising signals from four instruments using 8 and
16-bit A-to-D's resulting in data saved as int8, int16 and uint16. Time is stored as
uint16.
load integersignal

% Look at variables
whos Signal1 Signal2 Signal3 Signal4 Time1

Signal1 7550x1 7550 int8
Signal2 7550x1 7550 int8
Signal3 7550x1 15100 int16
Signal4 7550x1 15100 uint16

Time1 7550x1 15100 uint16

Plot Data
First we will plot two of the signals to see the signal ranges.
plot(Time1, Signal1, Time1, Signal2);
grid;
legend('Signal1','Signal2');

Here we see the values for int8. It is likely that these values would need to be
scaled to calculate the actual physical value that the signal represents e.g. Volts.
Process Data
We can perform standard arithmetic on integers such as +, -, *, and /. Let's say we
wished to find the sum of Signal1 and Signal2.
SumSig = Signal1 + Signal2; % Here we sum the integer signals.

Now let's plot the sum signal and see where it saturates.
cla;
plot(Time1, SumSig);
hold on;
Saturated = (SumSig == intmin('int8')) | (SumSig == intmax('int8')); % Find where
it has saturated
plot(Time1(Saturated),SumSig(Saturated),'rd');grid;
hold off;

The markers show where the signal has saturated.
Load Integer Image Data
Next we will look at arithmetic on some image data.
street1=imread('street1.jpg'); % Load image data
street2=imread('street2.jpg');
whos street1 street2

street1 480x640x3 921600 uint8
street2 480x640x3 921600 uint8

Here we see the images are 24-bit color, stored as three planes of uint8 data.
Display Images
Display first image.
cla;
image(street1); % Display image
axis equal; axis off

Display second image
image(street2); % Display image

Scale an Image
We can scale the image by a double precision constant but keep the image stored
as integers. For example,
duller = 0.5 * street2; % Scale image with a double constant but create an integer
whos duller

duller 480x640x3 921600 uint8

subplot(1,2,1);
image(street2);
axis off equal tight
title('Original'); % Display image

subplot(1,2,2);
image(duller);
axis off equal tight
title('Duller'); % Display image

Add the Images
We can add the two street images together and plot the ghostly result.
combined = street1 + duller; % Add |uint8| images
subplot(1,1,1)
cla;
image(combined); % Display image
title('Combined');

Single Precision Math
This gives some examples of performing arithmetic and linear algebra with single
precision data. It also shows an example where the results are computed
appropriately in single or double precision depending on the input.
create Double Precision Data
Let's first create some data, which is double precision by default.
Ad = [1 2 0; 2 5 -1; 4 10 -1]
Ad =

1 2 0

2 5 -1
4 10 -1

Convert to Single Precision
We can convert data to single precision with the single function.
A = single(Ad); % or A = cast(Ad,'single');
Create Single Precision Zeros and Ones
We can also create single precision zeros and ones with their respective functions.
n=1000;
Z=zeros(n,1,'single');
O=ones(n,1,'single');
Let's look at the variables in the workspace.
whos A Ad O Z n

A 3x3 36 single
Ad 3x3 72 double
O 1000x1 4000 single
Z 1000x1 4000 single
n 1x1 8 double

We can see that some of the variables are of type single and that the variable A (the
single precision version of Ad) takes half the number of bytes of memory to store
because singles require just four bytes (32-bits), whereas doubles require 8 bytes
(64-bits).
Arithmetic and Linear Algebra
We can perform standard arithmetic and linear algebra on singles.
B = A' % Matrix Transpose
B =

1 2 4
2 5 10
0 -1 -1

whos B

B 3x3 36 single

We see the result of this operation, B, is a single.
C = A * B % Matrix multiplication
C =

5 12 24
12 30 59
24 59 117

C = A .* B % Elementwise arithmetic
C =

1 4 0
4 25 -10
0 -10 1

X = inv(A) % Matrix inverse
X =

5 2 -2
-2 -1 1
0 -2 1

I = inv(A) * A % Confirm result is identity matrix
I =

1 0 0
0 1 0
0 0 1

I = A \ A % Better way to do matrix division than inv
I =

1 0 0
0 1 0
0 0 1

E = eig(A) % Eigenvalues
E =

3.7321
0.2679
1.0000

F = fft(A(:,1)) % FFT
F =

7.0000
-2.0000 + 1.7321i
-2.0000 - 1.7321i

S = svd(A) % Singular value decomposition
S =

12.3171
0.5149
0.1577

P = round(poly(A)) % The characteristic polynomial of a matrix
P =

1 -5 5 -1

R = roots(P) % Roots of a polynomial
R =

3.7321
1.0000
0.2679

Q = conv(P,P) % Convolve two vectors
R = conv(P,Q)
Q =

1 -10 35 -52 35 -10 1

R =

1 -15 90 -278 480 -480 278 -90 15 -1

stem(R); % Plot the result

A Program that Works for Either Single or Double Precision
Now let's look at a function to compute enough terms in the Fibonacci sequence so
the ratio is less than the correct machine epsilon (eps) for datatype single or
double.
% How many terms needed to get single precision results?
fibodemo('single')

% How many terms needed to get double precision results?
fibodemo('double')

% Now let's look at the working code.
type fibodemo

% Notice that we initialize several of our variables, |fcurrent|,
% |fnext|, and |goldenMean|, with values that are dependent on the
% input datatype, and the tolerance |tol| depends on that type as
% well. Single precision requires that we calculate fewer terms than
% the equivalent double precision calculation.
ans =

19

ans =

41

function nterms = fibodemo(dtype)
%FIBODEMO Used by SINGLEMATH demo.
% Calculate number of terms in Fibonacci sequence.

% $Revision: 1.1.4.1 $ $Date: 2004/03/22 23:53:55 $

fcurrent = ones(dtype);
fnext = fcurrent;
goldenMean = (ones(dtype)+sqrt(5))/2;
tol = eps(goldenMean);
nterms = 2;
while abs(fnext/fcurrent - goldenMean) >= tol
nterms = nterms + 1;
temp = fnext;
fnext = fnext + fcurrent;
fcurrent = temp;
end

Inverses of Matrices
This demo shows how to visualize matrices as images and uses this to illustrate
matrix inversion.
This is a graphic representation of a random matrix. The RAND command creates
the matrix, and the IMAGESC command plots an image of the matrix,
automatically scaling the color map appropriately.
n = 100;
a = rand(n);
imagesc(a);
colormap(hot);
axis square;

This is a representation of the inverse of that matrix. While the numbers in the
previous matrix were completely random, the elements in this matrix are anything

BUT random. In fact, each element in this matrix ("b") depends on every one of
the ten thousand elements in the previous matrix ("a").
b = inv(a);
imagesc(b);
axis square;

But how do we know for sure if this is really the correct inverse for the original
matrix? Multiply the two together and see if the result is correct, because just as
3*(1/3) = 1, so must a*inv(a) = I, the identity matrix. The identity matrix (almost
always designated by I) is like an enormous number one. It is completely made up
of zeros, except for ones running along the main diagonal.
This is the product of the matrix with its inverse: sure enough, it has the distinctive
look of the identity matrix, with a band of ones streaming down the main diagonal,
surrounded by a sea of zeros.
imagesc(a*b);
axis square;

Graphs and Matrices
This demo gives an explanation of the relationship between graphs and matrices
and a good application of SPARSE matrices.
A graph is a set of nodes with specified connections between them. An example is
the connectivity graph of the Buckminster Fuller geodesic dome (also a soccer ball
or a carbon-60 molecule).
In MATLAB, the graph of the geodesic dome can be generated with the BUCKY
function.
% Define the variables.
[B,V] = bucky;
H = sparse(60,60);
k = 31:60;
H(k,k) = B(k,k);

% Visualize the variables.
gplot(B-H,V,'b-');
hold on
gplot(H,V,'r-');
hold off
axis off equal

A graph can be represented by its adjacency matrix.
To construct the adjacency matrix, the nodes are numbered 1 to N. Then element
(i,j) of the matrix is set to 1 if node i is connected to node j, and 0 otherwise.
% Define a matrix A.
A = [0 1 1 0 ; 1 0 0 1 ; 1 0 0 1 ; 0 1 1 0];

% Draw a picture showing the connected nodes.
cla
subplot(1,2,1);
gplot(A,[0 1;1 1;0 0;1 0],'.-');
text([-0.2, 1.2 -0.2, 1.2],[1.2, 1.2, -.2, -.2],('1234')', ...
'HorizontalAlignment','center')

axis([-1 2 -1 2],'off')

% Draw a picture showing the adjacency matrix.
subplot(1,2,2);
xtemp=repmat(1:4,1,4);
ytemp=reshape(repmat(1:4,4,1),16,1)';
text(xtemp-.5,ytemp-.5,char('0'+A(:)),'HorizontalAlignment','center');
line([.25 0 0 .25 NaN 3.75 4 4 3.75],[0 0 4 4 NaN 0 0 4 4])
axis off tight

Here are the nodes in one hemisphere of the bucky ball, numbered polygon by
polygon.
subplot(1,1,1);
gplot(B(1:30,1:30),V(1:30,:),'b-');
for j = 1:30,
text(V(j,1),V(j,2),int2str(j),'FontSize',10);
end
axis off equal

To visualize the adjacency matrix of this hemisphere, we use the SPY function to
plot the silhouette of the nonzero elements.
Note that the matrix is symmetric, since if node i is connected to node j, then node
j is connected to node i.
spy(B(1:30,1:30))
title('spy(B(1:30,1:30))')

Now we extend our numbering scheme to the whole graph by reflecting the
numbering of one hemisphere into the other.
[B,V] = bucky;
H = sparse(60,60);
k = 31:60;
H(k,k) = B(k,k);
gplot(B-H,V,'b-');
hold on
gplot(H,V,'r-');
for j = 31:60
text(V(j,1),V(j,2),int2str(j), ...
'FontSize',10,'HorizontalAlignment','center');
end
hold off
axis off equal

Finally, here is a SPY plot of the final sparse matrix.
spy(B)
title('spy(B)')

In many useful graphs, each node is connected to only a few other nodes. As a
result, the adjacency matrices contain just a few nonzero entries per row.
This demo has shown one place where SPARSE matrices come in handy.
gplot(B-H,V,'b-');
axis off equal
hold on
gplot(H,V,'r-');
hold off

Sparse Matrices
This demonstration shows that reordering the rows and columns of a sparse matrix
S can affect the time and storage required for a matrix operation such as factoring
S into its Cholesky decomposition, S=L*L'.
Visualizing a Sparse Matrix
A SPY plot shows the nonzero elements in a matrix.
This spy plot shows a SPARSE symmetric positive definite matrix derived from a
portion of the Harwell-Boeing test matrix "west0479", a matrix describing
connections in a model of a diffraction column in a chemical plant.
load('west0479.mat')
A = west0479;
S = A * A' + speye(size(A));
pct = 100 / numel(A);

clf; spy(S), title('A Sparse Symmetric Matrix')
nz = nnz(S);
xlabel(sprintf('nonzeros=%d (%.3f%%)',nz,nz*pct));

Computing the Cholesky Factor
Now we compute the Cholesky factor L, where S=L*L'. Notice that L contains
MANY more nonzero elements than the unfactored S, because the computation of
the Cholesky factorization creates "fill-in" nonzeros. This slows down the
algorithm and increases storage cost.
tic, L = chol(S)'; t(1) = toc;
spy(L), title('Cholesky decomposition of S')
nc(1) = nnz(L);
xlabel(sprintf('nonzeros=%d (%.2f%%) time=%.2f sec',nc(1),nc(1)*pct,t(1)));

Reordering to Speed Up the Calculation
By reordering the rows and columns of a matrix, it may be possible to reduce the
amount of fill-in created by factorization, thereby reducing time and storage cost.
We will now try three different orderings supported by MATLAB.
reverse Cuthill-McKee
column count
minimum degree
Using the Reverse Cuthill-McKee
The SYMRCM command uses the reverse Cuthill-McKee reordering algorithm to
move all nonzero elements closer to the diagonal, reducing the "bandwidth" of the
original matrix.
p = symrcm(S);
spy(S(p,p)), title('S(p,p) after Cuthill-McKee ordering')
nz = nnz(S);


The fill-in produced by Cholesky factorization is confined to the band, so that
factorization of the reordered matrix takes less time and less storage.
tic, L = chol(S(p,p))'; t(2) = toc;
spy(L), title('chol(S(p,p)) after Cuthill-McKee ordering')
nc(2) = nnz(L);
xlabel(sprintf('nonzeros=%d (%.2f%%) time=%.2f sec', nc(2),nc(2)*pct,t(2)));

Using Column Count
The COLPERM command uses the column count reordering algorithm to move
rows and columns with higher nonzero count towards the end of the matrix.
q = colperm(S);
spy(S(q,q)), title('S(q,q) after column count ordering')
nz = nnz(S);

For this example, the column count ordering happens to reduce the time and
storage for Cholesky factorization, but this behavior cannot be expected in general.
tic, L = chol(S(q,q))'; t(3) = toc;
spy(L), title('chol(S(q,q)) after column count ordering')
nc(3) = nnz(L);

Using Minimum Degree
The SYMAMD command uses the approximate minimum degree algorithm (a
powerful graph-theoretic technique) to produce large blocks of zeros in the matrix.
r = symamd(S);
spy(S(r,r)), title('S(r,r) after minimum degree ordering')
nz = nnz(S);

The blocks of zeros produced by the minimum degree algorithm are preserved
during the Cholesky factorization. This can significantly reduce time and storage
costs.
tic, L = chol(S(r,r))'; t(4) = toc;
spy(L), title('chol(S(r,r)) after minimum degree ordering')
nc(4) = nnz(L);

Summarizing the Results
labels={'original','Cuthill-McKee','column count','min degree'};

subplot(2,1,1)
bar(nc*pct)
title('Nonzeros after Cholesky factorization')
ylabel('Percent');
set(gca,'xticklabel',labels)

subplot(2,1,2)
bar(t)
title('Time to complete Cholesky factorization')
ylabel('Seconds');
set(gca,'xticklabel',labels)

Graphical Representation of Sparse Matrices
This demo shows the finite element mesh for a NASA airfoil, including two
trailing flaps.
The data is stored in the file AIRFOIL.MAT. It holds 4253 pairs of (x,y)
coordinates of the mesh points. It also holds an array of 12,289 pairs of indices,
(i,j), specifying connections between the mesh points.
load airfoil
The Finite Element Mesh
First, scale x and y by 2^(-32) to bring them into the range [0,1]. Then form the
sparse adjacency matrix and make it positive definite.

% Scaling x and y
x = pow2(x,-32);
y = pow2(y,-32);

% Forming the sparse adjacency matrix and making it positive definite
n = max(max(i),max(j));
A = sparse(i,j,-1,n,n);
A = A + A';
d = abs(sum(A)) + 1;
A = A + diag(sparse(d));

% Plotting the finite element mesh
gplot(A,[x y])

Visualizing the Sparsity Pattern
SPY is used to visualize sparsity pattern. SPY(A) plots the sparsity pattern of the
matrix A.
spy(A)
title('The adjacency matrix.')

Symmetric Reordering - Reverse Cuthill-McKee
SYMRCM uses the Reverse Cuthill-McKee technique for reordering the adjacency
matrix. r = SYMRCM(A) returns a permutation vector r such that A(r,r) tends to
have its diagonal elements closer to the diagonal than A. This is a good preordering
for LU or Cholesky factorization of matrices that come from "long, skinny"
problems. It works for both symmetric and asymmetric A.
r = symrcm(A);
spy(A(r,r))
title('Reverse Cuthill-McKee')

Symmetric Reordering - COLPERM
Use j = COLPERM(A) to return a permutation vector that reorders the columns of
the sparse matrix A in non-decreasing order of non-zero count. This is sometimes
useful as a preordering for LU factorization: lu(A(:,j)).
j = colperm(A);
spy(A(j,j))
title('Column count reordering')

Symmetric Reordering - SYMAMD
SYMAMD gives a symmetric approximate minimum degree permutation. p =
SYMAMD(S), for a symmetric positive definite matrix A, returns the permutation
vector p such that S(p,p) tends to have a sparser Cholesky factor than S.
Sometimes SYMAMD works well for symmetric indefinite matrices too.
m = symamd(A);
spy(A(m,m))
title('Approximate minimum degree')

Matrix Exponentials
For background on the computation of matrix exponentials, see "Nineteen dubious
ways to compute the exponential of a matrix, twenty-five years later," SIAM
Review 45, 3-49, 2003.
Here are three of the 19 ways to compute the exponential of a matrix.

Start with the Matrix A
A = [0 1 2; 0.5 0 1; 2 1 0]
Asave = A;
A =

0 1.0000 2.0000
0.5000 0 1.0000
2.0000 1.0000 0

Scaling and Squaring
expmdemo1 is an implementation of algorithm 11.3.1 in Golub and Van Loan,
Matrix Computations, 3rd edition.
% Scale A by power of 2 so that its norm is < 1/2 .
[f,e] = log2(norm(A,'inf'));
s = max(0,e+1);
A = A/2^s;

% Pade approximation for exp(A)
X = A;
c = 1/2;
E = eye(size(A)) + c*A;
D = eye(size(A)) - c*A;
q = 6;
p = 1;
for k = 2:q
c = c * (q-k+1) / (k*(2*q-k+1));
X = A*X;
cX = c*X;
E = E + cX;
if p
D = D + cX;
else
D = D - cX;
end
p = ~p;
end
E = D\E;

% Undo scaling by repeated squaring
for k = 1:s
E = E*E;
end

E1 = E
E1 =

5.3091 4.0012 5.5778

2.8088 2.8845 3.1930
5.1737 4.0012 5.7132

Matrix Exponential via Taylor Series
expmdemo2 uses the classic definition for the matrix exponential. As a practical
numerical method, this is slow and inaccurate if norm(A) is too large.
A = Asave;

% Taylor series for exp(A)
E = zeros(size(A));
F = eye(size(A));
k = 1;
while norm(E+F-E,1) > 0
E = E + F;
F = A*F/k;
k = k+1;
end

E2 = E
E2 =

5.3091 4.0012 5.5778
2.8088 2.8845 3.1930
5.1737 4.0012 5.7132

Matrix Exponential via Eigenvalues and Eigenvectors
expmdemo3 assumes that the matrix has a full set of eigenvectors. As a practical
numerical method, the accuracy is determined by the condition of the eigenvector
matrix.
A = Asave;

[V,D] = eig(A);
E = V * diag(exp(diag(D))) / V;

E3 = E

E3 =

5.3091 4.0012 5.5778
2.8088 2.8845 3.1930
5.1737 4.0012 5.7132

Compare the Results
For this matrix, they all do equally well
E = expm(Asave);
err1 = E - E1
err2 = E - E2
err3 = E - E3
err1 =

1.0e-014 *

0.3553 0.1776 0
0.0444 0.0444 -0.0444
-0.0888 -0.0888 -0.3553

err2 =

1.0e-014 *

0 0 -0.2665
-0.0888 -0.0888 -0.0888
0.0888 -0.0888 0

err3 =

1.0e-013 *

-0.0711 -0.0444 -0.0888
-0.0799 -0.0666 -0.0844
-0.0622 -0.0533 -0.1155

Taylor Series Fails
Here is a matrix where the terms in the Taylor series become very large before they
go to zero. Consequently, expmdemo2 fails.
A = [-147 72; -192 93];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)
E1 =

-0.0996 0.0747
-0.1991 0.1494

E2 =

1.0e+006 *

-1.1985 -0.5908
-2.7438 -2.0442

E3 =

-0.0996 0.0747
-0.1991 0.1494

Eigenvalues and Vectors Fails
Here is a matrix that does not have a full set of eigenvectors. Consequently,
expmdemo3 fails.
A = [-1 1; 0 -1];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)
E1 =

0.3679 0.3679
0 0.3679

E2 =

0.3679 0.3679
0 0.3679

E3 =

0.3679 0
0 0.3679

Finite Difference Laplacian
This demo illustrates the computation and representation of the finite difference
Laplacian on an L-shaped domain.
The Domain
For this example, NUMGRID numbers points within an L-shaped domain. The
SPY function is a very useful tool for visualizing the pattern of non-zero elements
in a given matrix.
R = 'L'; % Other possible shapes include S,N,C,D,A,H,B
% Generate and display the grid.
n = 32;
G = numgrid(R,n);
spy(G)
title('A finite difference grid')
% Show a smaller version as sample.
g = numgrid(R,12)
g =

0 0 0 0 0 0 0 0 0 0 0 0
0 1 6 11 16 21 26 36 46 56 66 0
0 2 7 12 17 22 27 37 47 57 67 0
0 3 8 13 18 23 28 38 48 58 68 0
0 4 9 14 19 24 29 39 49 59 69 0
0 5 10 15 20 25 30 40 50 60 70 0
0 0 0 0 0 0 31 41 51 61 71 0
0 0 0 0 0 0 32 42 52 62 72 0
0 0 0 0 0 0 33 43 53 63 73 0
0 0 0 0 0 0 34 44 54 64 74 0
0 0 0 0 0 0 35 45 55 65 75 0
0 0 0 0 0 0 0 0 0 0 0 0

The Discrete Laplacian
Use DELSQ to generate the discrete Laplacian. The SPY function gives a
graphical feel of the population of the matrix.
D = delsq(G);
spy(D)

title('The 5-point Laplacian')
% Number of interior points
N = sum(G(:)>0)
N =

675

The Dirichlet Boundary Value Problem
Finally, we solve the Dirichlet boundary value problem for the sparse linear
system. The problem is setup as follows:
delsq(u) = 1 in the interior,
u = 0 on the boundary.
rhs = ones(N,1);
if (R == 'N') % For nested dissection, turn off minimum degree ordering.
spparms('autommd',0)
u = D\rhs;
spparms('autommd',1)
else

u = D\rhs; % This is used for R=='L' as in this example
end
The Solution
Map the solution onto the grid and show it as a contour map.
U = G;
U(G>0) = full(u(G(G>0)));
clabel(contour(U));
prism
axis square ij

Now show the solution as a mesh plot.
colormap((cool+1)/2);
mesh(U)
axis([0 n 0 n 0 max(max(U))])
axis square ij

Tessellation and Interpolation of Scattered Data
This demo describes convex hulls, Delaunay tessellations, and Voronoi diagrams
in 3 dimensions. It also shows how to interpolate three-dimensional scattered data.
Many applications in science, engineering, statistics, and mathematics use
structures like convex hulls, Delaunay tessellations and Voronoi diagrams for
analyzing data. MATLAB enables you to geometrically analyze data sets in any
dimension.
Here, in 3 dimensions, we show a set of 50 points with its convex hull.
% Create the data.
n = 50;
X = randn(n,3);
% Plot the points.
plot3(X(:,1),X(:,2),X(:,3),'ko','markerfacecolor','k');
% Compute the convex hull.

C = convhulln(X);
% Plot the convex hull.
hold on
for i = 1:size(C,1)
j = C(i,[1 2 3 1]);
patch(X(j,1),X(j,2),X(j,3),rand,'FaceAlpha',0.6);
end
% Modify the view.
view(3), axis equal off tight vis3d; camzoom(1.2)
colormap(spring)
rotate3d on

We can create a data set X of the 8 vertices of a cube plus its center.
X is a 9-by-3 matrix where each row is the 3-D coordinates of one point.
% Create X.
X = zeros(8,3);
X([5:8,11,12,15,16,18,20,22,24]) = 1; % Corners.
X(9,:) = [0.5 0.5 0.5]; % Center.
% Visualize X.

cla reset; hold on
d = [1 2 4 3 1 5 6 8 7 5 6 2 4 8 7 3];
plot3(X(d,1),X(d,2),X(d,3),'b:');
plot3(X(:,1),X(:,2),X(:,3),'b.','markersize',20);
t = text(X(:,1),X(:,2),X(:,3), num2str((1:9)'));
set(t,'VerticalAlignment','bottom','FontWeight','bold', 'FontSize',12);
view(3); axis equal tight off vis3d; camorbit(10,0);
rotate3d on

The convex hull of a data set is the smallest convex region that contains the data
set. The convex hull of the cube data set X can be computed by CONVHULLN.
For this data set X, the convex hull has 12 facets, each corresponding to a row in K
and plotted above with X. The cube is transparent so that you can see all the facets
and data points.
% Compute the convex hull.
tri = convhulln(X);
% Plot the data
cla reset;

% Plot the convex hull.
for i = 1:size(tri,1)
c = tri(i,[1 2 3 1]);
patch(X(c,1),X(c,2),X(c,3),i,'FaceAlpha', 0.9);
end
% Modify the view.
view(3); axis equal tight off vis3d
rotate3d on

A Delaunay tessellation in 3 dimensions is a set of tetrahedrons such that no data
points are contained in any tetrahedron's circumsphere. The Delaunay tessellation
of the data set X can be computed by DELAUNAYN.
The 12 rows of T represent the 12 tetrahedrons that partition the data set X.
% Compute the delaunay tessellation.
tri = delaunayn(X);
% Plot the data.
% Plot the tessellation.

y = tri(i,[1 1 1 2; 2 2 3 3; 3 4 4 4]);
x1 = reshape(X(y,1),3,4);
patch(x1,x2,x3,(1:4)*i,'FaceAlpha',0.8);
end
% Modify the view.
view(3); axis equal tight off vis3d; camorbit(10,0)
rotate3d on

A Voronoi diagram partitions the data space into polyhedral regions, with one
region for each data point. Anywhere within a region is closer to its data point than
any other in the set. The Voronoi diagram of the cube data set X can be computed
by VORONOIN.
V is the set of Voronoi vertices. C represents the set of Voronoi regions. For our
data set X, C has 9 Voronoi regions. Here we show one Voronoi region, the region
for the center point of the cube.
% Compute Voronoi diagram.
[c,v] = voronoin(X);

% Plot the data.
plot3(X(d,1),X(d,2),X(d,3),'b:.',X(9,1),X(9,2),X(9,3),'k.','markersize',20);
% Plot the Voronoi diagram.
nx = c(v{9},:);
tri = convhulln(nx);
patch(nx(tri(i,:),1),nx(tri(i,:),2),nx(tri(i,:),3),rand,'FaceAlpha',0.8);
end
% Modify the view.
view(3); axis equal tight off vis3d; camzoom(1.5); camorbit(20,0)
rotate3d on

GRIDDATAN interpolates multidimensional scattered data. It uses
DELAUNAYN to tessellate the data, and then interpolates based on the
tessellation. We start with a data set of 500 random points in 3 dimensions and
compute the values of a function, the squared distance from the origin, at each of
these points.
% Create the data.
n = 500;

X = 2*rand(n,3)-1;
v = sum(X.^2,2);
% Draw a picture to show how X is defined.
cla reset; hold on
plot3([0.02, 0.47],[0.02,0.57],[0.02,0.57],'k-','linewidth',2);
plot3(0,0,0,'bo','markerfacecolor','b');
cube = zeros(8,3);
cube([5:8,11,12,15,16,18,20,22,24]) = 1; % Corners
cube(9,:) = [0.5 0.5 0.5]; % Center.
plot3(cube(d,1),cube(d,2),cube(d,3),'r:');
plot3(0.5,0.6,0.6,'go','markerfacecolor','g');
text(0.02,-0.2,0,'the origin','fontsize',12,'fontweight','bold');
text(0.55,0.6,0.6,'a point in X','fontsize',12,'fontweight','bold');
text(0.28,0.3,0.35,'sqrt(v)','fontsize',12,'fontweight','bold');
view(3); axis equal tight off vis3d;
camorbit(20,-10);
rotate3d on

With GRIDDATAN we can interpolate X and the values v over a grid X0 to get
the function values v0 over this grid.

The black points are X and the red points are the grid X0.
% Grid the data.
d = -0.8:0.2:0.8;
[x0,y0,z0] = meshgrid(d,d,d);
X0 = [x0(:) y0(:) z0(:)];
v0 = reshape(griddatan(X,v,X0),size(x0));
% Plot results.
cla reset; hold on;
plot3(X(:,1),X(:,2),X(:,3),'k+','markerfacecolor','k');
plot3(X0(:,1),X0(:,2),X0(:,3),'r.','markerfacecolor','r');
view(3); axis equal tight off vis3d; camzoom(1.6);
rotate3d on

To visualize the surface at all points where the function takes on a constant value,
we can use ISOSURFACE and ISONORMALS. Since the function is the squared
distance from the origin, the surface at a constant value is a sphere.
c = 0.6; % constant value
cla reset; hold on;
h = patch(isosurface(x0,y0,z0,v0,c),'FaceColor','red','EdgeColor','none');

isonormals(x0,y0,z0,v0,h);
view(3); axis equal tight off vis3d; camzoom(1.6); camlight; lighting phong
rotate3d on

With more data points in X and a denser grid X0, the sphere is smoother but takes
longer to compute.
Here is a precomputed sphere generated using 5000 data points in X and a distance
between gridpoints of 0.05.
% Load saved results.
load qhulldemo
cla reset; hold on
d = -0.8:0.05:0.8;
[x0,y0,z0] = meshgrid(d,d,d);
h = patch(isosurface(x0,y0,z0,v0,0.6));
isonormals(x0,y0,z0,v0,h);
set(h,'FaceColor','red','EdgeColor','none');
view(3); axis equal tight off vis3d; camzoom(1.6)
camlight; lighting phong
rotate3d on

Creating and Editing Delaunay Triangulations
The Delaunay triangulation is the most widely used triangulation in scientific
computing. The properties associated with the triangulation provide a basis for
solving a variety of geometric problems. The following examples demonstrate how
to create, edit, and query Delaunay triangulations using the DelaunayTri class.
Construction of constrained Delaunay triangulations is also demonstrated, together
with an applications covering medial axis computation and mesh morphing.

Example One: Create and Plot a 2D Delaunay Triangulation
This example shows you how to compute a 2D Delaunay triangulation and how to
plot the triangulation together with the vertex and triangle labels.
x = rand(10,1);

y = rand(10,1);
dt = DelaunayTri(x,y)
dt =

DelaunayTri

Properties:
Constraints: []
X: [10x2 double]
Triangulation: [11x3 double]

triplot(dt);
%
% Display the Vertex and Triangle labels on the plot
hold on
vxlabels = arrayfun(@(n) {sprintf('P%d', n)}, (1:10)');
Hpl = text(x, y, vxlabels, 'FontWeight', 'bold', 'HorizontalAlignment',...
'center', 'BackgroundColor', 'none');
ic = incenters(dt);
numtri = size(dt,1);
trilabels = arrayfun(@(x) {sprintf('T%d', x)}, (1:numtri)');
Htl = text(ic(:,1), ic(:,2), trilabels, 'FontWeight', 'bold', ...
'HorizontalAlignment', 'center', 'Color', 'blue');
hold off

Example Two: Create and Plot a 3D Delaunay Triangulation
This example shows you how to compute a 3D Delaunay triangulation and how to
plot the triangulation.
X = rand(10,3)
X =

0.7509 0.1978 0.4959
0.8229 0.4084 0.2697
0.2251 0.3626 0.6842
0.3867 0.1619 0.7480
0.4187 0.8935 0.3916
0.9358 0.4405 0.5588
0.1056 0.6703 0.7322
0.5313 0.2104 0.6116
0.6212 0.4711 0.5812
0.2803 0.0499 0.0706

dt = DelaunayTri(X)

dt =

DelaunayTri

Properties:
Constraints: []
X: [10x3 double]

tetramesh(dt, 'FaceColor', 'cyan');
% To display large tetrahedral meshes use the convexHull method to
% compute the boundary triangulation and plot it using trisurf.
% For example;
% triboundary = convexHull(dt)
% trisurf(triboundary, X(:,1), X(:,2), X(:,3), 'FaceColor', 'cyan')

Example Three: Access the Triangulation Data Structure
There are two ways to access the triangulation data structure. One way is via the
Triangulation property, the other way is using indexing.

Create a 2D Delaunay triangulation from 10 random points.
X = rand(10,2)
X =

0.3742 0.5123
0.0336 0.4076
0.4061 0.2198
0.6580 0.2305
0.5752 0.4385
0.9760 0.9199
0.4416 0.3190
0.7778 0.2843
0.5695 0.3638
0.0211 0.4993

dt = DelaunayTri(X)
dt =

DelaunayTri

Properties:
Constraints: []
X: [10x2 double]

% The triangulation datastructure is;
dt.Triangulation
ans =

1 6 10
2 7 1
2 3 7
2 1 10
7 5 1
7 3 4
5 6 1
5 7 9
8 5 9

4 8 9
9 7 4
6 5 8

% Indexing is a shorthand way to query the triangulation. The format is
% dt(i, j) where j is the j'th vertex of the i'th triangle, standard
% indexing rules apply.
% The triangulation datastructure is
dt(:,:)
ans =

1 6 10
2 7 1
2 3 7
2 1 10
7 5 1
7 3 4
5 6 1
5 7 9
8 5 9
4 8 9
9 7 4
6 5 8

The second triangle is;
dt(2,:)
ans =

2 7 1

The third vertex of the second triangle is;
dt(2,3)
ans =

1

The first three triangles;
dt(1:3,:)
ans =

1 6 10
2 7 1
2 3 7

Example Four: Edit a Delaunay Triangulation to Insert or Remove Points
This example shows you how to use index-based subscripting to insert or remove
points. It is more efficient to edit a DelaunayTri to make minor modifications as
opposed to recreating a new DelaunayTri from scratch, this is especially true if the
dataset is large.
% Construct a Delaunay triangulation from
% 10 random points within a unit square
x = rand(10,1);
y = rand(10,1);
dt = DelaunayTri(x,y)
dt =

DelaunayTri

Properties:
Constraints: []
X: [10x2 double]

% Insert 5 additional random points
dt.X(end+(1:5),:) = rand(5,2)
dt =

DelaunayTri

Properties:
Constraints: []
X: [15x2 double]


Replace the fifth point
dt.X(5,:) = [0, 0]
dt =

DelaunayTri

Properties:
Constraints: []
X: [15x2 double]

Remove the fourth point
dt.X(4,:) = []
dt =

DelaunayTri

Properties:
Constraints: []
X: [14x2 double]

Example Five: Create a Constrained Delaunay Triangulation
This example shows you how to create a simple constrained Delaunay
triangulation and illustrates the effect of the constraints.
X = [0 0; 16 0; 16 2; 2 2; 2 3; 8 3; 8 5; 0 5];
C = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1];
dt = DelaunayTri(X, C);
subplot(2,1,1);

triplot(dt);
axis([-1 17 -1 6]);
xlabel('Constrained Delaunay triangulation', 'fontweight','b');
% Plot the constrained edges in red
hold on;
plot(X(C'),X(C'+size(X,1)),'-r', 'LineWidth', 2);
hold off;

% Now delete the constraints and plot the unconstrained Delaunay
dt.Constraints = [];
subplot(2,1,2);
triplot(dt);
axis([-1 17 -1 6]);
xlabel('Unconstrained Delaunay triangulation', 'fontweight','b');

Example Six: Create a Constrained Delaunay Triangulation of a
Geographical Map
Load a map of the perimeter of the conterminous United States. Construct a
constrained Delaunay triangulation representing the polygon. This triangulation
spans a domain that is bounded by the convex hull of the set of points. Filter out

the triangles that are within the domain of the polygon and plot them. Note: The
dataset contains duplicate datapoints; that is two or more datapoints have the same
location. The duplicate points are rejected and the DelaunayTri reformats the
constraints accordingly.
clf
load usapolygon
% Define an edge constraint between two successive
% points that make up the polygonal boundary.
nump = numel(uslon);
C = [(1:(nump-1))' (2:nump)'; nump 1];
dt = DelaunayTri(uslon, uslat, C);
io = dt.inOutStatus();
patch('faces',dt(io,:), 'vertices', dt.X, 'FaceColor','r');
axis equal;
axis([-130 -60 20 55]);
xlabel('Constrained Delaunay Triangulation of usapolygon', 'fontweight','b');
Warning: Duplicate data points have been detected and removed.
The Triangulation indices and Constraints are defined with
respect to the unique set of points in DelaunayTri property X.
Warning: Intersecting edge constraints have been split, this may have added
new points into the triangulation.

Example Seven: Curve Reconstruction from a Point Cloud
This example highlights the use of a Delaunay triangulation to reconstruct a
polygonal boundary from a cloud of points. The reconstruction is based on the
elegant Crust algorithm.
Reference: N. Amenta, M. Bern, and D. Eppstein. The crust and the beta-skeleton:
combinatorial curve reconstruction. Graphical Models and Image Processing,
60:125-135, 1998.
% Create a set of points representing the point cloud
numpts=192;
t = linspace( -pi, pi, numpts+1 )';
t(end) = [];
r = 0.1 + 5*sqrt( cos( 6*t ).^2 + (0.7).^2 );
x = r.*cos(t);
y = r.*sin(t);
ri = randperm(numpts);
x = x(ri);
y = y(ri);

% Construct a Delaunay Triangulation of the point set.
dt = DelaunayTri(x,y);
tri = dt(:,:);
% Insert the location of the Voronoi vertices into the existing
% triangulation
V = dt.voronoiDiagram();
% Remove the infinite vertex
V(1,:) = [];
numv = size(V,1);
dt.X(end+(1:numv),:) = V;
Warning: Duplicate data points have been detected and removed.
The Triangulation indices are defined with respect to
the unique set of points in DelaunayTri property X.
% The Delaunay edges that connect pairs of sample points represent the
% boundary.
delEdges = dt.edges();
validx = delEdges(:,1) <= numpts;
validy = delEdges(:,2) <= numpts;
boundaryEdges = delEdges((validx & validy), :)';
xb = x(boundaryEdges);
yb = y(boundaryEdges);
clf;
triplot(tri,x,y);
axis equal;
hold on;
plot(x,y,'*r');
plot(xb,yb, '-r');
xlabel('Curve reconstruction from a point cloud', 'fontweight','b');
hold off;

Example Eight: Compute an Approximate Medial Axis of a Polygonal
Domain
This example demonstrates the creation of an approximate Medial Axis of a
polygonal domain using a constrained Delaunay triangulation. The Medial Axis of
a polygon is defined by the locus of the center of a maximal disk within the
polygon interior.
% Construct a constrained Delaunay triangulation of a sample of points
% on the domain boundary.
load trimesh2d
dt = DelaunayTri(x,y,Constraints);
inside = dt.inOutStatus();
% Construct a TriRep to represent the domain triangles.
tr = TriRep(dt(inside, :), dt.X);

% Construct a set of edges that join the circumcenters of neighboring
% triangles; the additional logic constructs a unique set of such edges.
numt = size(tr,1);
T = (1:numt)';

neigh = tr.neighbors();
cc = tr.circumcenters();
xcc = cc(:,1);
ycc = cc(:,2);
idx1 = T < neigh(:,1);
neigh = [T(idx1) neigh(idx1,1); T(idx2) neigh(idx2,2); T(idx3) neigh(idx3,3)]';
% Plot the domain triangles in green, the domain boundary in blue and the
% medial axis in red.
clf;
triplot(tr, 'g');
hold on;
plot(xcc(neigh), ycc(neigh), '-r', 'LineWidth', 1.5);
axis([-10 310 -10 310]);
axis equal;
plot(x(Constraints'),y(Constraints'), '-b', 'LineWidth', 1.5);
xlabel('Medial Axis of a Polygonal Domain', 'fontweight','b');
hold off;

Example Nine: Morph a 2D Mesh to a Modified Boundary
This example shows how to morph a mesh of a 2D domain to accommodate a
modification to the domain boundary.
Step 1: Load the data. The mesh to be morphed is defined by trife, xfe, yfe, which
is a triangulation in face-vertex format.
load trimesh2d
clf; triplot(trife,xfe,yfe); axis equal;
axis([-10 310 -10 310]);
axis equal;
xlabel('Initial Mesh', 'fontweight','b');

Step 2: Construct a background triangulation - a Constrained Delaunay
triangulation of the set of points representing the mesh boundary. For each vertex
of the mesh, compute a descriptor that defines it's location with respect to the

background triangulation. The descriptor is the enclosing triangle together with the
barycentric coordinates with respect to that triangle.
dt = DelaunayTri(x,y,Constraints);
clf; triplot(dt); axis equal;
axis([-10 310 -10 310]);
axis equal;
xlabel('Background Triangulation', 'fontweight','b');
descriptors.tri = dt.pointLocation(xfe, yfe);
descriptors.baryCoords = dt.cartToBary(descriptors.tri, [xfe yfe]);

Step 3: Edit the background triangulation to incorporate the desired modification
to the domain boundary.
cc1 = [210 90];
circ1 = (143:180)';
x(circ1) = (x(circ1)-cc1(1))*0.6 + cc1(1);
y(circ1) = (y(circ1)-cc1(2))*0.6 + cc1(2);
tr = TriRep(dt(:,:),x,y);
clf; triplot(tr); axis([-10 310 -10 310]); axis equal;
xlabel('Edited Background Triangulation - Hole Size Reduced', 'fontweight','b');

Step 4: Convert the descriptors back to Cartesian coordinates using the deformed
background triangulation as a basis for evaluation.
Xnew = tr.baryToCart(descriptors.tri, descriptors.baryCoords);
tr = TriRep(trife, Xnew);
clf; triplot(tr);
axis([-10 310 -10 310]);
axis equal;
xlabel('Morphed Mesh', 'fontweight','b');

Differential Equations
MATLAB offers several numerical algorithms to solve a wide variety of
differential equations. This demo shows the formulation and solution for three
different types of differential equations using MATLAB.
Initial Value Problem
VANDERPOLDEMO is a function that defines the van der Pol equation.
type vanderpoldemo
function dydt = vanderpoldemo(t,y,Mu)
%VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO.

% $Revision: 1.2 $ $Date: 2002/06/17 13:20:38 $

dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

The equation is written as a system of two first order ODEs. These are evaluated
for different values of the parameter Mu. For faster integration, we choose an
appropriate solver based on the value of the parameter Mu.
For Mu = 1, any of the MATLAB ODE solvers can solve the van der Pol equation
efficiently. The ODE45 solver used below is one such example. The equation is
solved in the domain [0, 20].
tspan = [0, 20];
y0 = [2; 0];
Mu = 1;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode45(ode, tspan, y0);

% Plot of the solution
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')

For larger magnitudes of Mu, the problem becomes stiff. Special numerical
methods are needed for fast integration. ODE15S, ODE23S, ODE23T, and
ODE23TB can solve stiff problems efficiently.
Here is a solution to the van der Pol equation for Mu = 1000 using ODE15S.
tspan = [0, 3000];
y0 = [2; 0];
Mu = 1000;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode15s(ode, tspan, y0);

plot(t,y(:,1))
title('van der Pol Equation, \mu = 1000')
axis([0 3000 -3 3])
xlabel('t')

Boundary Value Problems
BVP4C solves boundary value problems for ordinary differential equations.

The example function TWOODE has a differential equation written as a system of
two first order ODEs.
type twoode
function dydx = twoode(x,y)
%TWOODE Evaluate the differential equations for TWOBVP.
%
% See also TWOBC, TWOBVP.

% Lawrence F. Shampine and Jacek Kierzenka
% $Revision: 1.5 $ $Date: 2002/04/15 03:37:29 $

dydx = [ y(2); -abs(y(1)) ];
TWOBC has the boundary conditions for TWOODE.
type twobc
function res = twobc(ya,yb)
%TWOBC Evaluate the residual in the boundary conditions for TWOBVP.
%
% See also TWOODE, TWOBVP.

% $Revision: 1.5 $ $Date: 2002/04/15 03:37:33 $

res = [ ya(1); yb(1) + 2 ];

Prior to using BVP4C, we have to provide a guess for the solution we want
represented at a mesh. The solver then adapts the mesh as it refines the solution.
BVPINIT assembles the initial guess in the form that the solver BVP4C will need.
For an initial mesh of [0 1 2 3 4] and a constant guess of y(x) = 1, y'(x) = 0, call
BVPINIT like this:
solinit = bvpinit([0 1 2 3 4],[1; 0]);
With this initial guess, we can solve the problem with BVP4C.

The solution sol (below) is then evaluated at points xint using DEVAL and plotted.
sol = bvp4c(@twoode, @twobc, solinit);

xint = linspace(0, 4, 50);
yint = deval(sol, xint);
plot(xint, yint(1,:),'b');
xlabel('x')
hold on

This particular boundary value problem has exactly two solutions. The other
solution is obtained for an initial guess of
y(x) = -1, y'(x) = 0
and plotted as before.
solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);

xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
hold off

Partial Differential Equations
PDEPE solves partial differential equations in one space variable and time.
The examples PDEX1, PDEX2, PDEX3, PDEX4, PDEX5 form a mini-tutorial on
using PDEPE. Browse through these functions for more examples.
This example problem uses functions PDEX1PDE, PDEX1IC, and PDEX1BC.
PDEX1PDE defines the differential equation.
type pdex1pde
function [c,f,s] = pdex1pde(x,t,u,DuDx)
%PDEX1PDE Evaluate the differential equations components for the PDEX1
problem.
%

% See also PDEPE, PDEX1.

% $Revision: 1.5 $ $Date: 2002/04/15 03:37:37 $

c = pi^2;
f = DuDx;
s = 0;
PDEX1IC sets up the initial conditions.
type pdex1ic
function u0 = pdex1ic(x)
%PDEX1IC Evaluate the initial conditions for the problem coded in PDEX1.
%

% $Revision: 1.5 $ $Date: 2002/04/15 03:37:41 $

u0 = sin(pi*x);
PDEX1BC sets up the boundary conditions.
type pdex1bc
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
%PDEX1BC Evaluate the boundary conditions for the problem coded in PDEX1.
%

% $Revision: 1.5 $ $Date: 2002/04/15 03:37:45 $

pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;

PDEPE requires x, the spatial discretization, and t, a vector of times at which you
want a snapshot of the solution. We solve this problem using a mesh of 20 nodes
and request the solution at five values of t. Finally, we extract and plot the first
component of the solution.
x = linspace(0,1,20);
t = [0 0.5 1 1.5 2];
sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
u1 = sol(:,:,1);
surf(x,t,u1);
xlabel('x');ylabel('t');zlabel('u');
hold on

u1 = sol(:,:,1);

surf(x,t,u1);
xlabel('x'); ylabel('t'); zlabel('u');

Graphical Approach to Solving Inequalities
Here is an interesting graphical approach to find out whether e^pi is greater than
piê or not.
The question is: which is greater, e^pi or piê? The easy way to find out is to type
it directly at the MATLAB command prompt. But it motivates a more interesting
question. What is the shape of the function z=x^y-y^x? Here is a plot of z.
%Define the mesh
x=0:0.16:5;
y=0:0.16:5;
[xx,yy]=meshgrid(x,y);

%The plot
zz=xx.^yy-yy.^xx;
h=surf(x,y,zz);

%Set the properties of the plot
set(h,'EdgeColor',[0.7 0.7 0.7]);
view(20,50);
colormap(hsv);
title('z=x^y-y^x'); xlabel('x'); ylabel('y');
hold on;

It turns out that the solution of the equation x^y-y^x=0 has a very interesting
shape. Because interesting things happen near e and pi, our original question is not
easily solved by inspection. Here is a plot of that equation shown in black.
c=contourc(x,y,zz,[0 0]);
list1Len=c(2,1);
xContour=[c(1,2:1+list1Len) NaN c(1,3+list1Len:size(c,2))];
yContour=[c(2,2:1+list1Len) NaN c(2,3+list1Len:size(c,2))];
% Note that the NAN above prevents the end of the first contour line from being
% connected to the beginning of the second line
line(xContour,yContour,'Color','k');

Here is a plot of the integer solutions to the equation x^y-y^x=0. Notice 2^4=4^2 is
the ONLY integer solution where x~=y. So, what is the intersection point of the
two curves that define where x^y=y^x?
plot([0:5 2 4],[0:5 4 2],'r.','MarkerSize',25);

Finally, we can see that e^pi is indeed larger than piê (though not by much) by
plotting these points on our surface.
e=exp(1);
plot([e pi],[pi e],'r.','MarkerSize',25);
plot([e pi],[pi e],'y.','MarkerSize',10);
text(e,3.3,'(e,pi)','Color','k', ...
'HorizontalAlignment','left','VerticalAlignment','bottom');
text(3.3,e,'(pi,e)','Color','k','HorizontalAlignment','left',...
'VerticalAlignment','bottom');
hold off;

Here is a verification of this fact.
e=exp(1);
e^pi
piê
ans =

23.1407

ans =

22.4592

Splines in Two Dimensions

This demonstration interpolates data with a cubic spline in 2 dimensions. It uses
the MATLAB SPLINE function. It does NOT use the Curve Fitting Toolbox
spline functions, which are a complete set of functions for B-splines and other
piecewise polynomials of any degree.
Randomly pick eight points. Plot them.
n = 7;
x = rand(n,1);
y = rand(n,1);
plot(x,y,'.')
axis([0 1 0 1])

Interpolate the points with two splines by evaluating them with a finer spacing.
Plot the interpolated curve with a red line.
t = 1:n;
ts = 1:1/10:n;
xs = spline(t,x,ts);
ys = spline(t,y,ts);
hold on

plot(xs,ys,'r');
hold off

Numerical Integration of Differential Equations
ODE23 and ODE45 are functions for the numerical solution of ordinary
differential equations. They employ variable step size Runge-Kutta integration
methods. ODE23 uses a simple 2nd and 3rd order pair of formulas for medium
accuracy and ODE45 uses a 4th and 5th order pair for higher accuracy. This demo
shows their use on a simple differential equation.
Consider the pair of first order ordinary differential equations known as the Lotka-
Volterra predator-prey model.
y1' = (1 - alpha*y2)*y1
y2' = (-1 + beta*y1)*y2
The functions y1 and y2 measure the sizes of the prey and predator populations
respectively. The quadratic cross term accounts for the interactions between the

species. Note that the prey population increases when there are no predators, but
the predator population decreases when there are no prey.
To simulate a system, create a function that returns a column vector of state
derivatives, given state and time values. For this example, we've created a file
called LOTKA.M.
type lotka
function yp = lotka(t,y)
%LOTKA Lotka-Volterra predator-prey model.

% $Revision: 5.7 $ $Date: 2002/04/15 03:33:21 $

yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y;

To simulate the differential equation defined in LOTKA over the interval 0 < t <
15, invoke ODE23. Use the default relative accuracy of 1e-3 (0.1 percent).
% Define initial conditions.
t0 = 0;
tfinal = 15;
y0 = [20 20]';
% Simulate the differential equation.
tfinal = tfinal*(1+eps);
[t,y] = ode23('lotka',[t0 tfinal],y0);
Plot the result of the simulation two different ways.
subplot(1,2,1)
plot(t,y)
title('Time history')

subplot(1,2,2)
plot(y(:,1),y(:,2))
title('Phase plane plot')

Now simulate LOTKA using ODE45, instead of ODE23. ODE45 takes longer at
each step, but also takes larger steps. Nevertheless, the output of ODE45 is smooth
because by default the solver uses a continuous extension formula to produce
output at 4 equally spaced time points in the span of each step taken. The plot
compares this result against the previous.
[T,Y] = ode45('lotka',[t0 tfinal],y0);

subplot(1,1,1)
title('Phase plane plot')
plot(y(:,1),y(:,2),'-',Y(:,1),Y(:,2),'-');
legend('ode23','ode45')

Loma Prieta Earthquake
This demo shows how to analyze and visualize real-world earthquake data.
Authors: C. Denham, 1990, C. Moler, August, 1992.
The file QUAKE.MAT contains 200Hz data from the October 17, 1989 Loma
Prieta earthquake in the Santa Cruz Mountains. The data are courtesy of Joel
Yellin at the Charles F. Richter Seismological Laboratory, University of
California, Santa Cruz. Start by loading the data.
load quake e n v
whos e n v

e 10001x1 80008 double
n 10001x1 80008 double
v 10001x1 80008 double

In the workspace now are three variables containing time traces from an
accelerometer in the Natural Sciences' building at UC Santa Cruz. The
accelerometer recorded the main shock of the earthquake. The variables n, e, v
refer to the three directional components measured by the instrument, which was
aligned parallel to the fault, with its N direction pointing in the direction of
Sacramento. The data is uncorrected for the response of the instrument.
Scale the data by the gravitational acceleration. Also, create a fourth variable, t,
containing a time base.
g = 0.0980;
e = g*e;
n = g*n;
v = g*v;
delt = 1/200;
t = delt*(1:length(e))';
Here are plots of the accelerations.
yrange = [-250 250];
limits = [0 50 yrange];
subplot(3,1,1), plot(t,e,'b'), axis(limits), title('East-West acceleration')
subplot(3,1,2), plot(t,n,'g'), axis(limits), title('North-South acceleration')
subplot(3,1,3), plot(t,v,'r'), axis(limits), title('Vertical acceleration')

Look at the interval from t=8 seconds to t=15 seconds. Draw black lines at the
selected spots. All subsequent calculations will involve this interval.
t1 = 8*[1;1];
t2 = 15*[1;1];
subplot(3,1,1), hold on, plot([t1 t2],yrange,'k','LineWidth',2); hold off

Zoom in on the selected time interval.
trange = sort([t1(1) t2(1)]);
k = find((trange(1)<=t) & (t<=trange(2)));
e = e(k);
n = n(k);
v = v(k);
t = t(k);
ax = [trange yrange];

subplot(3,1,1), plot(t,e,'b'), axis(ax), title('East-West acceleration')
subplot(3,1,2), plot(t,n,'g'), axis(ax), title('North-South acceleration')
subplot(3,1,3), plot(t,v,'r'), axis(ax), title('Vertical acceleration')

Focusing on one second in the middle of this interval, a plot of "East-West" against
"North-South" shows the horizontal acceleration.
subplot(1,1,1)
k = length(t);
k = round(max(1,k/2-100):min(k,k/2+100));
plot(e(k),n(k),'.-')
xlabel('East'), ylabel('North');
title('Acceleration During a One Second Period');

Integrate the accelerations twice to calculate the velocity and position of a point in
3-D space.
edot = cumsum(e)*delt; edot = edot - mean(edot);
ndot = cumsum(n)*delt; ndot = ndot - mean(ndot);
vdot = cumsum(v)*delt; vdot = vdot - mean(vdot);

epos = cumsum(edot)*delt; epos = epos - mean(epos);
npos = cumsum(ndot)*delt; npos = npos - mean(npos);
vpos = cumsum(vdot)*delt; vpos = vpos - mean(vpos);

subplot(2,1,1);
plot(t,[edot+25 ndot vdot-25]); axis([trange min(vdot-30) max(edot+30)])
xlabel('Time'), ylabel('V - N - E'), title('Velocity')

subplot(2,1,2);
plot(t,[epos+50 npos vpos-50]);
axis([trange min(vpos-55) max(epos+55)])
xlabel('Time'), ylabel('V - N - E'), title('Position')

The trajectory defined by the position data can be displayed with three different 2-
dimensional projections. Here is the first with a few values of t annotated.
subplot(1,1,1); cla; subplot(2,2,1)
plot(npos,vpos,'b');
na = max(abs(npos)); na = 1.05*[-na na];
ea = max(abs(epos)); ea = 1.05*[-ea ea];
va = max(abs(vpos)); va = 1.05*[-va va];
axis([na va]); xlabel('North'); ylabel('Vertical');

nt = ceil((max(t)-min(t))/6);
k = find(fix(t/nt)==(t/nt))';
for j = k, text(npos(j),vpos(j),['o ' int2str(t(j))]); end

Similar code produces two more 2-D views.
subplot(2,2,2)
plot(epos,vpos,'g');
for j = k; text(epos(j),vpos(j),['o ' int2str(t(j))]); end
axis([ea va]); xlabel('East'); ylabel('Vertical');

subplot(2,2,3)
plot(npos,epos,'r');
for j = k; text(npos(j),epos(j),['o ' int2str(t(j))]); end
axis([na ea]); xlabel('North'); ylabel('East');

The fourth subplot is a 3-D view of the trajectory.
subplot(2,2,4)
plot3(npos,epos,vpos,'k')
for j = k;text(npos(j),epos(j),vpos(j),['o ' int2str(t(j))]); end
axis([na ea va]); xlabel('North'); ylabel('East'), zlabel('Vertical');
box on

Finally, plot a dot at every tenth position point. The spacing between dots indicates
the velocity.
subplot(1,1,1)
plot3(npos,epos,vpos,'r')
hold on
step = 10;
plot3(npos(1:step:end),epos(1:step:end),vpos(1:step:end),'.')
hold off
box on
axis tight
xlabel('North-South')
ylabel('East-West')
zlabel('Vertical')
title('Position (cms)')

GRAPHICS
2-D Plots
Here are some examples of 2-D line plots in MATLAB as well as a few useful
tips and tricks for getting started.
You can find a gallery of all MATLAB plotting functions by following the
documentation link at the end of this demo. For more information on a particular
plotting function or property, type doc followed by the name of the function or
property at the command prompt.
Line Plot of a Chirp
This example shows a basic line plot of a chirp signal. It also shows how to enter
labels for the x and y axes using xlabel and ylabel as well as one method of
initializing x-values. Type doc linspace at the command prompt for information on
another method for initializing evenly spaced data sets. Here y-values are
computed as a function of x and stored in the variable y, but you could also plot
just the computed y-values in one command, for example, plot(sin((0:0.05:5).^2)).
x=0:0.05:5;
y=sin(x.^2);
plot(x,y);
xlabel('Time')
ylabel('Amplitude')

Bar Plot of a Bell Shaped Curve
As just suggested, you don't need a variable to hold y-values if you pass a function
that generates them to plot as its y-argument.
x = -2.9:0.2:2.9;
bar(x,exp(-x.*x));

Stairstep Plot of a Sine Wave
x=0:0.25:10;
stairs(x,sin(x));

Errorbar Plot
The errorbar function draws a line plot of x and y values and superimposes on each
observation a vertical error bar, the extent of which depends on an additional
argument. Here the error data is created in line 3 using the rand function, which
generates a sequence of random numbers.
x=-2:0.1:2;
y=erf(x);
e = rand(size(x))/10;
errorbar(x,y,e);

Polar Plot
This example plots the polar coordinates of the angle theta (t, in radians) versus the
radius, rho.
t=0:0.01:2*pi;
polar(t,abs(sin(2*t).*cos(2*t)));

Stem Plot
A stem plot draws a marker for each x, y value which is connected to a common
baseline by a vertical line.
x = 0:0.1:4;
y = sin(x.^2).*exp(-x);
stem(x,y)

Scatter Plot
This example shows the relationship between traffic counts on two streets hour by
hour using red stars for markers. The previous example used the default markers of
empty blue circles. For more information on changing line styles, colors, and
markers, type doc line, doc scatter, or doc line_props in the command line.
load count.dat
scatter(count(:,1),count(:,2),'r*')
xlabel('Number of Cars on Street A');
ylabel('Number of Cars on Street B');

3-D Plots
Here are some examples of surface plots in MATLAB.

Mesh Plot of Peaks
z=peaks(25);
mesh(z);
colormap(hsv)

Surface Plot of Peaks
z=peaks(25);
surf(z);
colormap(jet);

Surface Plot (with Shading) of Peaks
z=peaks(25);
surfl(z);
shading interp;
colormap(pink);

Contour Plot of Peaks
z=peaks(25);
contour(z,16);
colormap(hsv)

Quiver
x = -2:.2:2;
y = -1:.2:1;
[xx,yy] = meshgrid(x,y);
zz = xx.*exp(-xx.^2-yy.^2);
[px,py] = gradient(zz,.2,.2);
quiver(x,y,px,py,2);

Slice
[x,y,z] = meshgrid(-2:.2:2,-2:.25:2,-2:.16:2);
v = x.*exp(-x.^2-y.^2-z.^2);
xslice = [-1.2,.8,2]; yslice = 2; zslice = [-2,0];
slice(x,y,z,v,xslice,yslice,zslice)
colormap hsv

Earth's Topography
MATLAB can be used to create different kinds of maps. Here, we show ways of
representing the topography of the Earth. The topography data used in this demo is
available from the National Geophysical Data Center, NOAA US Department of
Commerce under data announcement 88-MGG-02.
The Data
The data is stored in a MAT file called topo.mat. The variable topo contains the
altitude data for the Earth. topomap1 contains the colormap for the altitude.
load('topo.mat','topo','topomap1');
whos topo topomap1

topo 180x360 518400 double
topomap1 64x3 1536 double

Contour Plot
CONTOUR creates a contour plot of the Earth from this data. Here, we are
producing a contour based on points on the map that have zero altitude. What you
thus see is an outline of the continents.
contour(0:359,-89:90,topo,[0 0],'b')
axis equal
box on
set(gca,'XLim',[0 360],'YLim',[-90 90], ...
'XTick',[0 60 120 180 240 300 360], ...
'Ytick',[-90 -60 -30 0 30 60 90]);

2-D Image Plot
IMAGE creates a 2-D image plot from the data in topo and topomap1. You can see
the elevation information in this figure. Here, altitudes correspond to shades of
green, while depths (below sea level) correspond to shades of blue.
hold on
image([0 360],[-90 90],topo,'CDataMapping', 'scaled');

colormap(topomap1);

3-D Plot
The globe!
The SPHERE function returns x,y,z data that are points on the surface of a sphere
(50 points in this case). Observe the altitude data in topo mapped onto the
coordinates of the sphere contained in x, y and z. Two light sources illuminate the
globe.
[x,y,z] = sphere(50);

cla reset
axis square off
props.AmbientStrength = 0.1;
props.DiffuseStrength = 1;
props.SpecularColorReflectance = .5;
props.SpecularExponent = 20;
props.SpecularStrength = 1;
props.FaceColor= 'texture';

props.EdgeColor = 'none';
props.FaceLighting = 'phong';
props.Cdata = topo;
surface(x,y,z,props);
light('position',[-1 0 1]);
light('position',[-1.5 0.5 -0.5], 'color', [.6 .2 .2]);
view(3)

Images and Matrices
For any matrix X, IMAGE(X) displays a graphical image with brightness or color
chosen from the elements of X used as indices into a colormap. This demo
illustrates this idea of representing a matrix as an image and in general displaying
images stored as matrices.
The Simple Spiral Matrix
SPIRAL stores a simple spiral pattern into a matrix. You can see the spiral pattern
of the matrix elements in the figure. The elements of the matrix spiral away from
the center, growing in magnitude linearly. Small numbers (center values) are
mapped to black and dark gray, while the larger values (around the edge of the

matrix) are mapped to light gray and white. The assignment of small values of the
matrix to black, large values of the matrix to white and intermediate values to
shades of gray determines a color map.
colormap(gray);
X = spiral(8);
image(X);

Colormaps
COLORMAP function is used to change the color mapping. The map had been set
with colormap(gray) in the previous screen. Here we change the colormap to hue-
saturation-value (hsv) color map. The colors begin with red, pass through yellow,
green, cyan, blue, magenta, and return to red.
colormap(hsv);

Another color map is 'hot'. The 'hot' colormap ranges from black through shades of
red and yellow to white.
colormap(hot);

The quantities 'hsv' and 'hot' used with the COLORMAP function are, of course,
matrices. (More precisely, they are the names of functions which return matrices.)
Color map matrices have three columns which specify intensities of the red, green
and blue video components. The number of rows depends upon the particular
image. In this example, the elements of X = spiral(8) range from 1 to 64, so we are
using 64 rows.
M = hot;
size(M)
ans =

64 3

The elements of X are used as indices into the color map and so X must have
positive, integer elements between 1 and the length of the map. To see how an
individual color is determined, pick one element of X, say X(7,1). The
corresponding color map entry is M(37,:). This has full intensity in the red gun, a
little over half intensity in the green gun, and no blue. It produces the shade of
orange in the cell in the (7,1) position near the lower left corner.

ColorMapIndex = X(7,1)
M(ColorMapIndex,:)
ColorMapIndex =

37

ans =

1.0000 0.5417 0

In general, the statements
image(X), colormap(M);
produce a display of colored cells where the RGB intensity of the (i,j)-th cell is the
3-vector M(X(i,j),:). The matrix X can be of any size, but its elements must be
positive integers between 1 and m. The matrix M should then have m rows, 3
columns, and elements between 0.0 and 1.0. COLORMAP(M) also sets the colors
used by PCOLOR(X), SURF(Z) and MESH(Z), but in these cases the data matrix,
X or Z, is rescaled to provide indices into the color map. A completely different
feature of our spiral example is revealed by the 'flag' color map. The 'flag'
colormap is simply m/4 copies of the matrix flag(4), shown below, stacked on top
of each other. The colors red, white, blue and black are used cyclically as the
elements of X vary and so finer details of the image data become apparent. In this
example, we can see the diagonal patterns in the matrix
colormap(flag);
flag(4)
ans =

1 0 0
1 1 1
0 0 1
0 0 0

Since color maps are matrices, it is possible to modify them, or create new ones,
with MATLAB array operations. For example the hot color map can be softened
by adding some gray.
S = (hot + gray)/2;
colormap(S)

A 'hot' colormap, softened by 'gray', can be brightened by raising the elements of
the color map to a power less than 1.
X = spiral(8);
image(X);
gamma = .6;
S = (hot + gray)/2;
S = S.^gamma;
colormap(S)

The command RGBPLOT, produces a plot of the color map. The x-axis is the map
index, which corresponds to the elements of X in IMAGE(X), and the y-axis is the
intensity of the red, green and blue components.
rgbplot(S)

Using SPY
A sparse matrix display function, SPY, is useful for displaying the location of
image elements which point to a particular color map entry. The following code
segment loads a file containing altitude data for eastern New England and displays
all the elements which use the second or third element of the color map. Locations
with X==1 correspond to sea level, so we see a crude representation of the coast
line.
load cape
spy((X==2) | (X==3))

Examples of Large Images
Our 8-by-8 spiral matrix is only a small, illustrative example. Larger matrices
resulting from extensive computations, or images obtained from photographs,
satellites, or scanners are more typical. The demos directory contains several
sample images with their own color maps and the color directory contains M-files
which generate other useful color maps. Below, you see a listing of MAT files that
have these images and their corresponding colormaps. You can choose to use any
of these images (from the structure imglist), and select any of the color maps listed
under colorlabels (including the default colormap 'map' that is loaded with the
image). To do this, you will have to select the following code and execute it. Make
changes to different parameters and see the changes.
clear X map;
imglist = {'flujet', ... Fluid Jet
'spine', ... Bone
'gatlin', ... Gatlinburg
'durer', ... Durer
'detail', ... Durer Detail
'cape', ... Cape Cod

'clown', ... Clown
'earth', ... Earth
'mandrill', ... Mandrill
'spiral'};

colorlabels = {'default', 'hsv','hot','pink',...
'cool','bone','prism','flag',...
'gray','rand'};

load(imglist{4},'X','map');
imagesc(X);
colormap(map);
%colormap(colorlabels{1});
axis off;

Viewing a Penny
The file PENNY.MAT contains measurements made at the National Institute of
Standards and Technology of the depth of the mold used to mint a U. S. penny,
sampled on a 128-by-128 grid.

Drawing a Contour Plot
Draw a contour plot with 15 copper colored contour lines.
load penny.mat
contour(P,15)
colormap(copper)
axis ij square

Drawing a Pseudocolor Plot
Draw a pseudocolor plot with brightness proportional to height.
pcolor(P)
axis ij square
shading flat

Drawing a Pseudocolor Plot With a Colormap
Draw a pseudocolor plot with brightness proportional to the Laplacian of the
height. A cell is bright if its height is greater than the average of its four neighbors
and dark if its height is less than the average of its four neighbors. This is an
unusual "lighting model", but it produces an image that looks like a photograph of
a penny.
D = -del2(P);
pcolor(D)
axis ij square
shading flat

Drawing a Surface Plot With a Colormap
Finally, produce a 3-D, copper colored, surface plot with the Laplacian lighting
model.
surf(P,D);
axis('ij','tight')
shading('flat')
view(-20,75)

Square Wave from Sine Waves
The Fourier series expansion for a square-wave is made up of a sum of odd
harmonics. We show this graphically using MATLAB.
We start by forming a time vector running from 0 to 10 in steps of 0.1, and take the
sine of all the points. Let's plot this fundamental frequency.
t = 0:.1:10;
y = sin(t);
plot(t,y);

Now add the third harmonic to the fundamental, and plot it.
y = sin(t) + sin(3*t)/3;
plot(t,y);

Now use the first, third, fifth, seventh, and ninth harmonics.
y = sin(t) + sin(3*t)/3 + sin(5*t)/5 + sin(7*t)/7 + sin(9*t)/9;
plot(t,y);

For a finale, we will go from the fundamental to the 19th harmonic, creating
vectors of successively more harmonics, and saving all intermediate steps as the
rows of a matrix.
These vectors are plotted on the same figure to show the evolution of the square
wave. Note that Gibbs' effect says that it will never really get there.
t = 0:.02:3.14;
y = zeros(10,length(t));
x = zeros(size(t));
for k=1:2:19
x = x + sin(k*t)/k;
y((k+1)/2,:) = x;
end
plot(y(1:2:9,:)')
title('The building of a square wave: Gibbs'' effect')

Here is a 3-D surface representing the gradual transformation of a sine wave into a
square wave.
surf(y);
shading interp
axis off ij

Functions of Complex Variables
MATLAB can help you perform some very interesting manipulations on
complex variable.
Let f(z) be a function of a complex variable. Consider the domain formed by the
unit disc (displayed below in polar coordinates). The height of the surface is the
real part, REAL(f(z)). The color of the surface is the imaginary part, IMAG(f(z)).
The color map varies the hue in the HSV color model.
CPLXMAP plots a function of a complex variable. It has the syntax
CPLXMAP(z,f(z),bound), where z is the domain, and f(z) is the mapping that
generates the range.
CPLXGRID generates a polar coordinate complex grid. Z = CPLXGRID(m) is an
(m+1)-by-(2*m+1) complex polar grid.
colormap(hsv(64))
z = cplxgrid(30);
cplxmap(z,z)

title('z')

f(z) = z^3. Three maxima at the cube roots of 1.
cplxmap(z,z.^3)
title('z^3')

f(z) = (z^4-1)^(1/4). Four zeros at the fourth roots of 1.
cplxmap(z,(z.^4-1).^(1/4));
title('(z^4-1)^(1/4)')

f(z) = 1/z. A simple pole at the origin.
cplxmap(z,1./(z+eps*(abs(z)==0)),5*pi);
title('1/z')

f(z) = atan(2*z). Branch cut singularities at +-i/2.
cplxmap(z,atan(2*z),1.9)
title('atan(2*z)')

f(z) = z^1/2. Viewed from the negative imaginary axis.
axis('auto')
cplxroot(2)
view(0,0)
title('sqrt(z)')

Another view for f(z) = z^1/2. The Riemann surface for the square root.
view(-37.5,30)
cplxroot(2)
title('sqrt(z)')

f(z) = z^1/3. The Riemann surface for the cube root.
cplxroot(3)
title('z^(1/3)')

3-D Visualization
Klein Bottle
A Klein bottle is a nonorientable surface in four-dimensional space. It is formed by
attaching two Mobius strips along their common boundary.
Klein bottles cannot be constructed without intersection in three-space. The figure
shown is an example of such a self-intersecting Klein bottle.
Thanks to Davide Cervone, University of Minnesota.
Generate The Klein Bottle
Define Klein bottle parameters
n = 12;
a = .2; % the diameter of the small tube
c = .6; % the diameter of the bulb
t1 = pi/4 : pi/n : 5*pi/4; % parameter along the tube
t2 = 5*pi/4 : pi/n : 9*pi/4; % angle around the tube
u = pi/2 : pi/n : 5*pi/2;
[X,Z1] = meshgrid(t1,u);
[Y,Z2] = meshgrid(t2,u);

% The handle
len = sqrt(sin(X).^2 + cos(2*X).^2);
x1 = c*ones(size(X)).*(cos(X).*sin(X) ...
- 0.5*ones(size(X))+a*sin(Z1).*sin(X)./len);
y1 = a*c*cos(Z1).*ones(size(X));
z1 = ones(size(X)).*cos(X) + a*c*sin(Z1).*cos(2*X)./len;
handleHndl=surf(x1,y1,z1,X);
set(handleHndl,'EdgeColor',[.5 .5 .5]);
hold on;

% The bulb

r = sin(Y) .* cos(Y) - (a + 1/2) * ones(size(Y));
x2 = c * sin(Z2) .* r;
y2 = - c * cos(Z2) .* r;
z2 = ones(size(Y)) .* cos(Y);
bulbHndl=surf(x2,y2,z2,Y);
set(bulbHndl,'EdgeColor',[.5 .5 .5])

colormap(hsv);
axis vis3d
view(-37,30);
axis off
light('Position',[2 -4 5])
light
hold off

Half of The Bottle
shading interp
c = X;
[row col] = size(c);
c(1:floor(row/2),:) = NaN*ones(floor(row/2),col);
set(handleHndl,'CData',c);

c = Y;
[row col] = size(c);
c(1:floor(row/2),:) = NaN*ones(floor(row/2),col);
set(bulbHndl,'CData',c);
set([handleHndl bulbHndl],'FaceAlpha',1);

Transparent Bottle
shading faceted;
set(handleHndl,'CData',X);
set(bulbHndl,'CData',Y);
set([handleHndl bulbHndl], ...
'EdgeColor',[.5 .5 .5], ...
'FaceAlpha',.5);

Changing Transparency
Modifying the transparency value for graphics objects reveals structure that is
obscured with opaque objects. For patches and surfaces, use the FaceAlpha and
EdgeAlpha properties to specify the transparency of faces and edges. The
following examples illustrate this.
A Transparent Isosurface
The flow function generates data for the speed profile of a submerged jet within an
infinite tank. One way to visualize this data is by creating an isosurface illustrating
where the rate of flow is equal to a specified value.
cla
[x y z v] = flow;

% Compute and create a patch object from data created from the isosurface with a
% isosurface scalar value of -3
p = patch(isosurface(x,y,z,v,-3));
% Get the normals to the isosurface based on the gradients of v - scalar values
% at the x,y,z coordinate
isonormals(x,y,z,v,p);
set(p,'facecolor','red','edgecolor','none');
daspect([1 1 1]);

Change view to 3-D.
view(3); axis tight; grid on;

Add a light source to the scene to give the surface more definition.
camlight;

Apply smooth rather than faceted shading by changing lighting from default 'flat'
to 'gouraud'.
lighting gouraud;

Adding transparency to the isosurface reveals that there is greater complexity in
the fluid flow than is visible using the opaque surface. The statement alpha(0.5)
sets the FaceAlpha value for the isosurface face to 0.5.
alpha(0.5)

A Texture Map
The MAT file topo.mat has color data representing the different continents on the
globe. This example shows how to map this data onto the surface of a sphere.
The sphere function generates coordinates of a sphere. Setting the cdata property of
the surface object to topo maps the color data contained in topo for each continent
onto the surface of the sphere. In addition, the facecolor property needs to be set to
texturemap if the size of the z-data is different from the size of the cdata
MATLAB then fills in the color at z-data points where the color is not specified
using interpolation. The software can also map the color data onto the surface if
you specify the facealpha transparency property to be of type texture. If this
property is not specified, the globe will not be transparent. Other options for
facealpha are available. Search for facealpha texture in the doc for more
information.
cla reset;
load topo;
[x y z] = sphere(45);
s = surface(x,y,z,'facecolor','texturemap','cdata',topo);

set(s,'edgecolor','none','facealpha','texture','alphadata',topo);
set(s,'backfacelighting','unlit');
colormap(topomap1);
alpha('direct');
alphamap([.1;1])
axis off vis3d;
campos([2 13 10]);
camlight;
lighting gouraud;

Volume Visualization
Run this M-file to see several examples of volume visualization in MATLAB.
See also ISOSURFACE, STREAMTUBE, STREAMRIBBON, STREAMLINE,
CONEPLOT.
Isosurface of MRI Data
figure;

colordef(gcf,'black')
load mri;
D = squeeze(D);
[x y z D] = subvolume(D, [nan nan nan nan nan 4]);
p = patch(isosurface(x,y,z,D, 5), 'FaceColor', 'red', 'EdgeColor', 'none');
p2 = patch(isocaps(x,y,z,D, 5), 'FaceColor', 'interp', 'EdgeColor', 'none');
isonormals(x,y,z,D,p);
view(3);
daspect([1 1 .4])
colormap(gray(100))
camva(9);
box on
camlight(40, 40);
camlight(-20,-10);
lighting gouraud

Coneplot of Wind Data
figure;
cla
load wind

[cx cy cz] = meshgrid(linspace(71,134,10),linspace(18,59,10),3:4:15);
daspect([1 1 1])
h=coneplot(x,y,z,u,v,w,cx,cy,cz,y,3);
set(h,'EdgeColor', 'none');
colormap(hsv);
box on;
axis tight
camproj perspective;
camva(35);
campos([175 10 85]);
camtarget([105 40 0])
camlight left;
lighting gouraud

Streamlines of Wind Data
figure;
cla
[sx sy sz] = meshgrid(80, 20:10:50, 0:5:15);
h=streamline(x,y,z,u,v,w,sx,sy,sz);
set(h, 'Color', 'cyan');

daspect([1 1 1])
box on;
camva(32);
axis tight
campos([175 10 85]);
camlight left;
lighting gouraud

Streamtubes of Wind Data
Useful for visualizing divergence of a vector field
figure;
cla
[sx sy sz] = meshgrid(80, [20 30 40], [5 10]);
daspect([1,1,1]);
h=streamtube(x,y,z,u,v,w,sx,sy,sz);
set(h,'facecolor','cyan','edgecolor','none');

box on;
axis([70 138 17 60 2.5 16]);
axis tight
camva(28);
campos([175 10 95]);
camlight left;
lighting gouraud

Streamribbons of Wind Data
Useful for visualizing curl of a vector field
figure;
cla
[sx sy sz] = meshgrid(80, [20 30 40], [5 10]);
daspect([1,1,1]);
h=streamribbon(x,y,z,u,v,w,sx,sy,sz);
set(h,'facecolor','cyan','edgecolor','none')

box on;
axis([70 138 17 60 2.5 16]);
axis tight
camva(28);
campos([175 10 85]);
camlight left;
lighting gouraud

Isosurface, Isocaps, Coneplot, and Streamlines of Wind Data
figure;
cla
spd = sqrt(u.*u + v.*v + w.*w);
p = patch(isosurface(x,y,z,spd, 40));
isonormals(x,y,z,spd, p)
set(p, 'FaceColor', 'red', 'EdgeColor', 'none');
p2 = patch(isocaps(x,y,z,spd, 40));
set(p2, 'FaceColor', 'interp', 'EdgeColor', 'none')
daspect([1 1 1]);

[f verts] = reducepatch(isosurface(x,y,z,spd, 30), .2);
h=coneplot(x,y,z,u,v,w,verts(:,1),verts(:,2),verts(:,3),2);
set(h, 'FaceColor', 'cyan', 'EdgeColor', 'none');
[sx sy sz] = meshgrid(80, 20:10:50, 0:5:15);
h2=streamline(x,y,z,u,v,w,sx,sy,sz);
set(h2, 'Color', [.4 1 .4]);
colormap(jet)
box on
axis tight
camva(34);
campos([165 -20 65]);
camtarget([100 40 -5])
camlight left;
lighting gouraud

PROGRAMMING
Manipulating Multidimensional Arrays
MATLAB supports arrays with more than two dimensions. Multidimensional
arrays can be numeric, character, cell, or structure arrays.
Multidimensional arrays can be used to represent multivariate data. MATLAB
provides a number of functions that directly support multidimensional arrays.
Creating Multi-Dimensional Arrays
Multidimensional arrays in MATLAB are created the same way as two-
dimensional arrays. For example, first define the 3 by 3 matrix, and then add a
third dimension.
A = [5 7 8;
0 1 9;
4 3 6];
A(:,:,2) = [1 0 4;
3 5 6;
9 8 7]
A(:,:,1) =

5 7 8
0 1 9
4 3 6

A(:,:,2) =

1 0 4
3 5 6
9 8 7

The CAT function is a useful tool for building multidimensional arrays. B =
cat(DIM,A1,A2,...) builds a multidimensional array by concatenating A1, A2 ...
along the dimension DIM.
B = cat( 3, [2 8; 0 5], [1 3; 7 9], [2 3; 4 6])
B(:,:,1) =

2 8
0 5

B(:,:,2) =

1 3
7 9

B(:,:,3) =

2 3
4 6

Calls to CAT can be nested.
A = cat(3,[9 2; 6 5], [7 1; 8 4]);
B = cat(3,[3 5; 0 1], [5 6; 2 1]);
C = cat(4,A,B,cat(3,[1 2; 3 4], [4 3; 2 1]));
Finding the Dimensions
SIZE and NDIMS return the size and number of dimensions of matrices.
SzA = size(A)
DimsA = ndims(A)
SzC = size(C)
DimsC = ndims(C)
SzA =

2 2 2

DimsA =

3

SzC =

2 2 2 3

DimsC =

4

Accessing Elements
To access a single element of a multidimensional array, use integer subscripts. For
example D(1,2,2,22), using D defined in the previous slide, returns 6.
Array subscripts can also be vectors. For example:
K = C(:,:,1,[1 3])
K(:,:,1,1) =

9 2
6 5

K(:,:,1,2) =

1 2
3 4

Manipulating Multi-Dimensional Arrays
RESHAPE, PERMUTE, and SQUEEZE are used to manipulate n-dimensional
arrays. RESHAPE behaves as it does for 2D arrays. The operation of PERMUTE is
illustrated below.

Let A be a 3 by 3 by 2 array. PERMUTE(A,[2 1 3]) returns an array with the row
and column subscripts reversed (dimension 1 is the row, dimension 2 is the
column, dimension 3 is the depth and so on). Similarly, PERMUTE(A,[3,2,1])
returns an array with the first and third subscripts interchanged.
A = rand(3,3,2);
B = permute(A, [2 1 3]);
C = permute(A, [3 2 1]);
Selecting 2D Matrices From Multi-Dimensional Arrays
Functions like EIG that operate on planes or 2D matrices do not accept multi-
dimensional arrays as arguments. To apply such functions to different planes of the
multidimensional arrays, use indexing or FOR loops. For example:
A = cat( 3, [1 2 3; 9 8 7; 4 6 5], [0 3 2; 8 8 4; 5 3 5], ...
[6 4 7; 6 8 5; 5 4 3]);
% The EIG function is applied to each of the horizontal 'slices' of A.
for i = 1:3
eig(squeeze(A(i,:,:)))
end
ans =

10.3589
-1.0000
1.6411

ans =

21.2293
0.3854 + 1.5778i
0.3854 - 1.5778i

ans =

13.3706
-1.6853 + 0.4757i
-1.6853 - 0.4757i

INTERP3, INTERPN, and NDGRID are examples of interpolation and data
gridding functions that operate specifically on multidimensional data. Here is an
example of NDGRID applied to an N-dimensional matrix.
x1 = -2*pi:pi/10:0;
x2 = 2*pi:pi/10:4*pi;
x3 = 0:pi/10:2*pi;
[x1,x2,x3] = ndgrid(x1,x2,x3);
z = x1 + exp(cos(2*x2.^2)) + sin(x3.^3);
slice(z,[5 10 15], 10, [5 12]); axis tight;

You can build multidimensional cell arrays and multidimensional structure arrays,
and can also convert between multidimensional numeric and cell arrays.
Structures
MATLAB supports specialized data constructs such as structures and cell arrays.
MATLAB structures are array-oriented data constructs. They provide a convenient
way to group related data of different types.

Structures are MATLAB data constructs with named "data containers" called
fields. The fields of a structure can contain any kind of data. For example, one field
might contain a text string representing a name, another might contain a scalar
representing a billing amount, a third might hold a matrix of medical test results,
and so on.
% Draw a visualization of a structure.
strucdem_helper(1);

You can construct a structure simply by assigning values to its fields. With these
commands, we create the structure we've depicted.
patient.name = 'John Doe';
patient.billing = 127.00;
patient.test = [79 75 73; 180 178 177.5; 172 170 169];
patient
patient =

name: 'John Doe'
billing: 127
test: [3x3 double]

You can also build an array of structures to collect similar items together. A
structure array has the following properties:
* All structures in the array have the same number of fields.
* All fields have the same field names.
You can build a structure array by adding subscripts after the structure name.
patient(2).name = 'Ann Lane';
patient(2).billing = 28.50;
patient(2).test = [68 70 68; 118 118 119; 172 170 169];
% Update the visualization.
strucdem_helper(2);

You can access any field in a structure as easily as you access a regular variable.
For example, we can draw a bar graph of the test data for patient(1).
bar(patient(1).test)

The FIELDNAMES function returns the field names for a structure array.
You can remove a given field from every structure within a structure array using
the RMFIELD function.
fnames1 = fieldnames(patient)
patient2 = rmfield(patient,'test');
fnames2 = fieldnames(patient2)
fnames1 =

'name'
'billing'
'test'

fnames2 =

'name'
'billing'

Structures can be nested. You can use the STRUCT function or direct assignment
statements to nest structures within existing structure fields.
A = struct( 'data', {[3 4 7; 8 0 1], [9 3 2; 7 6 5]}, ...
'nest', {...
struct( 'testnum', 'Test 1', ...
'xdata', [4 2 8], 'ydata', [7 1 6]), ...
struct( 'testnum', 'Test 2', ...
'xdata', [3 4 2], 'ydata', [5 0 9])});

% Update the visualization.
strucdem_helper(3)

Here are some more structure commands in action. For further information on
structures, please consult the HELPDESK or the MATLAB Manual.
anotherfield = 'myfield';
st = struct('yourfield','foo',anotherfield,'foo');
st.(anotherfield)='bar';
st = rmfield(st,anotherfield);
if isfield(st,anotherfield);

disp(st)
end

Function Functions
In MATLAB, one function take another as a parameter. This feature serves a
wide variety of purposes. Here we illustrate its use for finding zeros, optimization,
and integration.

The HUMPS Function
A MATLAB function is a file that starts with the keyword function. This is what
the function HUMPS looks like:
type humps
function [out1,out2] = humps(x)
%HUMPS A function used by QUADDEMO, ZERODEMO and FPLOTDEMO.
% Y = HUMPS(X) is a function with strong maxima near x = .3
% and x = .9.
%
% [X,Y] = HUMPS(X) also returns X. With no input arguments,
% HUMPS uses X = 0:.05:1.
%
% Example:
% plot(humps)
%
% See QUADDEMO, ZERODEMO and FPLOTDEMO.

% $Revision: 5.8 $ $Date: 2002/04/15 03:34:07 $

if nargin==0, x = 0:.05:1; end

y = 1 ./ ((x-.3).^2 + .01) + 1 ./ ((x-.9).^2 + .04) - 6;

if nargout==2,
out1 = x; out2 = y;
else
out1 = y;
end

Plot of HUMPS
This figure shows a plot of HUMPS in the domain [0,2] using FPLOT.
fplot(@humps,[0,2]);

Zero of HUMPS
The FZERO function finds a zeros of a function near an initial estimate. Our guess
here for HUMPS is 1.
z = fzero(@humps,1,optimset('Display','off'));

hold on;
plot(z,0,'r*');
hold off

Minimum of HUMPS
FMINBND finds the minimum of a function in a given domain. Here, we search
for a minimum for HUMPS in the domain (0.25,1).
m = fminbnd(@humps,0.25,1,optimset('Display','off'));
fplot(@humps,[0 2]);
hold on;
plot(m,humps(m),'r*');
hold off

Integral of HUMPS
QUAD finds the definite integral of HUMPS in a given domain. Here it computes
the area in the domain [0.5, 1].
q = quad(@humps,0.5,1);
title(['Area = ',num2str(q)]);

Nested Functions
This gives examples of how nested functions can be used for easy data sharing, as
well as providing a new way to create customized functions.
Example 1: Sharing Data
Let's first take a look at taxDemo.m, which contains a nested function.
type taxDemo.m
function y = taxDemo(income)
%TAXDEMO Used by NESTEDDEMO.
% Calculate the tax on income.

% $Revision: 1.1.6.2 $ $Date: 2004/03/02 21:47:05 $

AdjustedIncome = income - 6000; % Calculate adjusted income

% Call 'computeTax' without passing 'AdjustedIncome' as a parameter.

y = computeTax;

function y = computeTax

% This function can see the variable 'AdjustedIncome'
% in the calling function's workspace
y = 0.28 * AdjustedIncome;
end
end

The nested function computeTax can see the variables in the parent function's
workspace. This makes sharing of data between multiple nested functions easy and
particularly useful when processing large data sets. We can call the function in the
usual way.
y=taxDemo(80e3) % What is the tax on $80k income?
y =

2.0720e+004

For nested functions, the end statement is required at the end of a function. You
can also nest functions to any level.
Example 2: Creating Customized Functions
Nested functions allow the ability to create customized functions. Let's look at
makefcn.m which contains a nested function.
type makefcn.m
function fcn = makefcn(a,b,c)
%MAKEFCN Used by NESTEDDEMO.
% This function returns a handle to a customized version of 'parabola'.
% a,b,c specifies the coefficients of the function.

% $Revision: 1.1.6.2 $ $Date: 2004/03/02 21:46:56 $

fcn = @parabola; % Return handle to nested function

function y = parabola(x)
% This nested function can see the variables 'a','b', and 'c'
y = a*x.^2 + b.*x + c;
end
end

When you call makefcn, it returns a function handle to a customized function. For
example:
f = makefcn(3,2,10);
g = makefcn(0,5,25);
f and g are handles to two functions, each with different coefficients. We can
evaluate the functions by using their function handles and passing in parameters.
y=f(2)
y =

26

y=g(2)
y =

35

We can also pass the handle to function functions, such as optimization or
integration.
minimum=fminbnd(f,-5,5);
Or plot the function over a range.
ezplot(f); % Plot f over a range of x
hold on;
plot(2,f(2),'d'); % Plot a marker at (2,f(2))
plot(minimum,f(minimum),'s'); % Plot at minimum of f
text(minimum,f(minimum)-2,'Minimum');
h=ezplot(g); set(h,'color','red') % Plot g over a range of x

plot(2,g(2),'rd'); % Plot a marker at (2,g(2))
hold off;

Example 3: Creating Customized Functions with State
Let's look at makecounter.m which contains a nested function.
type makecounter.m
function countfcn = makecounter(initvalue)
%MAKECOUNTER Used by NESTEDDEMO.
% This function returns a handle to a customized nested function 'getCounter'.
% initvalue specifies the initial value of the counter whose's handle is returned.

% $Revision: 1.1.6.2 $ $Date: 2004/03/02 21:46:55 $

currentCount = initvalue; % Initial value
countfcn = @getCounter; % Return handle to getCounter

function count = getCounter
% This function increments the variable 'currentCount', when it

% gets called (using its function handle) .
currentCount = currentCount + 1;
count = currentCount;
end
end

When you call makecounter, it returns a handle to its nested function getCounter.
getCounter is customized by the value of initvalue, a variable it can see via nesting
within the workspace of makecounter.
counter1 = makecounter(0); % Define counter initialized to 0
counter2 = makecounter(10); % Define counter initialized to 10
Here we have created two customized counters: one which starts at 0 and one
which starts at 10. Each handle is a separate instance of the function and its calling
workspace. Now we can call the inner nested function via its handle. counter1 does
not take parameters but it could.
counter1Value=counter1()
counter1Value =

1

We can call the two functions independently as there are two separate workspaces
for the parent functions kept. They remain in memory while the handles to their
nested functions exist. In this case the currentCount variable gets updated when
counter1 is called.
counter1Value =

2

counter2Value =

11

Anonymous Functions
This demo shows some examples of how to define functions at the command line
in a new much simpler way than with the inline function.
Integrating a Function
Consider the function 10*x.

If we want to allow any multiplier of x, not just 10, we might create a variable g
(where g is initially set to 10), and create a new function

Let's do this in MATLAB by creating a function handle h.
g=10;
h=@(x) g*x;
You can integrate the function by passing its handle to quad.
quad(h,1,10)
ans =

495

Consider another function:

Create a function handle to this function where alpha = 0.9.
alpha=0.9;
f=@(x) sin(alpha*x);
Plot the function and shade the area under it.
x=0:pi/100:pi;
area(x,f(x)); % You can evaluate f without feval

title('f(x)= sin(\alpha x), \alpha =.9');

We can use quad to calculate the area under the function between a range of
values.
quad(f,0,pi)
ans =

2.1678

Minimizing a Function
Consider the function:

where a=1, b=-2, and c=1
Create a function handle for it.

a=1;b=-2;c=1;
f=@(x)(a*x.^2+b*x+c);
ezplot(f); % Plot the function
title('f(x)=ax^2+bx+c, a=1,b=-2,c=1');
hold on;

% Find and plot the minimum
minimum=fminbnd(f,-2,2); % We can pass the function handle direct to the
minimization routine
plot(minimum,f(minimum),'d'); % We can evaluate the function without using
feval
grid;
hold off;

2D Functions
We can create handles to functions of many variables
a=pi;b=15;
f=@(x,y) (a*x+b*y);
ezsurf(f);

title('f(x,y)=ax+by, a=\pi, b=15');

Function Composition
We can also create handles to functions of functions
f=@(x) x.^2;
g=@(x) 3*x;
h=@(x) g(f(x));
h(3)
ans =

27

Reading Arbitrary Format Text Files with TEXTSCAN
This example shows how to read an arbitrary format text file with textscan. This
function is similar to textread, however it also allows you to read the file one block
at a time, and each block can have a different format. The information in the text
file test80211.txt is the result from a wireless network communication quality test.

Each block is a different environment (e.g., mobile, indoor, outdoor). The
numerical results show the data error rate over a range of noise levels for a number
of independent tests.

File Format
After 4 lines of introduction, this particular file is made up of a number blocks of
data, each with the following format:
Two headerlines of description
A parameter m
A p x m table of data
All the information is read into cell arrays, allowing the storage of different size
blocks.
Open the Text File for Reading
fid = fopen('test80211.txt','r'); % Open text file
Read Introduction Lines
InputText=textscan(fid,'%s',4,'delimiter','\n'); % Read strings delimited by a
carriage return
Intro=InputText{1};
disp(Intro);
'*CCX'
'*CCX WiFi conformance test'
'*CCX BER Results'
'*CCX'

Read Each Block
For each block, we read a header, a table name, column headers for the data, then
the data itself.
Block = 1; % Initialize block index
while (~feof(fid)) % For each block...

sprintf('Block: %s', num2str(Block)); % Display block number
InputText=textscan(fid,'%s',2,'delimiter','\n'); % Read header line
HeaderLines{Block,1}=InputText{1};
disp(HeaderLines{Block});

InputText=textscan(fid,'Num SNR=%f'); % Read parameter value
NumCols=InputText{1};

FormatString=repmat('%f',1,NumCols); % Create format string based on
parameter
InputText=textscan(fid,FormatString,'delimiter',','); % Read data block

Data{Block,1}=cell2mat(InputText); % Convert to numerical array from cell
[NumRows,NumCols]=size(Data{Block}); % Size of table
disp(cellstr([xlate('Table data size: ') num2str(NumRows) ' x '
num2str(NumCols)]));
disp(' '); % New line

eob=textscan(fid,'%s',1,'delimiter','\n'); % Read and discard EOB marker ('EOF'
in this case)
Block = Block+1; % Increment block index
end
'* Mobile1'
'* SNR Vs test No'

'Table data size: 30 x 19'

'* Mobile2'
'* SNR Vs test No'


'* Mobile3'
'* SNR Vs test No'


'* Mobile4'
'* SNR Vs test No'


'* Mobile5'
'* SNR Vs test No'


'* Mobile6'
'* SNR Vs test No'


'* Mobile7'
'* SNR Vs test No'


'* Mobile8'
'* SNR Vs test No'


'* Indoor0'
'* SNR Vs test No'


'* Indoor1'
'* SNR Vs test No'


'* Indoor2'
'* SNR Vs test No'


'* Indoor3'
'* SNR Vs test No'


'* Outdoor1'
'* SNR Vs test No'


'* Outdoor2'
'* SNR Vs test No'


'* Outdoor3'
'* SNR Vs test No'


'* Outdoor4'
'* SNR Vs test No'


'* Outdoor5'
'* SNR Vs test No'


Close the Text File
fclose(fid);
How Many Blocks
How many blocks were there?
NrOfBlocks = Block-1
NrOfBlocks =

17

Look at Data
Let's take a look at the numerical data in one of the blocks.
Block=9;

% Headers and Data
disp(HeaderLines{Block});
disp(['SNR' sprintf(' %d',Data{Block,1}(1,2:end))])

user_format = get(0, 'format');
format short e % Use exponential format

disp(' ');
disp(Data{Block,1}(2:end,2:end));
set(0, 'format', user_format);
'* Indoor0'
'* SNR Vs test No'

SNR -7 -6

9.0600e-007 6.7100e-007
3.1700e-007 3.5400e-007

2.8600e-007 1.9600e-007
1.4800e-007 7.3400e-007
3.9500e-008 9.6600e-007
7.9600e-007 7.8300e-007
4.0000e-007 8.8100e-007
3.0100e-007 2.9700e-007

Creating Graphical User Interfaces
Displaying Matrix Data in a GUI
This is a demonstration of how to display matrix data in a GUI using the uitable
component. It shows how to display simple numeric matrix data, as well as more
complex data. It also shows how to modify the appearance of the uitable and how
to limit the changes that end users can make to the data in the uitable.
Using Simple Numeric Data
First, create a figure and position a uitable in it.
f = figure('Position', [100 100 752 350]);
t = uitable('Parent', f, 'Position', [25 25 700 200]);

Now, create some simple numeric data and set it into the uitable using the Data
property. For this example, use a magic square.
set(t, 'Data', magic(10));

Using Mixed Data Types
Next, replace the simple numeric data with a more complicated matrix. First, create
some new data that contains a mix of numeric, character, and logical values. Then,
set the new matrix into the uitable's 'Data' property.
complexData = { ...
'Andrew' 31 'Male' 76 230.75 false 'Somewhat Frequently' 10; ...
'Bob' 41 'Male' 77 235.5 true 'Not Very Frequently' 10; ...
'Clarice' 20 'Female' 74 130.25 false 'Very Frequently' 15; ...
'Debra' 35 'Female' 67 140.25 true 'Very Frequently' 15;
'Edward' 34 'Male' 73 247.5 true 'Very Frequently' 20};
set(t, 'Data', complexData);

Customizing the Display
Now that the complex data is in the uitable, customize the display to make the
information more meaningful. Start by using the ColumnName property to add
headings, or titles, to the top of each column. To create multi-line headings, use the
divider line symbol.
set(t, 'ColumnName', {'Name', 'Age', 'Gender', 'Height', 'Weight', ...
'Stress|Tests', 'Exercise|Frequency', 'Copay'});

Next, adjust the widths of some of the columns using the ColumnWidth property.
Set a specific width for two of the columns, leaving the rest of the columns to
autofit based on the contents.
set(t, 'ColumnWidth', {100 'auto' 'auto' 'auto' 'auto' 'auto' 150 'auto'});

Now, completely remove the row header by setting the RowName property to
empty, using [].
set(t, 'RowName', []);

Resize the uitable to remove any extra space.
set(t, 'Position', [25 25 702 119]);

Finally, change the foreground and background colors using the ForegroundColor
and BackgroundColor properties. As you can see, the RowStriping property is on
by default. This means the uitable shows multiple, in this case two, background
colors. To keep the striping effect, specify two different background colors. Note:
Turn off row striping by setting the RowStriping property to 'off'.
foregroundColor = [1 1 1];
set(t, 'ForegroundColor', foregroundColor);
backgroundColor = [.4 .1 .1; .1 .1 .4];
set(t, 'BackgroundColor', backgroundColor);

Allowing Restricted Editing
Now that the uitable looks the way it should, allow users to edit the data from the
GUI using the ColumnEditable property. However, limit the users' editing
capabilities. First, enable editing on all of the columns except the first one so users
will not be able to modify the values in the Name column.
set(t, 'ColumnEditable', [false true true true true true true true]);

Next, use the ColumnFormat property to change the format of the Gender and
Exercise Frequency columns to pop-up menus with a restricted set of options. In
the Gender column, force users to choose between the values Male and Female. In
the Exercise Frequency column, make users choose among the options: Very
Frequently, Somewhat Frequently, and Not Very Frequently.
set(t, 'ColumnFormat', {[] [] {'Male' 'Female'} [] [] [] ...
{'Very Frequently' 'Somewhat Frequently' 'Not Very Frequently'} []});

Finally, restrict the values in the Age column to be between 0 and 120. Do this by
attaching a CellEditCallback to validate edits made to the uitable data. Use the
following callback function:
function AgeVerificationCallback(o, e)
if (e.Indices(2) == 2 && ...
(e.NewData < 0 || e.NewData > 120))
tableData = get(o, 'data');
tableData{e.Indices(1), e.Indices(2)} = e.PreviousData;
set(o, 'data', tableData);
error('Age value must be between 0 and 120.')
end
By attaching this function as the CellEditCallback, if the edit to the value in the
Age column does not fall within the acceptable range (i.e., between 0 and 120), the
change will not be accepted. In other words, the value in the cell will revert back to
the original value.
set(t, 'CellEditCallback', @AgeVerificationCallback);

External Interfaces
Programming with COM on Windows
Component Object Model (COM), is a set of object-oriented technologies and
tools that enable software developers to integrate application-specific components
from different vendors into their own application solution.
COM helps in integrating significantly different features into one application in a
relatively easy manner. For example, using COM, a developer may choose a
database access component by one vendor, a business graph component by
another, and integrate these into a mathematical analysis package produced by yet
a third.
COM provides a framework for integrating reusable, binary software components
into an application. Because components are implemented with compiled code, the
source code may be written in any of the many programming languages that
support COM. Upgrades to applications are simplified, as components can be
simply swapped without the need to recompile the entire application. In addition, a
component's location is transparent to the application, so components may be
relocated to a separate process or even a remote system without having to modify
the application.
Automation is a method of communication between COM clients and servers. It
uses a single, standard COM interface called IDispatch. This interface enables the
client to find out about, and invoke or access, methods and properties supported by
a COM object. A client and server that communicate using IDispatch are known as
an Automation client and Automation server. IDispatch is the only interface
supported by MATLAB. Custom and dual interfaces are not supported.
MATLAB can communicate with both Automation servers and controls.
Demo Requirements
This demo runs on Windows systems only.
The MWSamp2 object is already registered during MATLAB installation.
However to get a better overview of how to work with COM components in

general, it is assumed in this demo that the user has to register the control. It is also
assumed that regsvr32.exe is located on the DOS path. The following are the steps
needed to register a component on your machine.
1. Run the command "regsvr32 < path >" where < path > indicates the full path to
the ocx/dll file supplied with the component.
2. Restart MATLAB.
cmd = sprintf('regsvr32 /s "%s"', ...
fullfile(matlabroot,'toolbox','matlab','winfun',computer('arch'),'mwsamp2.ocx'));

[s,c] = dos(cmd);
%
% This demo also requires Microsoft(R) Excel(R).

if ~ispc
errordlg('COM Demonstration is for PC only.')
return
end
Creating COM Objects in MATLAB
The following commands create an Automation control object and an Automation
server object in MATLAB:
% Create an Automation control object and put it in a figure.
hf = figure;
title('ActiveX Sample Control')
set(gca,'Xtick',[],'Ytick',[],'Box','on')
fp = get(hf,'Position');
mwsampPosition = get(hf,'DefaultAxesPosition').*fp([3 4 3 4]) ;
mwsamp = actxcontrol('MWSAMP.MwsampCtrl.2', mwsampPosition+1, hf)

% Create an Automation server object.
hExcel = actxserver('excel.application')

mwsamp =

COM.MWSAMP_MwsampCtrl_2

hExcel =

COM.excel_application

Displaying Properties of COM Objects
The properties of COM objects can be displayed to the MATLAB command
window using the GET function, and are displayed graphically using the property
inspector. For a demonstration of the property inspector, take a look at the
Graphical Interface section of this demo.
get(mwsamp)
Label: 'Label'
Radius: 20
Ret_IDispatch: [1x1
Interface.mwsamp2_ActiveX_Control_module._DMwsamp2]

Changing COM Object Properties
Properties of a COM object can be changed using the SET function.
% This makes the Excel(R) Automation server application visible.
set(hExcel,'Visible',1)
The SET function returns a structure array if only the handle to the COM Object is
passed as an argument.
out = set(mwsamp)
out =

Label: {}
Radius: {}
Ret_IDispatch: {}

You can also use the SET function to simultaneously change multiple properties of
COM objects.
set(mwsamp,'Label','Mathworks Sample Control','Radius',40)
Displaying and Changing Enumerated Property Types
You can display and change properties with enumerated values using the SET and
GET functions.
get(hExcel,'DefaultSaveFormat')
ans =

xlExcel8

The SET function can be used to display all possible enumerated values for a
specific property.
set(hExcel,'DefaultSaveFormat')
ans =

'xlAddIn'

'xlCSV'
'xlCSVMac'
'xlCSVMSDOS'
'xlCSVWindows'
'xlDBF2'
'xlDBF3'
'xlDBF4'
'xlDIF'
'xlExcel2'
'xlExcel2FarEast'
'xlExcel3'
'xlExcel4'
'xlExcel5'
'xlExcel5'
'xlExcel9795'
'xlExcel4Workbook'
'xlIntlAddIn'
'xlIntlMacro'
'xlWorkbookNormal'
'xlSYLK'
'xlTemplate'
'xlCurrentPlatformText'
'xlTextMac'
'xlTextMSDOS'
'xlTextPrinter'
'xlTextWindows'
'xlWJ2WD1'
'xlWK1'
'xlWK1ALL'
'xlWK1FMT'
'xlWK3'
'xlWK4'
'xlWK3FM3'
'xlWKS'
'xlWorks2FarEast'
'xlWQ1'
'xlWJ3'
'xlWJ3FJ3'
'xlUnicodeText'
'xlHtml'

'xlWebArchive'
'xlXMLSpreadsheet'
'xlExcel12'
'xlOpenXMLWorkbook'
'xlOpenXMLWorkbookMacroEnabled'
'xlOpenXMLTemplateMacroEnabled'
'xlTemplate'
'xlOpenXMLTemplate'
'xlAddIn'
'xlOpenXMLAddIn'
'xlExcel8'
'xlOpenDocumentSpreadsheet'
'xlOpenXMLWorkbook'

The SET function also enables you to set enumerated values for properties that
support enumerated types.
set(hExcel,'DefaultSaveFormat','xlWorkbookNormal');
Creating Custom Properties for a COM Object
You can create custom properties for a COM object in MATLAB. For instance,
you can make the handle to the Excel COM object a property of the MWSamp2
control and also make the handle to the MWSamp2 control a property of the Excel
COM Object.
addproperty(mwsamp,'ExcelHandle');
addproperty(hExcel,'mwsampHandle');
addproperty(mwsamp,'TestValue');
set(mwsamp,'ExcelHandle',hExcel);
set(mwsamp,'TestValue',rand);
set(hExcel,'mwsampHandle',mwsamp);
get(hExcel,'mwsampHandle')

ans =

COM.MWSAMP_MwsampCtrl_2

get(mwsamp,'ExcelHandle')

ans =

COM.excel_application

get(mwsamp,'TestValue')
ans =

0.5476

Custom properties that are created using the ADDPROPERTY function can also be
removed.
deleteproperty(mwsamp,'TestValue');
Displaying Methods of COM Objects
You can display methods of COM objects in MATLAB by using the INVOKE,
METHODS and METHODSVIEW functions. METHODSVIEW provides a way
to view the methods to the COM objects graphically. For a demonstration of the
METHODSVIEW function, take a look at the Graphical Interface section of this
demo.
invoke(hExcel)
Calculate = void Calculate(handle)
DDEExecute = void DDEExecute(handle, int32, string)
DDEInitiate = int32 DDEInitiate(handle, string, string)
DDEPoke = void DDEPoke(handle, int32, Variant, Variant)
DDERequest = Variant DDERequest(handle, int32, string)
DDETerminate = void DDETerminate(handle, int32)
Evaluate = Variant Evaluate(handle, Variant)
ExecuteExcel4Macro = Variant ExecuteExcel4Macro(handle, string)
Intersect = handle Intersect(handle, handle, handle, Variant(Optional))
Range = handle Range(handle, Variant, Variant(Optional))
Run = Variant Run(handle, Variant(Optional))
SendKeys = void SendKeys(handle, Variant, Variant(Optional))
Union = handle Union(handle, handle, handle, Variant(Optional))
ActivateMicrosoftApp = void ActivateMicrosoftApp(handle,
XlMSApplication)
AddCustomList = void AddCustomList(handle, Variant,
Variant(Optional))

CentimetersToPoints = double CentimetersToPoints(handle, double)
CheckSpelling = bool CheckSpelling(handle, string, Variant(Optional))
ConvertFormula = Variant ConvertFormula(handle, Variant,
XlReferenceStyle, Variant(Optional))
DeleteCustomList = void DeleteCustomList(handle, int32)
DoubleClick = void DoubleClick(handle)
GetCustomListContents = Variant GetCustomListContents(handle, int32)
GetCustomListNum = int32 GetCustomListNum(handle, Variant)
GetOpenFilename = Variant GetOpenFilename(handle, Variant(Optional))
GetSaveAsFilename = Variant GetSaveAsFilename(handle,
Variant(Optional))
Goto = void Goto(handle, Variant(Optional))
Help = void Help(handle, Variant(Optional))
InchesToPoints = double InchesToPoints(handle, double)
InputBox = Variant InputBox(handle, string, Variant(Optional))
MacroOptions = void MacroOptions(handle, Variant(Optional))
MailLogoff = void MailLogoff(handle)
MailLogon = void MailLogon(handle, Variant(Optional))
NextLetter = handle NextLetter(handle)
OnKey = void OnKey(handle, string, Variant(Optional))
OnRepeat = void OnRepeat(handle, string, string)
OnTime = void OnTime(handle, Variant, string, Variant(Optional))
OnUndo = void OnUndo(handle, string, string)
Quit = void Quit(handle)
RecordMacro = void RecordMacro(handle, Variant(Optional))
RegisterXLL = bool RegisterXLL(handle, string)
Repeat = void Repeat(handle)
SaveWorkspace = void SaveWorkspace(handle, Variant(Optional))
Undo = void Undo(handle)
Volatile = void Volatile(handle, Variant(Optional))
Wait = bool Wait(handle, Variant)
GetPhonetic = string GetPhonetic(handle, Variant(Optional))
CalculateFull = void CalculateFull(handle)
FindFile = bool FindFile(handle)
FileDialog = handle FileDialog(handle, MsoFileDialogType)
CalculateFullRebuild = void CalculateFullRebuild(handle)
CheckAbort = void CheckAbort(handle, Variant(Optional))
DisplayXMLSourcePane = void DisplayXMLSourcePane(handle,
Variant(Optional))

CalculateUntilAsyncQueriesDone = void
CalculateUntilAsyncQueriesDone(handle)
SharePointVersion = int32 SharePointVersion(handle, string)
methods(mwsamp)
Methods for class COM.MWSAMP_MwsampCtrl_2:

AboutBox GetVariantArray addproperty
AddDouble GetVariantVector constructorargs
Beep Redraw delete
FireClickEvent RetErrorInfo deleteproperty
FireEventArgs ReturnVTError events
FireMouseDownEvent SetBSTR get
Fire_Double_Click SetBSTRArray interfaces
GetBSTR SetI4 invoke
GetBSTRArray SetI4Array load
GetI4 SetI4Vector move
GetI4Array SetIDispatch propedit
GetI4Vector SetR8 release
GetIDispatch SetR8Array save
GetR8 SetR8Vector send
GetR8Array ShowVariant set
GetR8Vector VariantOfTypeHandle

Calling methods of COM objects can be done in one of the following ways:
Using the INVOKE function
hExcelWorkbooks = get(hExcel,'Workbooks');
hExcelw = invoke(hExcelWorkbooks, 'Add');
Using the method name
hExcelRange = Range(hExcel,'A1:D4');
set(hExcelRange,'Value',rand(4));
Passing Arguments by Reference
Certain COM Objects expose methods with arguments that are also used as output.
This is referred to as by-reference argument passing. In MATLAB, this is achieved
by sending the output as the return from calling the method.

The GetFullMatrix method of a MATLAB Automation server is an example of a
COM method that accepts arguments by reference. This example illustrates how
passing arguments by reference is achieved in MATLAB.
% Register MATLAB session as the automation server version.
regmatlabserver;

hmatlab = actxserver('matlab.application.single')

hmatlab =

COM.matlab_application_single

invoke(hmatlab)
GetFullMatrix = [SafeArray Pointer(double), SafeArray Pointer(double)]
GetFullMatrix(handle, string, string, SafeArray Pointer(double), SafeArray
Pointer(double))
PutFullMatrix = void PutFullMatrix(handle, string, string,
SafeArray(double), SafeArray(double))
Execute = string Execute(handle, string)
MinimizeCommandWindow = void MinimizeCommandWindow(handle)
MaximizeCommandWindow = void MaximizeCommandWindow(handle)
Quit = void Quit(handle)
GetCharArray = string GetCharArray(handle, string, string)
PutCharArray = void PutCharArray(handle, string, string, string)
GetWorkspaceData = Variant(Pointer) GetWorkspaceData(handle, string,
string)
PutWorkspaceData = void PutWorkspaceData(handle, string, string,
Variant)
Feval = Variant(Pointer) Feval(handle, string, int32, Variant(Optional))
GetVariable = Variant GetVariable(handle, string, string)
get(hmatlab)
Visible: 1

Interact with the MATLAB running as an Automation server using the
PutFullMatrix, Execute, and GetFullMatrix methods.
hmatlab.Execute('B2 = round(100*rand(1+round(10*rand)))');

In the next step, you can determine the size of the array to get from the MATLAB
Automation server without needing to check manually.
Execute(hmatlab,'[r,c] = size(B2); B2_size = [r,c];');
[B_size, z_none] = GetFullMatrix(hmatlab,'B2_size','base',[0 0],[0,0]);
Since the size has been determined, you can just get the B2 data using the
GetFullMatrix method.
[B, z_none] = GetFullMatrix(hmatlab,'B2','base',zeros(B_size),[0,0])
B =

91 16 96 66 32 45 50 89 24
13 97 66 17 95 65 96 96 93
91 96 4 71 3 71 34 55 35
63 49 85 3 44 75 59 14 20
10 80 93 28 38 28 22 15 25
28 14 68 5 77 68 75 26 62
55 42 76 10 80 66 26 84 47
96 92 74 82 19 16 51 25 35
96 79 39 69 49 12 70 81 83

z_none =

0 0

delete(hmatlab)
Event Handling
Events associated with Automation controls can be registered with event handler
routines, and also unregistered after the Automation control object has been
created in MATLAB.
events(hExcel)
NewWorkbook = void NewWorkbook(handle Wb)
SheetSelectionChange = void SheetSelectionChange(handle Sh, handle
Target)
SheetBeforeDoubleClick = void SheetBeforeDoubleClick(handle Sh,
handle Target, bool Cancel)

SheetBeforeRightClick = void SheetBeforeRightClick(handle Sh, handle
Target, bool Cancel)
SheetActivate = void SheetActivate(handle Sh)
SheetDeactivate = void SheetDeactivate(handle Sh)
SheetCalculate = void SheetCalculate(handle Sh)
SheetChange = void SheetChange(handle Sh, handle Target)
WorkbookOpen = void WorkbookOpen(handle Wb)
WorkbookActivate = void WorkbookActivate(handle Wb)
WorkbookDeactivate = void WorkbookDeactivate(handle Wb)
WorkbookBeforeClose = void WorkbookBeforeClose(handle Wb, bool
Cancel)
WorkbookBeforeSave = void WorkbookBeforeSave(handle Wb, bool
SaveAsUI, bool Cancel)
WorkbookBeforePrint = void WorkbookBeforePrint(handle Wb, bool
Cancel)
WorkbookNewSheet = void WorkbookNewSheet(handle Wb, handle Sh)
WorkbookAddinInstall = void WorkbookAddinInstall(handle Wb)
WorkbookAddinUninstall = void WorkbookAddinUninstall(handle Wb)
WindowResize = void WindowResize(handle Wb, handle Wn)
WindowActivate = void WindowActivate(handle Wb, handle Wn)
WindowDeactivate = void WindowDeactivate(handle Wb, handle Wn)
SheetFollowHyperlink = void SheetFollowHyperlink(handle Sh, handle
Target)
SheetPivotTableUpdate = void SheetPivotTableUpdate(handle Sh, handle
Target)
WorkbookPivotTableCloseConnection = void
WorkbookPivotTableCloseConnection(handle Wb, handle Target)
WorkbookPivotTableOpenConnection = void
WorkbookPivotTableOpenConnection(handle Wb, handle Target)
WorkbookSync = void WorkbookSync(handle Wb, Variant
SyncEventType)
WorkbookBeforeXmlImport = void WorkbookBeforeXmlImport(handle
Wb, handle Map, string Url, bool IsRefresh, bool Cancel)
WorkbookAfterXmlImport = void WorkbookAfterXmlImport(handle Wb,
handle Map, bool IsRefresh, Variant Result)
WorkbookBeforeXmlExport = void WorkbookBeforeXmlExport(handle
Wb, handle Map, string Url, bool Cancel)
WorkbookAfterXmlExport = void WorkbookAfterXmlExport(handle Wb,
handle Map, string Url, Variant Result)

WorkbookRowsetComplete = void WorkbookRowsetComplete(handle
Wb, string Description, string Sheet, bool Success)
AfterCalculate = void AfterCalculate()
The following command registers five of the supported events for MWSamp2 to
the event handler, e_handler.m.
dbtype e_handler.m 1:3
1 function e_handler(varargin)
2
3 disp(['Event ',varargin{end},' triggered!!!'])

registerevent(mwsamp, {'Click' 'e_handler';...
'DblClick' 'e_handler';...
'MouseDown' 'e_handler';...
'Event_Args' 'e_handler'})
eventlisteners(mwsamp)
ans =

'Click' 'e_handler'
'DblClick' 'e_handler'
'MouseDown' 'e_handler'
'Event_Args' 'e_handler'

Another way of doing this would be to first register all the events, and then
unregister the events that are not needed. First, restore the Automation control to
its original state before any events were registered.
unregisterallevents(mwsamp)
ans =

{}

Now register all the events that this COM object supports to the event handler,
e_handler.m.
registerevent(mwsamp,'e_handler')

ans =

'Click' 'e_handler'
'MouseDown' 'e_handler'
'Event_Args' 'e_handler'

Next unregister any events you will not be needing.
unregisterevent(mwsamp,{'Event_Args' 'e_handler';...
'MouseDown' 'e_handler'})
ans =

'Click' 'e_handler'

Error Handling
If there is an error when invoking a method, the error thrown shows the source, a
description of the error, the source help file, and help context ID, if supported by
the COM Object.
set(hExcelw,'Saved',1);
invoke(hExcelWorkbooks,'Close')
try
Open(hExcelWorkbooks,'thisfiledoesnotexist.xls')
catch e
disp(e.message)
end
Invoke Error, Dispatch Exception:
Source: Microsoft Office Excel
Description: 'thisfiledoesnotexist.xls' could not be found. Check the spelling of the
file name, and verify that the file location is correct.

If you are trying to open the file from your list of most recently used files, make
sure that the file has not been renamed, moved, or deleted.
Help File: C:\Program Files (x86)\Microsoft
Office\Office12\1033\XLMAIN11.CHM

Help Context ID: 0
Destroying COM Objects
COM objects are destroyed in MATLAB when the handle to the object or the
handle to one of the object's interfaces is passed to the DELETE function. The
resources used by a particular object or interface are released when the handle of
the object or interface is passed to the RELEASE function.
By displaying the contents of the MATLAB workspace using the WHOS
command, you can observe the COM object and interface handles before and after
using the RELEASE and DELETE functions.
whos mwsamp hExcel

hExcel 1x1 COM.excel_application
mwsamp 1x1 COM.MWSAMP_MwsampCtrl_2

release(hExcelw)
whos mwsamp hExcel

hExcel 1x1 COM.excel_application
mwsamp 1x1 COM.MWSAMP_MwsampCtrl_2

Quit(hExcel)
delete(hExcel);
delete(mwsamp);
close
whos mwsamp hExcel

hExcel 1x1 handle
mwsamp 1x1 handle

Matlab Notes

Uploaded by

Copyright:

Available Formats

Matlab Notes

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Matlab Notes

Uploaded by

Copyright:

Available Formats

Overview of the MATLAB Environment:

MATLAB is a high-level technical computing language and interactive

platforms, start the MATLAB program by double-

platforms, start MATLAB by double-clicking the MATLAB

platforms, start MATLAB by typing matlab at the operating system

You might also like