The Mesmerizing World

Of

MATLAB
Second Edition

Sagren Munsamy
For all those students who bothered me with
M ATLAB questions. . . and those that hate the
inventors of M ATLAB.
The Mesmerizing World of M ATLAB:
An Introduction

Second Edition

Sagren Munsamy

c
Copyright °Sagren Munsamy 2005

All rights reserved. No part of this publication may be reproduced, stored in a retrieval
system, or transmitted in any form or by any means, electronic, mechanical, photocopy-
ing, recording or otherwise, without the prior written permission of the author.

2
Contents

I Introducing M ATLAB 9
1 The Essentials of m-file Programming 13
1.1 The Basic Structure of an m-file . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 Variable Naming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.3 Data Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.4 Assigning Values to Variables . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.5 Matrix, Scalar and Array Operations . . . . . . . . . . . . . . . . . . . . . 16
1.6 Flow Control: Loops . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.7 Flow Control: Decision Structures . . . . . . . . . . . . . . . . . . . . . . . 18
1.8 Essential M ATLAB Functions . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2 Input and Output: The Essentials and More. . . 23

2.1 Basic Input and Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.1.1 Using disp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.1.2 Getting User Input: input . . . . . . . . . . . . . . . . . . . . . . . 24
2.1.3 A Word on Prompts . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.2 Advanced Input and Output . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.1 Basic Uses of fprintf . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.2 Command Window Output: Creating Tables of Data . . . . . . . . 25
2.2.3 Exporting and Importing Data: Text Files . . . . . . . . . . . . . . . 26
2.3 Exporting Numeric Data Only . . . . . . . . . . . . . . . . . . . . . . . . . 28

3 Basic Data Visualisation: Graphics 31

3.1 plot and plot3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Log Plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3 surf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.4 Multiple Plots on a Figure . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.5 Axis Limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.6 Adding Legends to Graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.7 Annotating Graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

II Intermediate Programming 39
4 Cell Arrays and Structures 43
4.1 Structure or Cell Array? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.2 Cell Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3
 4.2.1 Creating Cell Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.2 Referencing Cell Array Items . . . . . . . . . . . . . . . . . . . . . 44
4.2.3 Removing and Clearing Contents of Cell Arrays . . . . . . . . . . . 45
4.2.4 Reshaping Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.2.5 Nested Cell Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.2.6 Visualising Cell Arrays . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.3 Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.3.1 Creating Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.3.2 Types of Arrangements . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.3.3 Dynamic Fields: A Way to Search for Data . . . . . . . . . . . . . . 49
4.3.4 Setting Up a Searchable Structure . . . . . . . . . . . . . . . . . . . 49
4.3.5 Removing Fields from Structures . . . . . . . . . . . . . . . . . . . 50
4.3.6 Operations on Structures . . . . . . . . . . . . . . . . . . . . . . . . 50
4.3.7 Nested Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

5 Logical Subscripting 53
5.1 Logical Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.2 Applications of Logical Vectors . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.2.1 The find command . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.2.2 Discontinuous Graphs . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.2.3 Avoiding Division by Zero . . . . . . . . . . . . . . . . . . . . . . . 56
5.2.4 Avoiding Infinity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.2.5 Counting Random Numbers . . . . . . . . . . . . . . . . . . . . . . 57
5.3 Replacing Stacked elseif Ladders . . . . . . . . . . . . . . . . . . . . . . 57

6 Function Handles, feval, inline & Function Overloading 59

6.1 The inline Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.2 The Requirements of Iterative Procedures and Symbolic Variables . . . . . 61
6.3 Function Handles and feval . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.4 Making Subfunctions Global . . . . . . . . . . . . . . . . . . . . . . . . . . 63

III Creating GUIs 65

7 Creating GUIs in M ATLAB 69
7.1 M ATLAB GUI Components . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
7.2 Getting Started . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
7.3 An Easy Example: Sum of Two Numbers . . . . . . . . . . . . . . . . . . . 71
7.4 UIControls: Useful Properties and Related Stuff . . . . . . . . . . . . . . . 77
7.4.1 Push Button . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
7.4.2 Radio Button . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
7.4.3 Sliders . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
7.4.4 List Box . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
7.4.5 Pop Up Menu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7.4.6 Axes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7.4.7 Accessing Workspace Variables . . . . . . . . . . . . . . . . . . . . 79
7.5 A More Complex Example . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

4
IV Numerical Methods 87
8 Basic Integration and Differentiation, and Data Regression 91
8.1 Why Numerical Methods? . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
8.2 Basic Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . . . 92
8.3 Basic Numerical Differentiation . . . . . . . . . . . . . . . . . . . . . . . . 94
8.4 Data Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
8.5 When the Derivative Does Not Exist at a Point. . . . . . . . . . . . . . . . . 96

9 First Order Differential Equations and Runge-Kutta Methods 99

9.1 Euler’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
9.2 Vector-valued Functions and Systems of Differential Equations . . . . . . 102
9.3 Using M ATLAB Built-in ODE Solvers: ODE45 . . . . . . . . . . . . . . . . . 103
9.4 Using the events Option of ODE Solvers . . . . . . . . . . . . . . . . . . 105

5
 Acknowledgements

I would like to thank

1. Dr. Randhir Rawatlal for answering my M ATLAB related questions, and for his
excellent notes An Introduction to Programming in M ATLAB With Applications in
Chemical Engineering prepared in April 2003. The chapter in this book Basic Inte-
gration and Differentiation, and Data Regression is almost the same as the equivalent
sections in Dr. Rawatlal’s notes, albeit that the code is presented differently here.
2. Brian D. Hahn for his very good book Essential M ATLAB for Scientists and Engi-
neers.
3. All those unknown authors that posted M ATLAB related information and tutori-
als on the Internet.
4. Donald E. Knuth (the creator of TEX) and Leslie Lamport(the creator of LATEX 2ε ),
without whom I would not be able to write this book. This entire book is typeset
in LATEX.
5. KS (this person wanted to remain anonymous!) who read the final draft of this
book on the day before New Year’s Eve of 2005 and made valuable suggestions
as to what should be included in a possible future edition of this book.
6. . . . and of course, all those students who unknowingly helped me decide what to
include in this book by asking me many questions.

6
 Preface to the Second Edition

This book was made available to the students of the University of KwaZulu-Natal
during the first semester of 2005. Some students found the book useful, while others
were not aware that the book was available. After having had a chance to gain some
feedback on the material in the book, I decided to review some parts of the book. The
changes in the book are minor. Typographical errors were corrected and so were some
of the programs. A new cover was also designed.
The intended audience for the book remains the same . . . those that are beginning
M ATLAB, and intermediate M ATLAB programmers.
I sincerely hope that this book proves to be useful.

– Sagren Munsamy
[email protected]
BScEng(Chemical), 3rd Year
July 2005

7
 Preface to the First Edition

This book took just under one and a half months to complete. My vision was to
write a book that would serve as a simple introduction to using M ATLAB as a computa-
tional environment. The end-result and my initial vision are totally different. This book
quickly covers the basics of m-file programming, and then goes on to introduce more
advanced topics such as cell arrays and structures. I found myself including information
that was not so basic, very early in the writing of this book. I have tried to keep the
language and explanations as simple as possible. I intended including a section with
examples and exercises at the end of this book, but time constraints have not made that
possible. In a future edition, it shall be included.
This book is not a true introductory guide to using M ATLAB. The aim of this book
is to provide you with sufficient information, tools and explanations to broaden your
knowledge of M ATLAB programming. Some experience of using M ATLAB is thus re-
quired.
In my opinion, this book and Brian D. Hahn’s book Essential Matlab for Scientists
and Engineers complement each other quite well. His book assumes no knowledge of
M ATLAB programming and also has a greater volume of information regarding basic
m-file programming. Graphics is not covered in detail in this book, but it is in Hahn’s.
The very interesting topic of logical vectors is presented very differently here, albeit that
the same examples have been used for that section. Some of the topics that have been
omitted will be posted as separate articles on the Engineering Server at UKZN. The path
on the server is S:\COURSES\mecheng\matlab book.
I have attempted to explain numerical methods and the use of ODE solvers in more
detail than that covered in Hahn’s book. I think you will find it easier to understand the
sections in this book. GUIs, structures and cell arrays have also been covered in greater
detail in this book.
This book alone will not be sufficient to develop your M ATLAB programming skills.
You need to practice, and to practice you need examples. There are many tutorials
available on the Internet and on the engineering server, and great questions in Hahn’s
book.
I do not by any means consider myself to be a great M ATLAB programmer, or some-
one that knows a lot of M ATLAB. I have just tried to write down some of the stuff that
I think will prove useful to undergraduate engineering students. It is my sincere hope
that this book will help unleash your broaden your repertoire of programming skills
and aid in the development of structured thinking and prove useful to you.

– Sagren Munsamy
[email protected]
BScEng(Chemical), 3rd Year
December 2004

8
 Part I

Introducing M ATLAB

9
 This part introduces the basics of m-file programming and also includes some ad-
vanced programming techniques and information. Input, output and displaying of data
(including graphics) is covered.

11
12
Chapter 1

The Essentials of m-file

Programming

This chapter will deal with

• The basic structure of an m-file.
• Naming of variables
• Matrix, scalar and element-by-element (array) operations
• Flow control in m-files
• Decision structures
• Essential M ATLAB functions

Those that find this introductory chapter useful have KS(see the preface) to thank,
as this individual told me in a very nice manner that this book would not be complete
without this chapter. There are many great introductory M ATLAB tutorials around and
that is why I did not initially want to include this. I will do my best to introduce the
basics of M ATLAB in this chapter, but this introduction is by no means a definitive intro-
duction. I still assume that you have had some experience using the M ATLAB command
window.

1.1 The Basic Structure of an m-file

M ATLAB m-files are so named because they have an extension of .m. There are two
types of m-files viz., script m-files and function m-files. Script m-files are just collections
of M ATLAB statements. A function m-file has a definition line and function body1 . The
structure is thus as follows:
1 Function m-files can also have multiple functions in a single m-file. The first function is then

called the primary function and the rest are called subfunctions. This matter is dealt with in an-
other chapter.

13
function outputargs = function_name(inputargs)

FUNCTION BODY

outputargs is the variable via which the function returns results, and inputargs
is the variable via which the function receives inputs. For those that have done some
other programming language–inputargs is also known as a parameter. You can have
multiple input and/or output arguments in the function definition line. Examples of
function definition lines are:

function meanvalue=mv(x,y)

function [meanvalue,sum]=mv(x,y)

As you can see, multiple output arguments are enclosed by square brackets and
multiple input arguments are enclosed by round brackets. In both cases, multiple ar-
guments are separated by commas. Also note that the m-file name should be the same
as the name of the function in the case of a function m-file. For the first example of a
function definition line, the name of the m-file should be mv.m
Of course, if one has many arguments, then it becomes an irritation to use the above
notation for multiple arguments. In that case, one can resort to cell arrays or structures.
Since it is quite an advanced topic, it is presented later on in this book.
To access the function mv, we need only type, provided that a and b are defined in
the base workspace

mv(a,b); %returns the resuts in the variable "ans"

mean=mv(a,b) % stores the result in "mean"

As you can see, the arguments passed to the function need not have the same name
as the input arguments of the function.

1.2 Variable Naming

Type the following in the command window,

x=[0 1 2 3];
y=[0 2 4 6];
plot(x,y);

You should get a plot of a straight line.

Now type,

close all;
plot=6;
plot(x,y);

14
 This time, you will not get a plot of a straight line.M ATLAB has assigned the numeric
value of 6 to the function name(which it now thinks is a variable!). So what does one
do if one accidently assigns a value to an already defined function? All you need to do
is type clear all to clear all variables in the workspace and restore everything to its
normal state. Variable naming problems can be sorted out by,

1. Typing the proposed variable name in the command window. If M ATLAB returns

???undefined function or variable ’proposed_name’

where proposed name represents the proposed variable name; and step 2 re-
turns the value one, then the variable name can be used.
2. type isvarname proposed name

1.3 Data Types

There are 15 fundamental data types in M ATLAB. Since M ATLAB is matrix-based, it
treats all data types as arrays2 . In this book we shall concern ourselves with:

• logical arrays
• character arrays3
• numeric data type double
• cell arrays
• structures
• function handles

In this introductory chapter we will restrict ourselves to numeric data type double.

1.4 Assigning Values to Variables

Assignment of values takes place as one would expect.

EXAMPLES

• assign the value 2 to x

x=2
µ ¶
1 2
• assign the matrix to x
3 4
x=[1 2; 3 4]
or alternatively,
x=[1,2; 3,4]

To prevent the assignment result from being echoed in the command window, state-
ments should be terminated with semicolons.
2 M ATLAB refers to vectors, matrices and n-dimensional matrices as arrays.
3 Character arrays are dealt with in the chapter on string handling

15
1.5 Matrix, Scalar and Array Operations
Matrix Operations M ATLAB interprets all operators in a matrix-sense unless oth-
erwise specified by the user. In other words, * and / mean multiply and divide in the
matrix sense. This of course means that the dimensions of the matrices that are involved
must agree. Matrix operations are covered in Applied Math I(a), so I will just present
matrix multiplication and matrix division.

TOOL 1.1 Let A ∗ B = C where A, B and C are matrices with dimensions m-by-n, n-by-m
and m-by-m respectively. The matrix product of A and B is defined by the matrix C where,

n X
X n

Cij = Aij Bji

i=1 j=1

TOOL 1.2 Let A be an m-by-n matrix; x and b be of dimension n-by-1. If Ax = b then,

x = A−1 b 6= bA−1

In M ATLAB 4 this can be represented as either x=inv(A)*b or A\b

Scalar Operations Scalar operations are the same as scalar operations in the matrix
sense i.e., each element is scaled by the scalar. An example would be x=3*A

Array Operations Array operations can be considered as being element by element

operations. An example would be if we wanted to multiply the elements of two ma-
trices and store the result in another matrix. If the dimension of one of the matrices is
m-by-n, then the dimension of the other two matrices is necessarily m-by-n. Note the
difference from matrix multiplication. Note that addition and subtraction are defined
as being element-by-element operators only i.e., nothing special happens to these op-
erators when matrices are involved. Array operations are indicated by prefixing the
operator(either * or /) by a decimal point.

TOOL 1.3 Let C = A ∗ B. If this is an array operation, then the dimension of A, B and C is
m-by-n. Then,
Cij = Aij Bij
This can be represented in M ATLAB as C=A.*B
Element-wise division is defined analogously.

4 The latter is often said as “A under b”, to distinguish it from “A over b”. Remember that b is

brought underneath A, therefore in M ATLAB, b appears at the bottom.

16
1.6 Flow Control: Loops
Horror Story I remember the first experience I had with loops in M ATLAB. It was
quite funny, but at the time it did not seem amusing. We had an assignment to hand
in and I had finished the assignment and was confident that everything was logically
correct, and that the syntax was also correct. But, when I ran my program I got the error
??? Attempt to execute SCRIPT for as a function.
I was horrified. I went over the code a million times and could find nothing wrong.
I went to one of the lecturers (I won’t say who–it was not Dr. Rawatlal) and he said that
there was nothing wrong. I did not know what to do. I retyped the code and it worked.
I compared the two m-files and saw that the only difference was one letter—the “f” in
“for” was capitalized. It was supposed to be a small letter! That was what resulted in
the error. From then on, I never used capital letters in programming. I use underscores
if required e.g., rather than MeanValue I use mean value, for purposes of readability.
Anyway. . .

Loop Structures M ATLAB has two loop structures, for and while. The syntax for
for is

for loop_counter=start_val:increment:end_val

STATEMENTS

end;

It is rather similar to PASCAL in that the loop is terminated with an end, but dis-
similar in the sense that it does not have a do. Note that the increment is optional. If it
is not specified, then it is assumed to be one. For example to run a loop from 1 to 10 in
increments of 2, we would say,

for j=1:2:10

blah blah blah ....

end;

Coming from a PASCAL background, I am obsessed with i and j as counters in

definite looping structures. There is some sort of history behind why these two variables
are used, but I don’t recall what. Anyway, try not to use both i and j in a program, as
both represent the i in complex numbers.
One can exit a loop prematurely by using break. If there are nested loops, then
break exits to the nearest outer loop. By the way, m-files can also be terminated pre-
maturely by using return. return returns to the caller of the m-file that is currently
active. Let’s say a was called from the command window and b from a. Having return
in a will make M ATLAB return to the command window, and having return in b will
make M ATLAB return to a. There is also continue, which skips the body of the loop
and resumes the next iteration.
The syntax of a while loop is,

17
while decision

STATEMENTS

end

For example to display “hello” while n is not equal to 4,

while n˜=4
disp(’hello’)
end;

The boolean operators not, and and or are represented by ˜, & and | respectively.
The first symbol is called a tilde and the last one is usually shift+backslash (for each |)
on computer keyboards. Note that these are element-wise operators, and are the ones
that are usually intended when thinking in terms of Boolean algebra.

1.7 Flow Control: Decision Structures

Once again there are two structures. These are if and switch. The syntax for if is,

if decision

STATEMENTS

else

end

The else part is optional, but each if needs a corresponding end. The syntax for
switch is,

switch switch_expr
case case_expr,

STATEMENTS

case {case_expr1, case_expr2, case_expr3,...}

STATEMENTS

otherwise

STATEMENTS

end

The otherwise part is optional. Note that case expr and switch expr can be
either a string or scalar, making switch quite restrictive. Another important thing to

18
note is that switch short circuits. In other words, if one of the case expressions is true,
then those corresponding statements are executed, and the remainder of the switch
structure is bypassed. If none of the case expressions are true but the otherwise is
true, then the statements for the otherwise part are executed and program execution
is resumed after the end of the switch structure.

Implications of Stacking and Nesting If you have stacked if’s then it could end
up looking like,

if expr1
STATEMENTS
end
if expr2
STATEMENTS
end
if expr3
STATEMENTS
end

This is termed a stacked if ladder and implies that the decisions are independent of
each other. This sort of code can be replaced by,
if expr1
STATEMENTS
elseif expr2
STATEMENTS
elseif expr3
STATEMENTS
end

It looks neater and still implies that the decisions are independent of one another. The
else part can be incorporated into the code with elseif, by putting the else in its
usual place—before the end.
Note that the primary if of nested if’s cannot be turned into elseif’s as the
nested structure implies dependence of the decisions.

1.8 Essential M ATLAB Functions

For those students that are using this book as an aid for the Applied Computer Methods
course at UKZN, it is important that you know how use certain very important func-
tions. Most of your tasks will use these functions in some way. These functions and
methods are,
1. transpose
2. linspace
3. vector initialization
4. repmat
5. size

19
transpose A matrix can be transposed by using the ’ operator, or the transpose
function. transpose(a) is the same as a’, provided that a is non-complex. transpose(a)
is the pure transpose of a and is always the same as a.’.

linspace This function is used to generate linearly spaced points over an interval.
The general syntax is var=linspace(start val,end val,numpoints). Note that
the number of points is specified, not the step size.

Vector initialization Here you specify the step size and the start and end values. The
general syntax is var=start val:step size:end val. Let a represent the starting
value, b the end value, h the step size and n the number of points. The relationship
amongst these variables is h = (b − a)/(n − 1). It is quite useful to keep this in mind.

repmat This function replicates and tiles a matrix or vector. For example let’s say
that x=[1 2 3]. If we say,

newx=repmat(x,[4 1])

the value of newx is,

1 2 3
1 2 3
1 2 3
1 2 3

So the general syntax is repmat(VectorToCopy,[Rows Columns]).

The major use for this is to generate linear combinations of variables. Let’s say that
x varies from 1 to 10 and y varies from 100 to 200, and we want to see how z = x./y
behaves over the specified ranges for all possible combinations of x and y. We proceed
as follows:
1. Set one vector up as a column vector, and the other as a row vector. The choice is
arbitrary. We shall choose x to be set up as rows i.e., a column vector.

x=linspace(0,10,100)’ %take 100 points arbitrary choice

y=linspace(100,200,60) %could have chosen the same number of points

2. For the variable that has been set up in rows, replicate in columns i.e., the first
argument of the square bracket in repmat should be 1 and the other should be
n, where n is the number of elements in the other vector. This is done to make
sure that the matrix dimensions agree. For the variable that has been set up in
columns replicate in rows.

newx=repmat(x,[1 60])
newy=repmat(y,[100 1])

3. Evaluate the function based on the new matrices

z=x./y %don’t forget the dot!
After doing this we may plot the data as a surface by saying surf(newx,newy,z),
or plot x versus y by plot(x,y).
An alternative to using repmat is meshgrid. M ATLAB help states:

20
 MESHGRID X and Y arrays for 3-D plots.
[X,Y] = MESHGRID(x,y) transforms the domain specified by vectors
x and y into arrays X and Y that can be used for the evaluation
of functions of two variables and 3-D surface plots.
The rows of the output array X are copies of the vector x and
the columns of the output array Y are copies of the vector y.

[X,Y] = MESHGRID(x) is an abbreviation for [X,Y] = MESHGRID(x,x).

[X,Y,Z] = MESHGRID(x,y,z) produces 3-D arrays that can be used to
evaluate functions of three variables and 3-D volumetric plots.

For example, to evaluate the function x*exp(-xˆ2-yˆ2) over the

range -2 < x < 2, -2 < y < 2,

[X,Y] = meshgrid(-2:.2:2, -2:.2:2);

Z = X .* exp(-X.ˆ2 - Y.ˆ2);
mesh(Z)
The code for the example presented can be simplified to,

%100 points in x and 60 in y

[x,y]=meshgrid(0:10/(100-1):10,100:(200-100)/(60-1):200);
z=x./y;
surf(x,y,z); %if you want a surface
shading interp; %makes the surface look pretty!

Note that shading interp is not required, and that the number of points have been
specified in meshgrid. As you can see, the code is much simpler, and you don’t have to
worry about the dimensions of x and y. A word of warning . . . round–off errors can lead
to the end point not being included, if the number of points in an interval are specified,
but you decide to specify the interval, as in the meshgrid statement above. This issue
is not as serious as I make it sound, though.

size This function returns the number of rows and columns of an one-dimensional
or two dimensional array, as a matrix of size 1-by-2. If you want to know the number of
columns there are in a matrix then you can execute size(m,2). Alternatively, to find
the number of rows, you can execute size(m,1). To find the number of elements in a
one-dimensional array (let’s say x), you can also use length(x). length is useful in
string handling.

21
22
Chapter 2

Input and Output: The

Essentials and More. . .

This chapter will present

• input and output using input and disp
• controlling command window output using fprintf
• exporting data to text files using fprintf
• exporting data to Microsoft EXCEL using fprintf
• inputting data into M ATLAB from text files

Most people will be happy knowing how to output and input data using disp and
input, respectively. This chapter will deal with the use of the basic input and output
commands and the more versatile(thus, more complicated) fprintf command. You
may consider Graphical User Interfaces as being advanced input/output structures1 , but
since there is so much of information that one needs to know about GUIs, a separate
chapter is devoted to the creation and use of GUIs.

2.1 Basic Input and Output

2.1.1 Using disp
disp displays a string or array, together with a linefeed at the current line in the com-
mand window.

⊲ E XAMPLE
Let x=[1 2 3 4]. Use disp to display,

x(1) is 1
1 In this instance structures do not refer to M ATLAB structures, although structures are used

widely in GUI applications.

23
Solution
Since disp only accepts a single argument, an appropriate argument needs to be cre-
ated. The code is,

x=[1 2 3 4];
s=[’x(1) is ’, num2str(x(1))]
disp(s);

2.1.2 Getting User Input: input

The general syntax is,

i=input(’prompt’) % for numeric input

i=input(’prompt’,’s’) % for string input

Note that when the optional argument s is used, the user does not have to enclose a
string input within quotes. Quotes are necessary when the optional argument is omit-
ted.

2.1.3 A Word on Prompts

I have have seen many students that use prompts like:

Please enter the value of k

It is bad programming practice. Programs are supposed to be user-friendly. This

does not mean that the prompts have to sound friendly. Prompts need only be appro-
priate and informative. A better version of the above prompt would be,

Enter k:

An even better, but still informative prompt is,

k?

Now if k was, say, thermal conductivity in units of W/mK, we could have some-
thing like,

k (thermal conductivity W/mK)?

It is better than saying,

Please enter the value of thermal conductivity in units of W/mK

Long prompts irritate people as most people do not like reading much! I know that
I have found it irritating to debug a program that requires more than 5 inputs, with
prompts that are a mile long.

24
2.2 Advanced Input and Output
2.2.1 Basic Uses of fprintf
fprintf allows exact control over command window output. It can be used to:
1. control command window output
2. write tables of data to text files

2.2.2 Command Window Output: Creating Tables of Data

fprintf makes it easy to create tables with headings in the command window. The
output is not as pretty as that created in a spreadsheet, but it is good enough. In its
simplest form, the syntax for fprintf is,

fprintf(’format_string’,list_of_variables);

The basic procedure to be followed for displaying data as columns with headings
is,
1. The column headings must be set up as rows i.e., a column vector.
2. The data to be displayed must be set up in rows. Each column of data must be
represented by a separate row.
3. Apply fprintf to the headings then the data.

Format Strings
An example of a format string is,

%-10.2e\n

If you have done some C programming, then this will look very familiar. The format
string always starts with a %. The minus is called a flag and indicates that the data should
be left justified. In output of tables, headings and data should be right justified. To right
justify, omit the minus. The next number specifies the width of the field. It is analogous
to a column width in spreadsheets. The number after the decimal point indicates the
number of decimal places to be included in the output. The e, tells M ATLAB to output
the data in exponential (scientific) notation. g specifies that M ATLAB should choose the
more compact of fixed point and scientific notation. To output a sting using fprintf,
the e is replaced by an s. The \n is called an escape code. The n indicates that the cursor
should be moved to the next line.

⊲ E XAMPLE
Lets say we have the following data and want to output only the run number and
height:
Run Time/s Height/m
1 3 4
2 8 8
3 6 6
4 10 9
We can use the code,

25
r=1:4;
h=[4 8 6 9];
fprintf(’%10s %10s\n’,’Run’,’Height/m’);
fprintf(’%10.2g %10.2g\n’,[r;h]);

This returns the output,

Run Height/m
1 4
2 8
3 6
4 9

Note that each format string is coupled with a corresponding row in the second
fprintf statement i.e., the first format string is coupled with the first row of the matrix
etc.
When using M ATLAB to output tables, it is easy to output only the data that you
want to see. The command window output can then be copied and pasted into EXCEL.
You can then click on the clipboard that appears and select Use Text Import Wizard2 . You
can then set up the data as required. Thereafter, you can make the spreadsheet look
pretty!

2.2.3 Exporting and Importing Data: Text Files

Data can be saved in .mat file format using save and loaded back into M ATLAB using
load. This format is native to M ATLAB. It is useful to know how to save data to text
files, as almost any software package can import text files.
fprintf can also be used to export data to text files. The syntax is very similar to
that used for output of formatted data to the command window. The basic syntax is,

fprintf(file_id,’format_string’,list_of_variables)

fileid is an integer value returned by the function fopen. The procedure to write
data to a text file is:
1. Arrange the data as for output of formatted data to the command window.
2. Use fopen to open a file, and assign the file ID to a variable.
3. Use fprintf to write the data to the file.
4. Use fclose to close the file.
This is best explained with an example.

⊲ E XAMPLE
Let us take the same data used in the previous example and keep the specification the
same. This time we will output the data to a file data.txt3 using fprintf. The code
is,
2 This applies to Office XP and later. If you are using some other version of Office then type Text

Import in the help box.

3 This file will reside in the current directory at the time of running this code. You can of course

change the path by using DOS syntax.

26
r=1:4;
h=[4 8 6 9];
file_id=fopen(’data.txt’,’w’);
fprintf(file_id,’%10s %10s\n’,’Run’,’Height/m’);
fprintf(file_id,’%10.2g %10.2g\n’,[r;h]);
fclose(file_id);

An explanation of fopen is required. The general syntax is,

fopen(’file_name’,’access_mode’);

access mode can be,

• w overwrite existing data and write access only.
• r read access only.
• a append data to end of file and write access only.
• a+ append data to end of file and read/write access.
• r+ read/write access
• w+ overwrite existing data and read/write access.
Note that the access modes that contain a plus, cannot be used as simply as the other
access modes. See M ATLAB help for details.
This text file can then be imported into EXCEL using the Text Import Wizard. I sup-
pose that one of the uses of this technique is that you can do your calculations at at a
computer that has M ATLAB and then use EXCEL to complete the calculations at home.
For the sake of completeness, here’s how to import all the data back into M ATLAB,

[run,data]=textread(’data.txt’,’%s %s’)

This imports the entire columns i.e., header and data as strings, and stores them in
variables run and data. The output is,

run =

’Run’
’1’
’2’
’3’
’4’

data =

’Height/m’
’4’
’8’
’6’
’9’

The headers and data can be retrieved via normal matrix operations, and the strings
can be converted to numerical values by using char and then str2num i.e.,

27
r=run(2:end,1);
d=data(2:end,1);
r=char(r);
d=char(d);
r=str2num(r);
d=str2num(d);

If you are going to use this sort of importing a lot or require automated importing of
data for say, a presentation or post practical interview, I suggest that you write up a
function m-file that can do this for you. It is an interesting task! Of course, you can just
use the import wizard to do this for you.
An easier alternative to textread is importdata. All you need to do in this case
is to specify the file name. The conversion to numerical values is not required. The code
is,

expdata=importdata(’data.txt’);

This returns a structure,

expdata =

data: [4x2 double]

textdata: {’Run’ ’Height/m’}
colheaders: {’Run’ ’Height/m’}

The data can then be retrieved by,

expdata.data

which returns,

ans =

1 4
2 8
3 6
4 9

This is probably the preferred way of programatically import this sort of text file
into M ATLAB.There are tonnes of other functions that are available for the import of
data. Type help fileformats in the command window to see a list.

2.3 Exporting Numeric Data Only

After reading all of this, you must be thinking that there must be an easier way to export
say, columns of data to a text file so that you can use it later (say in another application).
It turns out that there is. Let’s say that we had the set of (arbitrary) data,

28
x y
10 40
30 200
80 400
90 350

Now assume that the first column is stored in a variable x, and the second column
stored in the variable y(these must be set up as column vectors for the code below to
work!). To export these two columns to the file mydata.txt, we first need to set up a
single variable that has the data in the arrangement that we want, and then export the
data i.e.,

data=[x y]
save mydata d -ascii -double -tab

-ascii tells M ATLAB to save the data in ASCII format, -double causes the data
to be save with 16 bit precision, and -tab indicates that the data should be tabulator
delimited. The contents of the file should be, (you may view the file with a text-editor)

1.0000000000000000e+001 4.0000000000000000e+001
3.0000000000000000e+001 2.0000000000000000e+002
8.0000000000000000e+001 4.0000000000000000e+002
9.0000000000000000e+001 3.5000000000000000e+002

The drawback of using save is that you cannot have column headings exported
with the numeric data. You will have to add these in if you require them. Oh yes, the
data can be read into M ATLAB by using load.

29
30
Chapter 3

Basic Data Visualisation:

Graphics

This chapter will cover:

• Plotting graphs.
• Generating surfaces.
• Annotation of graphs using TEX commands.

M ATLAB has many functions that allow excellent visualisation of data. There are far
too many functions, so only those necessary and useful to undergraduate students will
be covered in this chapter. One of the great things about using the graphics functions of
M ATLAB is that the output is computed extremely quickly. It is quite an amazing thing.

3.1 plot and plot3

A good explanation of plot is from the online help for M ATLAB,

PLOT Linear plot.

PLOT(X,Y) plots vector Y versus vector X. If X or Y is a matrix,
then the vector is plotted versus the rows or columns of the matrix,
whichever line up. If X is a scalar and Y is a vector, length(Y)
disconnected points are plotted.

PLOT(Y) plots the columns of Y versus their index.

If Y is complex, PLOT(Y) is equivalent to PLOT(real(Y),imag(Y)).
In all other uses of PLOT, the imaginary part is ignored.

Various line types, plot symbols and colors may be obtained with
PLOT(X,Y,S) where S is a character string made from one element
from any or all the following 3 columns:

31
 b blue . point - solid
g green o circle : dotted
r red x x-mark -. dashdot
c cyan + plus -- dashed
m magenta * star
y yellow s square
k black d diamond
v triangle (down)
ˆ triangle (up)
< triangle (left)
> triangle (right)
p pentagram
h hexagram

For example, PLOT(X,Y,’c+:’) plots a cyan dotted line with a plus

at each data point; PLOT(X,Y,’bd’) plots blue diamond at each data
point but does not draw any line.

PLOT(X1,Y1,S1,X2,Y2,S2,X3,Y3,S3,...) combines the plots defined by

the (X,Y,S) triples, where the X’s and Y’s are vectors or matrices
and the S’s are strings.

For example, PLOT(X,Y,’y-’,X,Y,’go’) plots the data twice, with a

solid yellow line interpolating green circles at the data points.

The PLOT command, if no color is specified, makes automatic use of

the colors specified by the axes ColorOrder property. The default
ColorOrder is listed in the table above for color systems where the
default is blue for one line, and for multiple lines, to cycle
through the first six colors in the table. For monochrome systems,
PLOT cycles over the axes LineStyleOrder property.

PLOT returns a column vector of handles to LINE objects, one

handle per line.

The X,Y pairs, or X,Y,S triples, can be followed by

parameter/value pairs to specify additional properties
of the lines.

If we had, say,

x=[1 2 3 4];
y=[5 6 7 8];

y could be plotted against x by using,

plot(x,y);

32
 If we want a green line with the points plotted as circles1 , we could use,

plot(x,y,’go’)

or alternatively,

plot(x,y,’og’)

If you love using properties, then you could say,

plot(x,y,’linestyle’,’none’,’color’,’g’,’marker’,’o’);

Other properties can be set in this manner. See line in the search tab of help. Of
course, one can also use handles, but that topic is a bit too much for this chapter.
plot3 is the three-dimensional version of plot. It is much the same as plot.
M ATLAB also has easy-to-use plotting functions e.g., ezplot and ezsurf. Look
them up in M ATLAB help. They are easy to use!

3.2 Log Plots

If you are an engineering undergrad, you will have probably used logarithmic plots
quite a bit. If you have tried to do log plots in EXCEL, you would have probably found
that the output is not of an excellent quality. . . and of course you can’t do all those nifty
things you can do in M ATLAB in EXCEL!
To generate a plot with a linear x axis and a logarithmic y axis use,

semilogy(x_values,y_values);

For a linear y axis and a logarithmic x axis use,

semilogx(x_values,y_values);

For logarithmic x and y axes simultaneously use,

loglog(x_values,y_values);

3.3 surf
surf is used to plot surfaces. The general syntax is,

surf(x_values,y_values,z_values);

The arguments are usually two-dimensional arrays. See help for more details.

1 When points are emphasized by using circles, dots etc., these circles, dots etc., are referred to

as markers—they mark the points.

33
3.4 Multiple Plots on a Figure
If you want to plot more than one set of data on the same set of axes, you can use hold.
If we use,

hold on

The current figure is not cleared and further plots are appended to the current figure.
To replace the plots on the current figure, use hold off. You may think of hold as
being an instruction to hold all current plots on the current figure.
When using hold, the axis limits are automatically redefined to fit the data. You
can also have multiple axes on the same figure. This is accomplished by using subplot.
The syntax is,

subplot(r,c,p)

This creates r by c axes on the current figure and makes axes with index p the current
axis. M ATLAB counts cumulatively. Rows are counted, and then M ATLAB moves to the
next column. See find in the index. Indexing is done the same way.

3.5 Axis Limits

M ATLAB automatically calculates axis limits. Sometimes you will want to set the axis
limits manually. This can be done by using axis, whose syntax is,

axes([xmin xmax ymin ymax])

3.6 Adding Legends to Graphs

It is likely that you will want to put a legend on a graph. Look at the following:

x=[1 2 3 4];
y=[5 6 7 8];
plot(x,y,’r’);
plot(x,y,’go’);
legend(’data points’,’line plot’);

This gives,

34
 8
data points
line plot

7.5

6.5

5.5

5
1 1.5 2 2.5 3 3.5 4

How does legend know whether the green plot with markers takes the string data
points or whether the red plot does? We did get the correct result though. It’s not
magic, unfortunately! If we change the code to,

figure; %new figure

x=[1 2 3 4];
y=[5 6 7 8];
plot(x,y,’r’);
hold on;
plot(x,y,’go’);
legend(’line plot’,’data points’);

we get,

35
 8
line plot
data points

7.5

6.5

5.5

5
1 1.5 2 2.5 3 3.5 4

which gives the wrong legend.

legend associates the first string with the last plot. Although this usually works
it is better if you use the children property to make sure that the correct string is
associated with the correct plot. For our example, this can be done by,

plotfig=figure; %new figure

clf; % clear the current figure
x=[1 2 3 4];
y=[5 6 7 8];
redplot=plot(x,y,’r’); %handle returned in "redplot"
hold on;
greenplot=plot(x,y,’go’); %ditto
figure(plotfig); %set the current figure
plotaxis=gca %get current axes
set(plotaxis,’children’,fliplr([redplot greenplot]))
legend(’line plot’,’data points’);

This may seem confusing at first, but it isn’t after you try it on a few plots.

3.7 Annotating Graphs

You can attach a label for any of the coordinate axes by using xlabel, ylabel or
zlabel. A title can be added to a plot by using title. All of these commands take
string arguments.
Let’s say that you have the following label for the x-axis,

α = 1/3β 2/3 y2
How to do this in M ATLAB? Answer:

36
figure; %generate a blank figure
xlabel(’\alpha=1/3\betaˆ{2/3} y_2’,’interpreter’,’tex’)

In M ATLAB Greek letters can be generated by prefixing the English spelling with
a backslash. Other symbols can also be generated. Examples can be found in M ATLAB
help. Other useful things are:

\times multiplication sign

\div division sign
\nabla "upside down delta"
\neq not equal to
\partial partial derivative symbol
\propto proportional symbol
\prime prime or complement
\surd root sign
\ldots ellipses

Note that when using super- and subscripts, all the text that you want to appear
as super- or subscript needs to be surrounded by curly braces. You can also select the
fontsize and fontname for the text labels. For example,

xlabel(’\alpha=1/3\betaˆ{2/3} y_2’,’interpreter’,’tex’, ...

’fontname’,’euclid’,’fontsize’,’10’)

Now I don’t really like typing out the TEX commands by hand, so I use MathType
5.0 to generate the TEX commands that result in the label that I require. It is much easier.
Of course, you could spend a day learning some LATEX 2ε , if you really want to.
You can also place text at any point on an axis by specifying the x and y coordinates
and the string. For example,

text(3,7,’\alpha=1/3\betaˆ{2/3} y_2’,’interpreter’,’tex’)

puts the formula above at coordinate (3,7) on the axis. This of course means that
you need to know how the graph is going to look before issuing such a command, or
at least have some way of calculating the coordinate that you want to place the text at.
Note that the default alignment is center. You can change it. . . and I am sure that you
know how to do that by now!
You can place the same formula interactively at a point by using the mouse. To do
this, use

gtext(’\alpha=1/3\betaˆ{2/3} y_2’,’interpreter’,’tex’)

Try it.
By the way, you can use ginput to get the x and y values of the points you click on
using the mouse. You can do this by,

[x,y]=ginput

To terminate selection of points, press Enter.

That completes the introduction to graphics. I don’t think that you will require any
other advanced graphics functions, but if you do M ATLAB help is always just a few key
presses away!

37
38
 Part II

Intermediate Programming

39
 This part covers the very interesting topics of
1. cell arrays
2. structures
3. logical vectors and logical subscripting
4. function handles
These topics are not essential topics, but familiarity with these topics will make some
programming tasks easier. Function handles are required for the chapter on using M AT-
LAB ODE solvers and structures are useful for designing GUIs.

41
42
Chapter 4

Cell Arrays and Structures

This chapter covers:

• Initializing Cell Arrays and Structures
• Uses of Cell Arrays and Structures

Cell arrays and structures can be considered to be conceptually similar. Both these
data types allow one to store variables of different data types in a single variable. This
cannot be done in matrices—all the elements of the matrix need to be of the same data
type.
Another use of cell arrays is in the programming of function m-files. Sometimes, it is
cumbersome to specify all the input arguments to a function in the function definition
line. A structure can be used to simplify the situation where there are many input
arguments to a function or where there are a variable number of input arguments to a
function.1

4.1 Structure or Cell Array?

As mentioned, these data types are conceptually similar. Cell arrays are a collection
of independent(because the data type of each cell does not depend on the data type of
any other cell in the cell array.) cells. A structure can be considered to be a cell ar-
ray where each cell is a named field. For example if we had a structure student
and wanted to have the name of a student in the structure at level 2, we would say
student.name=’Sagren’. Cell arrays are just variables that can hold any data type.
From this it is apparent that one should use structures when named fields are available.
According to M ATLAB help, cell arrays should be used instead of structures when,
• You need to access multiple fields of data with one statement.
1 Of course, the programmers of M ATLAB have thought of exactly this situation when one

programs function m-files. They have provided us with the reserved variables varargin and
varargout, which are to be used in this situation. This shall be presented later in this chapter.

43
 • You want to access subsets of the data as comma-separated variable lists.
• You don’t have a fixed set of field names.
• You routinely remove fields from the structure.

4.2 Cell Arrays

4.2.1 Creating Cell Arrays
Using assignment statements Cell arrays may be created by explicitly specifying
the index and item to be stored in the cell array. M ATLAB provides two ways of initial-
ising cell arrays in this manner viz.,
• via content indexing
• via cell indexing
Of course, there needs to be some way of distinguishing between cell arrays and normal
arrays2 . Cell arrays use curly brackets in place of the usual square brackets for normal
arrays. Let’s say that we wanted to create a cell array that has, as its first element, the
matrix,
1 2 3
à !
4 5 6
7 8 9
This can be accomplished via,
• c{1,1}=[1 2 3;4 5 6;7 8 9]
• c(1,1)={[1 2 3;4 5 6;7 8 9]}
The first method is an example of content indexing and the second an example of cell
indexing. The choice of method is up to you. I prefer the first method.
I have found cell arrays useful in creating tables of data3 You cannot use a string
array to store distinctly separate strings. M ATLAB treats all the elements of a string
array as a single unit.

Preallocation Just as there are functions to preallocate normal arrays, the cell func-
tion is used to preallocate cells. It’s use is analogous to the usual preallocation functions
zeros and ones. For example, to preallocate a cell array of dimension 3-by-4, one need
only say, c=cell(3,4).

4.2.2 Referencing Cell Array Items

Assume that c(1,1)={[1 2 3;4 5 6;7 8 9]}. To reference the third element in
the third row, we just say c{1,1}(3,4). I hope that referencing elements in this way
is fairly obvious. It is if you evaluate the expression from left to right—the first part
returns a numeric array, and to reference a particular element that array, we would
need to use the more familiar square bracket referencing.
Nested cell arrays4 can also be created using the referencing above.
2 When I refer to normal arrays, I am referring to arrays that consist of item(s) of a single data

type e.g., string arrays and numeric arrays.

3 The chapter on input and output presents this use of cell arrays.
4 Nested Cell Arrays are cell arrays that contain other cells. They are thus a cell array of cells.

44
4.2.3 Removing and Clearing Contents of Cell Arrays
You should have guessed by now that to clear an element of the cell array, you would
do the same thing as you would with a numeric array, except that you would use the
curly bracket notation. For example, assume that a=’what’,’is’;1, [1 2; 3 4].
To clear the contents of the cell that contains what, we would say,

a{1,1}={}

If you are curious as to the accuracy of the information presented in this book, then you
probably have looked at the online M ATLAB help and seen that it says,
You can delete an entire dimension of cells using a single statement. Like
standard array deletion, use vector subscripting when deleting a row or
column of cells and assign the empty matrix to the dimension.
A(cell_subscripts) = []
When deleting cells, curly braces do not appear in the assignment state-
ment at all.
This is also a valid command, but also causes some peculiar things to occur. Let’s say
that a=cell(3,4);a{3,4}=[1 2; 3 4]. M ATLAB outputs,

a =
[] [] [] []
[] [] [2x2 double] []
[] [] [] []

Let’s say that we wanted to delete the first element of the cell. One may be inclined
to say (after reading the M ATLAB help quite quickly and not carefully enough!) that,

a(1,1)=[]

This does not work and results in the error,

??? A null assignment can have only one non-colon index.

If you reread the quote from M ATLAB help, you will notice that it mentions entire di-
mension. When we speak of dimension in matrix terms, we are referring to rows or
columns. Then, entire dimension refers to an entire row or entire column. The error
returned should now make sense—if you specify the row or column number you are
specifying the position of the entire row or column. Okay, now we know that we cannot
use this command as we only want to remove a single element, not an entire dimension.
But what if we said,

a(1)=[]

This does not give an error, but M ATLAB does the following:
1. M ATLAB counts cumulatively starting at one at position (1,1) and removes the
element with the specified index. In other words, for the cell array above, the
counting(count is represented in open-close chevrons) is as follows:

45
 a =
<1> [] <4> [] <7> [] <10> []
<2> [] <5> [] <8> [2x2 double] <11> []
<3> [] <6> [] <9> [] <12> []
2. Since the structure cannot exist as a matrix, M ATLAB reshapes the cell array into
a one dimensional array with the element at count position 1 removed.
Note the subtle distinction between clearing the contents of a cell and deleting a cell.
It is often quite easy to confuse the two. So, be careful.
So then, how do we clear the contents of the first element without reshaping the
array? We can say,

a{1,1}=[]

or alternatively,

a{1,1}={}

The above two statements seem to accomplish the same thing. They do not, tech-
nically. One returns an empty cell in the position (1,1) and the other returns an empty
matrix in position (1,1) i.e., they return two different data types. Another subtlety. I
cannot think of any situations where using one of the above methods would cause a
problem. So so not worry too much about this.

4.2.4 Reshaping Arrays

If you have had some experience in a language like Turbo Pascal, you will have un-
doubtedly come across the situation where a two dimensional array needs to be con-
verted to a one dimensional array. The task is usually accomplished by the use of a
for loop.M ATLAB provides reshape to accomplish this task. Since M ATLAB supports
multidimensional arrays, reshape can be applied to multidimensional arrays as well.
The general syntax of reshape is,

reshape(matrix_to_reshape,size_1st_dim, size_2nd_dim,...)

Note that the reshaped matrix can only have the same number of elements as the
matrix to reshape i.e., the product of the size of the matrix to reshape in each dimension
must be the same as the number of elements in the new matrix. The matrix that is
returned by reshape is of dimension

size_1st_dim * size_2nd_dim * size_3rd_dim ... * size_nth_dim

4.2.5 Nested Cell Arrays

These are easily created by using the cell command or direct assignment. For example,
to create a cell array that is of dimension 4-by-3, within which the item in position (2,3)
is a cell array of 1-by-2, we would say,

a=cell(4,3);
a{2,3}=cell(1,2);

46
 The cell array of dimension 1-by-2 is a nested cell array. We may assign a value of 3
to the element in position (1,2) of the nested cell array by saying,

a{2,3}{1,2}=3;

We could also assign a numeric array to position (1,1) of the nested cell array. Let’s
assume that we want to create a matrix with zeros except for the element in position
(1,3), which will be 100. Further, the size of the numeric array is 10-by-3. This can be
accomplished by,

a{2,3}{1,1}=zeros(10,3);
a{2,3}{1,1}(1,3)=100;

4.2.6 Visualising Cell Arrays

If you have tried working with cell arrays, or have looked carefully at the material pro-
vided thus far, you will notice that M ATLAB does not display the actual contents of a
cell array—instead a summary of the data types contained in a cell array is displayed.
celldisp and cellplot can be used visualise the cell array. To display the contents
of a cell array explicitly, use celldisp. To get a graphical representation of the ar-
rangement and structure of a cell array, use cellplot. Try them out.

4.3 Structures
Structures are particularly useful in cases where you have a set of fixed field names. An
example of such a case is in the creation of a database. In doing calculations for practi-
cals (labs) such as the one in various courses in Chemical Engineering, using M ATLAB
with structures in an appropriate manner will lead to a very elegant, easily editable set
of calculations5 .

4.3.1 Creating Structures

A structure consists of the parent structure name and successive sub layers, called fields.
Successive sub layers are indicated by prefixing the field name with a period i.e., ‘‘.’’.
One would instinctively expect to be able to create structures by using direct assignment
and preallocation. These are presented next.

Good practices In creating structures, it is good to sit down and think about how
you want to arrange the structure and its fields. It is good practice to draw a rough
diagram of the intended structure.

Assignment statements To store the character string Jeffry in the field name of
the parent structure Characters, we would say,

Characters.name=’Jeffry’
5 I have also written an article Practical Applications of Cell Arrays and Strucures: Spreadsheeting

that explains how to go about using M ATLAB to do such calculations, and also includes an actual
calculation done by me.

47
This is implicitly interpreted by M ATLAB as being,

Characters(1).name=’Jeffry’

Of course we may add to the structure. Say we wanted to store the name Kitana
in the appropriate field. We would then just say,

Characters(2).name=’Kitana’

Preallocation Preallocation can be achieved by using struct. I don’t think that

it is absolutely necessary to know how to do this since structures are automatically
enlargened to accommodate data. But for those who are really interested and also for
those that are obsessed with performance considerations. . . an excerpt from M ATLAB
help:
The basic syntax of the struct function is
str_array = struct(’field1’,val1,’field2’,val2, ...)
where the arguments are field names and their corresponding values.
A field value can be:
• A single value, represented by any MATLAB data construct
• A cell array of values.
Note: All field values in the argument list must be of the same scale i.e.,
data type (single value or cell array).
This of course makes the use of struct restrictive, as only one data type can be
used as arguments to the function. Some examples of the use of struct are provided
in M ATLAB help.

4.3.2 Types of Arrangements

Plane-field arrangements In this type of arrangement, each field is kept together
as a single unit.

Element-by-element/distributed arrangement Here the each field is distributed

through the structure.

⊲ E XAMPLE
For example, say we have a database, debtors, which we shall assume to be the name
of the parent structure. We wish to have the name and outstanding balance of the
debtor. We could arrange the data in two ways viz.,
1. We could store all the names in a single field and all the balances in a single field.
Related data in fields share common indices, usually. This is an example of plane-
field arrangement.
2. We could treat debtors as an array, where the first debtor’s data is represented
by debtors(1) and so forth. We could then store the name and outstanding
balance of each debtor as a unit, as an item of the array. For example, we could
say, for the first two debtors,

48
 debtors.name = ’Mike Boswell’
debtors.balance = 100000
debtors(2).name = ’Charlie Chaplin’
debtors(2).balance = 4000

I think that the easiest way to remember the two basic arrangement types is to think
of the database here. In plane-field arrangement, names are stored together as a single
field, and so are account balances. In distributed arrangement, all information regarding
a specific person (say) is stored together.

4.3.3 Dynamic Fields: A Way to Search for Data

It took three readings of the M ATLAB help on this section to eventually realise what was
actually meant by dynamic fields.
If you have done computer programming, then you will have probably done some
sort of database project using a specific programming language. Consider what we have
done so far regarding structures. We have seen that we may reference structures using
indices. In database management, if you want to update the record for a particular
individual, using indices to refer to that particular record becomes a bother if you have
more than three records. It would be nice to enter a search string and update all records
that contain that search string as required. This search string task is one of the things
that need to be programmed when doing a database project, as mentioned above.
Dynamic fields allow us to do something similar to the above. We may specify that
the record for John Boswell be updated by specifying the string John Boswell, pro-
vided that John Boswell is a field name. M ATLAB also necessitates the specification of
the field that is to be made dynamic. The next section explains this in detail.

4.3.4 Setting Up a Searchable Structure

Consider the case where we have a database of people’s names and ages. Dynamic
fields allow us to search for a particular field. This of course means that if we want to
be able to search for a record by name, then the name of each person must be a field.
For example, say we have the following data:

Name Age
Sagren 21
Preshodin 17
Keith 20

To create a structure(assume the parent structure to be student) that is searchable

by name, we have to arrange the structure as,

student.Sagren.age=21
student.Preshodin.age=17
student.Keith.age=20

We may now retrieve any of the above student’s age by using a single function m-
file that has the name of the student as a dynamic field.
The m-file would be something like,

49
function getage(studentstr,name)
disp([name, ’: ’, num2str(studentstr.(name).age), ’ years.’])
To retrieve the age of student Keith, we would say,
getage(student,’Keith’)
which returns
Keith: 20 years.
Note that in database management, the dynamic field is usually unique to each
record. The reason for this is obvious.

4.3.5 Removing Fields from Structures

To remove a field and its all its dependent fields, you can use rmfield. The basic syntax
is,
rmfield(parent_struc_name,’field_to_remove’)
As an example, take the structure student created above. If we said,
rmfield(student,’Sagren’)
The field Sagren and all children of this field are removed. This is important to re-
member, although it may seem trivial.

4.3.6 Operations on Structures

It is very likely, that after setting up a structure you may want to do operations on the
fields in the structure. Operations on fields and structures are usually carried out using
for loops. M ATLAB allows us to do these calculations in a more succinct way.

⊲ E XAMPLE

Assume that some sort of experiment was conducted, where values of discharge
coefficients we calculated for a Venturi Meter, and we want to find the average discharge
coefficient for each run. The data is,
Run 1 Run 2 Run 3
0.93 0.86 0.96
0.89 0.84 0.94
0.91 0.91 0.93
We could set up the structure by,
exp.run.cd(1)=0.93;
exp.run.cd(2)=0.89;
exp.run.cd(3)=0.91;
exp.run(2).cd(1)=0.86;
exp.run(2).cd(2)=0.84;
exp.run(2).cd(3)=0.91;
exp.run(3).cd(1)=0.96;
exp.run(2).cd(2)=0.94;
exp.run(2).cd(3)=0.93;

50
 This sort of set up is not an efficient way of doing this. A better way would be,

cd=[0.93,0.89,0.91;0.86,0.84,0.91;0.96,0.94,0.93]
for r=1:3
for c=1:3
exp.run(r).cd(c)=cd(r,c);
end
end

This is easier to type, and lends itself to easy editing of the data.
Let us find the average for each run. We could say,

avgcd=zeros(1,3);
for r=1:3
avgcd(1,r)= mean([exp.run(r).cd])
end

This is equivalent to,

avgcd=zeros(1,3);
for r=1:3
avgcd(1,r)= mean([exp.run(r).cd(1), exp.run(r).cd(2),exp.run(r).cd(3)])
end

This then returns,

avgcd =

0.9100 0.8700 0.9433

4.3.7 Nested Structures

Since structures are very general data types, you can have fields that are of different
data types. This, we have already established. Just as we had nested cell arrays, we may
have structures within structures i.e., nested structures. The creation of nested structures
shall not be presented as a separate section. If you worked through this chapter, it is my
belief that you will already know how to set up nested structures.

51
52
Chapter 5

Logical Subscripting

This chapter will cover:

• Introduce the concept of logical subscripting.
• Using find to accomplish varied tasks.
• Using logical thinking to replace elseif ladders.

Brian D. Hahn says in his book that,

Those of us that grew up on more conventional programming languages
in the last century may find it difficult to think in terms of logical vectors
when solving general problems.
Coming from a PASCAL background, I do find it quite difficult to think in terms of
logical vectors. I found the chapter on logical vectors in Hahn’s book a little confusing.
Using the approach presented, I could not find a way of generalizing the approach.
Anyways, I later found out that the method presented by Hahn, was not always the best
choice. Then I came across the find command. It was just what I had been looking for!
It solved most of the problems that I had with logical vectors. In this chapter, I present
an approach that I view as being easier to understand. Of course, the choice of approach
is a matter of programming style.

5.1 Logical Vectors

I suppose you are wondering what logical vectors are. Quite simply, logical vectors are
vectors are vectors that consist of only zeroes and ones. Note that logical vectors are
distinctly different from numeric arrays that consist of zeroes and ones. To declare an
array as a logical vector, use the logical command e.g.,

a=logical([1 0 1 1 0]);

which declares a to be a logical vector.

Logical vectors can be used as subscript vectors. An example is,

53
p=logical([1 0 1 0]);
a=[1 2 3 4];

Now typing,

a(p)

returns the elements of a that correspond to the ones in p. Note that p and a do not
have to be the same dimension.
Relational operators can also be used in combination with logical vectors e.g.,

a>2

returns

ans =

0 0 1 1

So, typing,

a(a>2)

should return the elements of a that are greater than two. The result returned is,

ans =

3 4

The majority of the chapter in Hahn’s book presents various uses of this sort of
strategy.

5.2 Applications of Logical Vectors

Most of these examples are the same as those in Essential M ATLAB for Scientists and
Engineers, Third ed., Ch 5.
In each example that follows, Hahn’s method is presented first. My interpretation
shall then follow.

5.2.1 The find command

find returns the item numbers corresponding to a logical statement e.g.,

a=[1 2 3 4];
find(a>2)

returns

ans =

3 4

54
Now say, we had

a=[1 2 3;4 5 6];

find(a>2)

returns

ans =

2
4
5
6

Writing the matrix as,

a =

1 2 3
4 5 6

we see that find counts cumulatively in rows and then moves to the next column i.e.,
the item numbers for the matrix a are,

1 3 5
2 4 6

In the case of a one-dimensional array, item numbers like this are fine. When dealing
with matrices, we usually want row and column indices. This can be accomplished by,

[r,c]=find(a>2);

which returns the row index in variable r and column index in variable c i.e.,

r =

2
2
1
2

c =

1
2
3
3

The first element the find comes across that is greater than two is a(r(1),c(1)), and
so forth.
That’s all you need to know about find.

55
5.2.2 Discontinuous Graphs
Plot the following graph,
½
sin(x), sin(x) > 0
y(x) =
0, sin(x) ≤ 0

for x ∈ [0, 3π]

1. x=linspace(0,3*pi,100); %A slight change from the text
y=sin(x);
y=y.*y(y>0);
plot(x,y);
2. x=linspace(0,3*pi,100);
y=sin(x);
y(find(y<0))=0;
plot(x,y);
Note that the logical statement appears on the left hand side of the assignment state-
ment in my approach; and on the right hand side of the assignment statement in Hahn’s
approach. This shall be the case in many of the examples that follow.
Had the function been,
½
sin(x), sin(x) > 0
y(x) =
2, sin(x) ≤ 0

it would be rather difficult to use Hahn’s approach. Try it!

5.2.3 Avoiding Division by Zero

Plot the graph of sin(x)/x for x ∈ [−4π, 4π].
1. x=linspace(-4*pi,4*pi,100);
x=x+(x==0).*eps;
y=sin(x)./x;
plot(x,y);
2. x=linspace(-4*pi,4*pi,100);
x(find(x==0))=eps;
y=sin(x)./x;
plot(x,y);
Here eps has been used. eps is smallest absolute difference that can be represented by
M ATLAB. It is almost zero and is positive.

5.2.4 Avoiding Infinity

In this example we shall use y = tan(x). We know that this graph approaches infinity
as π/2 is approached. Thus, in plotting this graph, the large values should be filtered
out. All values with magnitude greater than 109 will be removed.
1. x=linspace(-3/2*pi,3/2*pi,100);
y=tan(x);
y=y.*(abs(y)<1e10);
plot(x,y);

56
 2. x=linspace(-3/2*pi,3/2*pi,100);
y=tan(x);
y(find(abs(y)>=1e10))=0;
plot(x,y);
Note that the logical statement was reversed in order to use my approach. Also, I
would prefer replacing,

y(find(abs(y)>=1e10))=0;

with

y(find(abs(y)>=1e10))=NaN;

NaN stands for not a number, and is useful in the cropping of plots. M ATLAB ignores
NaNs, and if plotting ordered pairs, ignores the element at the corresponding position
in the other vector. Surfaces may also be cropped using NaN.

5.2.5 Counting Random Numbers

Thus far, I have presented the use of logical vectors as filters. In this application, logical
vectors do not function as filters.
For this example, we shall generate a 1-by-7 random matrix and count the number
of elements that are less than 0.5.
1. r=rand(1,7);
sum(r<0.5);
2. r=rand(1,7);
size(find(r<0.5));
In this case, Hahn’s approach is more direct and easier to follow. My approach is a
little abstract and uses a dirty trick!

5.3 Replacing Stacked elseif Ladders

Consider the problem presented in Hahn’s book:
Income Tax
R10000 and less: 10% of the income.
Between R10000 and R20000: R1000 plus 20% of the amount by which
R10000 is exceeded.
More than R20000: R3000 plus 50% of the amount by which R20000 is ex-
ceeded.
Take incomes of R5000, R10000, R15000, R30000 and R50000.
Using an if construct, which is the usual approach,

inc=[5000 10000 15000 30000 50000];

tax=zeros(1,5);
for j=1:size(inc,2)
if inc(j)<=10000
tax(j)=0.10*inc(j);

57
elseif inc(j)>10000 && inc(j)<=20000
tax(j)=1000+0.20*(inc(j)-10000);
else
tax(j)=3000 + 0.50*(inc(j)-20000);
end;
end;

which returns,

tax =

500 1000 2000 8000 18000

Now, using Hahn’s logical vector approach,

inc=[5000 10000 15000 30000 50000];

tax=0.1.*inc.*(inc<=10000);
tax=tax+(inc>10000 & inc <= 20000).*(0.2.*(inc-10000)+1000);
tax=tax+(inc>20000).*(0.5.*(inc-20000)+3000);

What did not make sense at first was why tax was added to itself in the third and
fourth lines.
My approach is to use,

inc=[5000 10000 15000 30000 50000];

tax(find(inc<=10000))=0.1.*inc(find(inc<=10000));

tax(intersect(find(inc>10000),find(inc<=20000)))=1000 + ...
0.2.*(inc(intersect(find(inc>10000),find(inc<=20000)))-10000);

tax(find(inc>20000))=3000 + 0.5.*(inc(find(inc>20000))-20000);

Once again, a dirty trick is used. intersect finds the common elements of two
vectors. Can you figure out why this needs to be done. What would happen if you had
three logical statements joined by &?
After looking at this bit of code you are probably thinking that logical vectors are
not such a good idea. My advice is that you should use logical vectors where their use
is obvious to you. . . otherwise stick to elseif ladders!

58
Chapter 6

Function Handles, feval,

inline & Function
Overloading

This chapter covers:

• Uses of the function handle data type
• Uses of feval
• inline functions
• Iterative procedure requirements
• Symbolic variables

The chapter introduces function handles. Some individuals consider function han-
dles as being something quite abstract, when in fact they are quite simple and find use
in many areas of M ATLAB programming.
If you have a copy of Essential M ATLAB for Scientists and Engineers, you will notice
that the solution of differential a system of differential equations over a finite time inter-
val is accomplished by the use of an m-file that contains the definition of the system of
DEs. Also, in the section that provides the coding for Newton’s Method, two function
m-files are used viz., one for the function definition and one for the definition of the
derivative of the function. The use of inline functions reduces the number of m-files
required in Newton’s Method, as both the definition of the function and its derivative
can appear in the same m-file. For a system of FODEs1 , a single m-file can be used,
where an inline function represents the FODE system, and also contains the commands
that actually initialize solution of the system.
1 There are cases where it is easier to give function definitions in seperate m-files such as in

solving a system of FODEs with a constraint equation. See the section on First Order Differential
Equations and Runge-Kutta Methods for information.

59
 The question does then arise — how does one choose whether to use m-files or
inline functions to solve systems of DEs? As it turns out, one cannot use inline func-
tions when the system has non-constant variable coefficients. In simpler terms, one
may only use inline functions if the system has explicitly specified numeric coefficients (i.e.,
numbers must appear as coefficients) for the variables.

6.1 The inline Function

You may consider the inline function as being a M ATLAB function that allows the
creation of standard mathematical functions, an example of which is f (x) = x2 . We
know that we may substitute a valid value for x and get out an answer for f (x). This
is exactly the functionality that inline inparts to M ATLAB i.e. the ability to substitute
(valid) values into an expression. Enough said. Here’s an example. . .

⊲ E XAMPLE
Generate an inline function f (x, y) = x2 + y 2 and find f (x, y) for x = [1, 2, 3, 4] and
y = [5, 6, 7, 8].

Of course the task can be solved very easily using the following code:

x=[1 2 3 4]; y=[5 6 7 8];

f=x.ˆ2+y.ˆ2;

Okay, I know that you’re probably thinking that if this problem can be solved so
easily without inline functions, why should we use inline functions? Let’s say that the
problem was modified to read,
Given f (x, y) = x2 +y 2 . Find f (x, y) for x1 = [1, 2, 3, 4] and y1 = [5, 6, 7, 8];
and x2 = [10, 12, 13, 14] and y2 = [15, 26, 57, 78].

Possible approaches to solve the modified problem are:

1. Use the code above and to run the program, and then manually change the values
of x and y in the program and rerun the program.
2. Rewrite the above program so that it could accept input via the command win-
dow i.e.

x=input(’x ? ’);
y=input(’y ? ’);
disp([’xˆ2+yˆ2= ’, ’ ’, num2str(x.ˆ2+y.ˆ2)]);

The first approach is not the cleverest as if one had 100 sets of data, then it would
take a million years to generate the required output. The second method is good, as it
can be used for any number of data sets —it is probably the method that you thought
of. The problem with this is that you have to create an m-file. Another way of solving
the problem is to use the inline function. We could simply say, (for the first set of
data)

60
f=inline(’x.ˆ2+y.ˆ2’,’x’,’y’);
x=[1 2 3 4];
y=[5 6 7 8];
f(x,y);

Writing f(x,y) imparts a “mathematical sense” to the solution of the problem.

One may consider that inline generates an m-file similar to the second approach to
the solution of the problem.

TOOL 6.1 Let f (x) be any mathematical expression where x is a vector of variables, and
x = [x1 , x2 , . . . , xn ]. A substitutable mathematical function f (x) may be created by the
use of, f (x1 , x2 , . . . , xn )=inline(f (x, y), x1 , x2 , . . . , xn )

6.2 The Requirements of Iterative Procedures and

Symbolic Variables
Iterative procedures require some substitutable mathematical function(SMFs). This re-
quirement also holds for the use of the numerical differential equation solvers in M AT-
LAB . In M ATLAB SMFs can be created using m-files or the inline function. I know that
some individuals are probably thinking that I have forgotten that mathematical func-
tions may be created using symbolic variables in M ATLAB. I have not! As an example,
the solution the problem in the previous section may also be written as (this is especially
for those students that had a problem evaluating a symbolic expression in M ATLAB!),

syms x y %alternatively x=sym(’x’);y=sym(’y’);

f=x.ˆ2+y.ˆ2;
x=[1 2 3 4];y=[5 6 7 8];
eval(f);

I have found symbolic variables very useful. You can do many things with sym-
bolic expressions e.g. differentiate, integrate and plot symbolic expressions using
ezplot. Look up the relevant sections in M ATLAB help. The only downside of using
symbolic variables is that the Symbolic Math Toolbox is not part of the standard M AT-
LAB package. The University of Kwazulu Natal, for example, does not have the toolbox
installed on the computers in the student LANs. Also, sometimes lecturers will insist
that you not use the Symbolic Math Toolbox in your assignments. The most important
drawback of the Symbolic Math Toolbox is that it cannot be used to generate the SMFs
that the solvers in M ATLAB require. This is because feval cannot be used on symbolic
variables—only eval is valid. I guess that you’re now wondering what feval does.
The next section explains function handles and feval.

6.3 Function Handles and feval

It so happens that the term function handle is very appropriate. A function handle can
be thought of as being a link to an m-file or a link to a built-in M ATLAB function. feval
is then used to evaluate the function linked by the function handle. As I mentioned

61
that an in depth knowledge of function handles is not required at this stage, I will not
present much in this section. By the way, function handles are created with the @ oper-
ator.

⊲ E XAMPLE
Let’s say that you wanted to calculate the square root of a number. Let the number be
x. You could say sqrt(x). The more complex way of doing this is to say,

s=@sqrt;
feval(s,x);

If you know the function that you want to create a handle to, then you can just use the
name of the function prefixed with @, together with feval. An alternative to the code
above is,

feval(@sqrt,x);

Note that if sqrt required two input arguments, then you would have to say some-
thing like feval(s,x,y), where x and y are the two input arguments.
I know that it looks kind of stupid, but function handles are quite useful as ex-
plained next. You probably never will use them, but it is good to know.

⊲ C ASE IN U SE : O VERLOADED F UNCTIONS AND S NAPSHOTTING

Functions handles allow one to overload functions. Overloading just means providing
more than one definition of. Let’s take an example. Say your friend told you to travel to
the museum today on your own. Let’s assume that your drivers’ test is tomorrow. Let
us also say that the only available mode of transport was the bus. So, today you would
have taken the bus to the museum. Now, let’s say that you were successful in passing
your drivers’ test and were allowed to use the family car. If your friend had now told
you to travel to the museum on your own, you would use the family car (I think!).
Let’s look at the situation. On the two separate days, you were performing the same
task, viz., travelling to the museum on your own. What did change over the two days
was your definition of travelling. How you accomplished the task of travelling to the
museum was determined by your definition of travelling at that point in your life. In
effect by redefining travelling, you now have two definitions of travelling.
From a M ATLAB point of view, what is being said is that if you have an m-file and
for some reason decide to redefine the m-file by creating a new m-file with the same name in
a different directory, you can specify that M ATLAB should evaluate the m-file that existed
prior to supplying the redefinition. This is what I term snapshotting, as M ATLAB takes
a snapshot of the state of the current session of M ATLAB. This sort of approach may be
necessary if a program edits an existing m-file that is used later on in the program.

⊲ E XAMPLE
In the current directory, create isodd.m as follows,

%ISODD(x) tests if x is odd and returns an appropriate

%message.

62
function msg=isodd(x)
if mod(x,2)˜=0
msg=’odd’;
else
msg=’even’;
end;

• Now type f=@isodd in the command window.

• Create a directory within the current directory called fhandles. Change the
current directory to fhandles.
• Now, edit isodd.m as follows:

%ISODD(x) tests if x is odd and returns an appropriate

%message.

function msg=isodd(x)
if mod(x,2)˜=0
msg=[num2str(x), ’ is odd’ ];
else
msg=[num2str(x), ’is even’];
end;

Lastly,
• type isodd(3) in the command window. Take careful note of the result.
• type feval(f,3) in the command window. Is this result what you expect? Can
you explain what M ATLAB has done?

6.4 Making Subfunctions Global

An important use of function handles is to make subfunctions pseudo-global2 . A sub-
function is a function that appears below the main function of a function m-file. Let’s
say that you had the following hypothetical m-file (assume that x and y are integers),

function sumnums = nsum(x,y)

sumnums=x+y;

%the subfunction follows...

function meannums = nmean(x,y)

meannums=(x+y)/2;

nmean can only be called by nsum according to rules of scope. To make nmean
accessible from outside the function m-file, we edit it as follows,
2 This term is used as using function handles in this way does not make the subfunction global—

the subfunction is merely accessible beyond its usual scope. In other words, function handles allow
subfunctions to supersede their usual scope.

63
function sumnums = nsum(x,y)
sumnums=x+y;
nmean=@nmean;

%the subfunction follows...

function meannums = nmean(x,y)

meannums=(x+y)/2;

nmean can now be called from anywhere within M ATLAB, provided that the function han-
dle is passed as a parameter to the function or script file in which you want to make use of the
subfunction or declared as global3 , while the subfunction is in focus. We have effectively pro-
moted a subfunction to pseudo-global scope. This technique is useful if you have many
function m-files for a particular task. They can be condensed into fewer m-files. Let’s
say that we wanted to know if the mean of x and y is even or not. Now assume we have
another m-file that checks if the mean is even or odd, and is as follows:

function ismeaneven(mhandle,x,y)

if mod(feval(mhandle,x,y),2)==0
disp(’mean is even’)
else
disp(’mean is odd’);
end

Another advantage of using function handles in this way should be apparent—

inadvertently evaluating an incorrect overloaded function m-file is avoided4 .

3 The function handle will need to be declared as global from within the function m-file and in

the base workspace. evalin could be used here.

4 See the section on Overloaded Functions for an explanation as to why this is true.

64
 Part III

Creating GUIs

65
 This part covers the creation of GUIs in M ATLAB. Two examples (one easy and
the other slightly more involved) show how you can quickly create GUIs in M ATLAB.
Makes for good reading, for those that have always wanted to know how to create GUIs
in M ATLAB.

67
68
Chapter 7

Creating GUIs in M ATLAB

This chapter will cover:

• Common M ATLAB GUI components.
• A simple example on the creation of GUIs.
• Useful properties of GUI components.
• Accessing workspace variables from within a GUI.
• A more complex example on the creation of GUIs.

I remember the first time that I wanted to use GUIs in M ATLAB . . . it was when I
had to the Oil and Mineral Processing assignment in 2003. The calculations had been
completed in M ATLAB, and the output in the command window was acceptable, but
my friend FJ Nadaraju suggested that we make our solutions to the assignment more
interesting by the using GUIs. So off we went to Venter LAN, where I spent five hours
trying to figure out how to create a GUI in M ATLAB. I knew that GUIDE existed, and
initialised GUIDE. The first thing that I tried to do was to create a command button1
with text on it and then change the text when the user clicked on another command
button that was on the GUI. What a task it proved to be. I eventually managed to do
it, but then I realised that it would probably take me the entire semester to create a GUI
for the assignment. I asked lots of students and a couple of the lecturers to explain how
to create GUIs. Some did not know how to explain, and some simply said that it is
not really necessary to know how to create GUIs. The problem at that time was that I
did not understand function handles and structures. Anyway, we eventually handed our
assignments in without the GUI. Not even the final year students who are excellent at
programming in M ATLAB, knew how to create a GUI application in M ATLAB. Is there
hope for us students that know just a little2 M ATLAB programming? It is my hope that
this chapter helps students to create GUIs in M ATLAB.
1 I will sometimes refer to the M ATLAB uicontrols by their Visual Basic equivalents. The corre-

spondences between VB controls and M ATLAB uicontrols are presented later in this chapter.
2 Even though I have written this book, I still consider myself as being a beginner to program-

ming in M ATLAB. There are billions of things that one can do in M ATLAB, and I don’t think that
anyone can actually become a M ATLAB M ASTER!

69
7.1 M ATLAB GUI Components
M ATLAB has all the usual components that are required for the creation of GUIs. The
axes component is unique to M ATLAB 3 . In M ATLAB, GUI components are referred to
as uicontrols, which is short for user-interface controls.
M ATLAB has the following uicontrols4 :
Push Button The M ATLAB equivalent of a command button in VB.
Toggle Button A push button that looks as if it is pressed when it is selected, and looks
like a normal push button when it is not selected.
Radio Button Note that these are mutually exclusive i.e., only one from a group can be
selected at any one time.
Check Box Not mutually exclusive.
Static Text The equivalent of a label.
Edit Text equivalent of a textbox. Note that scrollable textboxes are not available in
M ATLAB.
Slider
Frame Note that axes cannot be framed with other uicontrols.
List Box
Pop Up Menu
Axes Used to plot data, or to load images
I assume that you have some idea of the uses and the common uses of the above uicon-
trols. If not, then ask someone or take a look at a text book on Visual Basic.

7.2 Getting Started

GUIDE, the GUI layout tool for M ATLAB, can be initialized by typing guide in the
command window. You can then drag and drop controls as you would do in Visual
Basic. You then need to save the file. The GUI consists of two files viz., and m-file that
contains the callbacks(which are the equivalents of Subs in VB or Procedures in Pascal),
and the .fig file that contains the layout information for the GUI. Note that both files
have the same name. Also, if you change the tag(same as Name in VB) property for
a component in the GUI, then you cannot just save the file. You need to replace the
existing file i.e., choose “save as” and type in the same file name and confirm that the
existing file should be replaced.
As previously mentioned, structures are used widely in GUIs. It is not necessary
to have an in-depth knowledge of structures. You may think of a structure as being a
file-folder and the fields of the structure as being the stuff in the file-folder. In simple
terms, a structure is the thing before the dot and the field is the thing after the dot. As
an example, in VB we would say something like

text1.text="Someting"
3 When programming in Visual Basic, one can insert an Excel Chart Control into a GUI. I sup-

pose that this could be regarded as the VB equivalent of the axes component of M ATLAB.
4 A description will follow where necessary. If a description does not follow, then assume that

the operation of the uicontrol is the same as the corresponding control in Visual Basic.

70
to change the text in text box. In this case, the structure is text1 and the field is text.
In M ATLAB, however, the automatically generated structure handles contains links
to all information associated with the GUI e.g., the text that appears in a textbox. It
is therefore, what I like to call a global structure. Note that although I said global, the
structure and its fields are restricted to the GUI i.e., the fields cannot be accessed from
outside the GUI. It therefore makes sense to store variables that we want to keep global
(with respect to the GUI) in the handles structure. By doing this, we are able to share
information between callbacks. Note that properties (e.g., the text in a textbox) is set
using set, and the current state of the property can be retrieved using get. The usage
of these functions shall become apparent via the examples in this chapter.

7.3 An Easy Example: Sum of Two Numbers

⊲ E XAMPLE
Design a GUI in M ATLAB, that accepts two numbers and finds the sum of the numbers.
The input must be via edit text uicontrols and the output should appear in an edit text
uicontrol. A push button should be used to initialise the calculation.
Solution The GUI was constructed as follows,

Note that
1. The properties of each UIC is edited via the Property Editor (it’s in the Tools
menu). Note that the properties of the currently selected UIC is displayed in
the property editor.
2. The edit text UIC for the first number has its tag value set to txtnum1, in keeping
with what one would do in Visual Basic. The edit text UIC for the 2nd number
has a tag of txtnum2
3. The edit text UIC for the output of the sum has a tag value of txtsum
4. The push button has a tag value of cmdsum.
5. All edit text UICs have the string value set to blank i.e., nothing

71
 In the examples that follow, I will not explicitly state what the property values for
each uicontrol on the GUI are. These shall be quite obvious. The corresponding m-file
for the callbacks is as follows,

72
function varargout = eg1(varargin)
% EG1 M-file for eg1.fig
% EG1, by itself, creates a new EG1 or raises the existing
% singleton*.
%
% H = EG1 returns the handle to a new EG1 or the handle to
% the existing singleton*.
%
% EG1(’CALLBACK’,hObject,eventData,handles,...) calls the local
% function named CALLBACK in EG1.M with the given input arguments.
%
% EG1(’Property’,’Value’,...) creates a new EG1 or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before eg1_OpeningFunction gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to eg1_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE’s Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help eg1

% Last Modified by GUIDE v2.5 30-Aug-2003 09:42:22

% Begin initialization code - DO NOT EDIT

gui_Singleton = 1;
gui_State = struct(’gui_Name’, mfilename, ...
’gui_Singleton’, gui_Singleton, ...
’gui_OpeningFcn’, @eg1_OpeningFcn, ...
’gui_OutputFcn’, @eg1_OutputFcn, ...
’gui_LayoutFcn’, [] , ...
’gui_Callback’, []);
if nargin & isstr(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before eg1 is made visible.

function eg1_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure

73
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to eg1 (see VARARGIN)

% Choose default command line output for eg1

handles.output = hObject;

% Update handles structure

guidata(hObject, handles);

% UIWAIT makes eg1 wait for user response (see UIRESUME)

% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = eg1_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure

varargout{1} = handles.output;

% --- Executes during object creation, after setting all properties.

function txtnum1_CreateFcn(hObject, eventdata, handles)
% hObject handle to txtnum1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.
if ispc
set(hObject,’BackgroundColor’,’white’);
else
set(hObject,’BackgroundColor’,get(0,’defaultUicontrolBackgroundColor’));
end

function txtnum1_Callback(hObject, eventdata, handles)

% hObject handle to txtnum1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,’String’) returns contents of txtnum1 as text

% str2double(get(hObject,’String’)) returns contents of txtnum1 as a double

74
% --- Executes during object creation, after setting all properties.
function txtNum2_CreateFcn(hObject, eventdata, handles)
% hObject handle to txtNum2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.
if ispc
set(hObject,’BackgroundColor’,’white’);
else
set(hObject,’BackgroundColor’,get(0,’defaultUicontrolBackgroundColor’));
end

% --- Executes on button press in cmdsum.

function cmdsum_Callback(hObject, eventdata, handles)
% hObject handle to cmdsum (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

%%%%%% THE LINES THAT ARE DELIMITED IS WHAT HAS BEEN ENTERED IN THE M-FILE%%%%%%%%%%
x=get(handles.txtnum1,’string’);
x=str2double(x);
y=get(handles.txtnum2,’string’);
y=str2double(y);
s=x+y;
set(handles.txtsum,’string’,s);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% --- Executes during object creation, after setting all properties.

function txtsum_CreateFcn(hObject, eventdata, handles)
% hObject handle to txtsum (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.
if ispc
set(hObject,’BackgroundColor’,’white’);
else
set(hObject,’BackgroundColor’,get(0,’defaultUicontrolBackgroundColor’));
end

function txtsum_Callback(hObject, eventdata, handles)

% hObject handle to txtsum (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB

75
% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,’String’) returns contents of txtsum as text

% str2double(get(hObject,’String’)) returns contents of txtsum as a double

% --- Executes during object creation, after setting all properties.

function txtnum2_CreateFcn(hObject, eventdata, handles)
% hObject handle to txtnum2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.
if ispc
set(hObject,’BackgroundColor’,’white’);
else
set(hObject,’BackgroundColor’,get(0,’defaultUicontrolBackgroundColor’));
end

function txtnum2_Callback(hObject, eventdata, handles)

% hObject handle to txtnum2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,’String’) returns contents of txtnum2 as text

% str2double(get(hObject,’String’)) returns contents of txtnum2 as a double

Everything not delimited has been automatically generated by GUIDE. As you can
see, everything is accessed and set using the handles structure {handles. I suppose that
this example is too simple. . . so on to something more interesting. By the way can you
make the push button disappear when it is clicked? All you have to do is say,

set(handles.cmdsum,’visible’,’off’);

In essence, to retrieve the current state of a property of a component5 , all one needs
to say is,

get(handles.uicontrol,’property_to_retrieve’);

and to set a property value,

set(handles.uicontrol,’property’,’value_property_must_be_set_to’);
Quite simple isn’t it? After looking at this, you’re probably thinking that GUI creation
in M ATLAB is quite similar to GUI creation in VB. If you’re not thinking this, then GUI
creation has not clicked as yet. With a little practice(say, 2 hours), and a couple of
tricks(which I shall mention), GUI creation should become second nature. By the way,
the hints in the autogenerated m-file are very helpful.
5 Of course, the retrieved property value can be assigned to a variable for later use.

76
7.4 UIControls: Useful Properties and Related Stuff
I have already mentioned the UICs that are available in M ATLAB. Most people will only
need to know about a few properties for each UIC. The aim of this section is to mention
the properties that are most commonly used.

7.4.1 Push Button

The most commonly used property is string. It can be retrieved and set by the get
and set functions respectively. You most probably will not find it necessary to change
this property during run time, but it is good to know.

7.4.2 Radio Button

To retrieve the current state i.e., selected or not selected do,

get(handles.radio_button_tag,’value’)

which returns either a zero(not selected) or one. A similar thing can be done for check
boxes.

7.4.3 Sliders
The maximum and minimum value of the slider can be set by saying,

set(handles.slider_tag, ’max’, maxvalue);

set(handles.slider_tag, ’min’, minvalue);

Since sliders can be either vertical or horizontal, the increment for both these cases
can be set. When the triangles on the slider control is clicked, it increments or decre-
ments by a fraction of the range defined by the maximum and minimum slider values.
This increment or decrement is known as the slider step As an example,

set(handles.slider_tag, ’max’, 10);

set(handles.slider_tag, ’min’, 0);
set(handles.slider_tag, ’sliderstep’,[0.1 0.2]);

will cause the slider value to increment or decrement by 0.1 × 10 = 1 each time the
arrow on the slider is clicked, if the slider is horizontal. A increment/decrement of 2 for
a vertical slider is defined by the above commands. It is useful to have a edit text UIC
linked to the slider that indicates the current slider value.

7.4.4 List Box

List boxes are used to display a set of results, or to retrieve options selected by the user.
In order to allow the user to select more than one option simultaneously, the max and
min values need to be set such that the maximum minus the minimum value is greater
than one. This can be done by,

set(handles.list_box_tag,’max’,2);
set(handles.list_box_tag,’min’,0);

77
 The issuing of the above commands also defines that the user may select 0, 1 or 2
of the options in the list box. This is fairly obvious as we are setting the maximum and
minimum values.
Data may be output to the list box UIC by,

set(handles.list_box_tag,’string’,output_string,’value’,1)

Note that ’value’,1 must always appear when you want to output to a list box.
This clears the list box and populates the list box with the string that you wish to output.
The items indices of the items currently selected can be retrieved during the list box call
back by,

get(gcbo,’value’);

which is equivalent to,

get(handles.list_box_tag,’value’)

gcbo means get current callback object, which from within the callback for the list box
refers to the list box! A trick to get the actual item that has been selected is to get the
items as a cell array and to then reference the correct position in the cell array6 i.e.,

listitems=get(handles.list_box_tag,’string’) %returns a cell array

indices=get(handles.list_box_handle,’value’);
itemselected=listitems{indices}; %note the curly braces

7.4.5 Pop Up Menu

Pop up menus should actually be called drop down menus. It makes sense if you look
at the appearance of the UIC. Pop up menus are used to list options that can be selected
by the user. It is quite useful if you want to make a GUI more compact, in the case of
having many user options that would usually be listed in a list box. Pop up menus can
be regarded as a fancy list box, thus their operations are equivalent.

7.4.6 Axes
If you have more than one axes7 on a GUI, you need to specify the axes on which you
wish to plot data. This can be done by,

Axes(handles.axes_tag)

which sets the active axes.

6 Note that to create and reference cell arrays, curly braces are used. These braces distinguish

them from numeric and string arrays. Just as a review, recall that cell arrays allow one to store
mixed data types in a single structure i.e., you can store numerics and strings together without any
problem. Cell arrays can also be used to display multiple lines of text in edit text boxes, and to
display multiple items in pop up menus(drop down menus). Note that the cell arrays must be row
vectors i.e., set up in columns to force multiple lines.
7 This is not a grammatical error. A single UIC of this type is referred to as axes in M ATLAB .

78
7.4.7 Accessing Workspace Variables
You will at some point want to access variables that are in the main (base) M ATLAB
workspace i.e., variables that are visible from the command window scope. You can do
this by using evalin8 . Let’s say that we have variables x and y in the workspace and
want to plot these on an axes(with tag set to axes1) on a GUI. This can be accomplished
by,

GUIx=evalin(’base’,’x’);
GUIy=evalin(’base’,’y’);
axes(handles.axes1); %make sure its the right axes
plot(GUIx,GUIy);

This function is very useful as you can do calculations based on variables that are
in the base workspace from within the scope of the GUI and push the results into the
main workspace so that you can access these results later from the command window.
For example, say we wanted to calculate z where z = x2 + y 2 and store the result in
the main workspace. It is assumed that x and y are defined in the main workspace. We
would need to say,

evalin(’base’,’z=x.ˆ2+y.ˆ2’);

Now, if you have variables within the GUI and want to make these accessible from
the workspace, then you will have to,
1. declare them as global9 from within the GUI.
2. declare them as global from the command window.
Say we had a variable, x, only visible from inside the GUI. To accomplish making x
visible from the command window, we could say,

global x;
evalin(’base’,’global x’);

This stores x in the command window workspace.

7.5 A More Complex Example

This example makes use of axes, edit text boxes, static text, and list box.
I still consider this example to be quite basic, but it is informative. If you understand
this example, then more complex GUIs should not pose a problem.

⊲ E XAMPLE
Design a GUI that,
8 Note that evalin cannot have nested arguments. Also, the expression must be defined in the

scope defined by the first argument i.e., if the expression contains x, y and z and the first argument
is ’base’, then these variables need to be already defined in the base workspace.
9 If you are within a function or GUI and want to make the variable visible from the command

line, then you need to declare the variable as global both in the function or GUI and in the com-
mand window. Note that once a variable is global its value can be changed from within any scope.
This can lead to unexpected errors. If you must use global variables, be very careful.

79
 • accepts a vector of x values and a vector of y values
• allows the user to choose the order of polynomial fit via a list box
• plots the input data and shows the fit graphically on an axes UIC on the GUI
Solution
The GUI layout and sample run is as follows,

It should be fairly obvious from the m-file for the GUI, what property values have been
set. The m-file is,

function varargout = eg2(varargin)

% EG2 M-file for eg2.fig
% EG2, by itself, creates a new EG2 or raises the existing
% singleton*.
%
% H = EG2 returns the handle to a new EG2 or the handle to
% the existing singleton*.
%
% EG2(’CALLBACK’,hObject,eventData,handles,...) calls the local
% function named CALLBACK in EG2.M with the given input arguments.
%
% EG2(’Property’,’Value’,...) creates a new EG2 or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before eg2_OpeningFunction gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to eg2_OpeningFcn via varargin.

80
%
% *See GUI Options on GUIDE’s Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help eg2

% Last Modified by GUIDE v2.5 04-Jan-2005 12:54:41

% Begin initialization code - DO NOT EDIT

gui_Singleton = 1;
gui_State = struct(’gui_Name’, mfilename, ...
’gui_Singleton’, gui_Singleton, ...
’gui_OpeningFcn’, @eg2_OpeningFcn, ...
’gui_OutputFcn’, @eg2_OutputFcn, ...
’gui_LayoutFcn’, [] , ...
’gui_Callback’, []);
if nargin & isstr(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before eg2 is made visible.

function eg2_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to eg2 (see VARARGIN)

% Choose default command line output for eg2

handles.output = hObject;

% Update handles structure

guidata(hObject, handles);

% UIWAIT makes eg2 wait for user response (see UIRESUME)

% uiwait(handles.figure1);

%%%%%%%%%%%%Set the items that are displayed in the popup menu%%%%%%%%%%

set(handles.puorder,’string’,{’1’;’2’;’3’;’4’})
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

81
% --- Outputs from this function are returned to the command line.
function varargout = eg2_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure

varargout{1} = handles.output;

% --- Executes during object creation, after setting all properties.

function txtx_CreateFcn(hObject, eventdata, handles)
% hObject handle to txtx (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.
if ispc
set(hObject,’BackgroundColor’,’white’);
else
set(hObject,’BackgroundColor’,get(0,’defaultUicontrolBackgroundColor’));
end

function txtx_Callback(hObject, eventdata, handles)

% hObject handle to txtx (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,’String’) returns contents of txtx as text

% str2double(get(hObject,’String’)) returns contents of txtx as a double

% --- Executes during object creation, after setting all properties.

function txty_CreateFcn(hObject, eventdata, handles)
% hObject handle to txty (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.
if ispc
set(hObject,’BackgroundColor’,’white’);
else
set(hObject,’BackgroundColor’,get(0,’defaultUicontrolBackgroundColor’));
end

82
function txty_Callback(hObject, eventdata, handles)
% hObject handle to txty (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,’String’) returns contents of txty as text

% str2double(get(hObject,’String’)) returns contents of txty as a double

% --- Executes during object creation, after setting all properties.

function puorder_CreateFcn(hObject, eventdata, handles)
% hObject handle to puorder (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.

% See ISPC and COMPUTER.
if ispc
set(hObject,’BackgroundColor’,’white’);
else
set(hObject,’BackgroundColor’,get(0,’defaultUicontrolBackgroundColor’));
end

% --- Executes on selection change in puorder.

function puorder_Callback(hObject, eventdata, handles)
% hObject handle to puorder (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,’String’) returns puorder contents as cell array

% contents{get(hObject,’Value’)} returns selected item from puorder

% --- Executes on button press in pbfitdata.

function pbfitdata_Callback(hObject, eventdata, handles)
% hObject handle to pbfitdata (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

%%%%%%%%%%%%%%%% Fitting the Data %%%%%%%%%%%%%%%%%%

handles.fitdata.x=str2num(get(handles.txtx,’string’));
handles.fitdata.y=str2num(get(handles.txty,’string’));
itemlist=get(handles.puorder,’string’);
indices=get(handles.puorder,’value’);
itemsel=itemlist{indices};

83
n=str2double(itemsel);
p=polyfit(handles.fitdata.x,handles.fitdata.y,n);
handles.fitdata.xi=linspace(min(handles.fitdata.x), ...
max(handles.fitdata.x),100);
handles.fitdata.yi=polyval(p,handles.fitdata.xi);
cla;
hold on;
axes(handles.axes);
plot(handles.fitdata.xi,handles.fitdata.yi,’g-’);
plot(handles.fitdata.x,handles.fitdata.y,’ro’);
hold off;
legend(’fitted data’,’original data’,4);
fresult=cell(4,1);
for j=1:n
fresult{j,1}=num2str(p(j));
end
set(handles.lstp,’string’,fresult,’value’,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% --- Executes during object creation, after setting all properties.

function txtp_CreateFcn(hObject, eventdata, handles)
% hObject handle to txtp (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.
if ispc
set(hObject,’BackgroundColor’,’white’);
else
set(hObject,’BackgroundColor’,get(0,’defaultUicontrolBackgroundColor’));
end

function txtp_Callback(hObject, eventdata, handles)

% hObject handle to txtp (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,’String’) returns contents of txtp as text

% str2double(get(hObject,’String’)) returns contents of txtp as a double

% --- Executes during object creation, after setting all properties.

function lstp_CreateFcn(hObject, eventdata, handles)
% hObject handle to lstp (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

84
% Hint: listbox controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,’BackgroundColor’,’white’);
else
set(hObject,’BackgroundColor’,get(0,’defaultUicontrolBackgroundColor’));
end

% --- Executes on selection change in lstp.

function lstp_Callback(hObject, eventdata, handles)
% hObject handle to lstp (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,’String’) returns lstp contents as cell array

% contents{get(hObject,’Value’)} returns selected item from lstp

Comments on Code
The following should be noted:
1. The items in the popup menu have been created in the opening function. This
function is called just before the GUI is made visible. See comments in the m-file.
This could have been done a bit differently—you could have used the property
editor to set the string of the popup menu.
2. All variables are stored in the handles structure. This is not required in this
case, as the variables are only used within a single function. However, when data
needs to be made available to all functions in the scope of the GUI, then the data
must be stored in the handles structure. You could have also more simply said
something like,
handles.xi=str2num(get(handles.txtx,’string’));
instead of the more elaborate structure presented in the example.
3. By using num2str, the user can input the x and y vectors as comma-separated
lists of space-delimited lists. If you use str2double the code would not work as
expected. See M ATLAB help on these two functions, and as to why they are used
here.
4. cla clears the current axes.

That’s it for GUIs. . .

85
86
 Part IV

Numerical Methods

87
 This part of introduces some basic numerical methods and the use of M ATLAB as
a tool in numeric calculations. Topics such as integration, differentiation, solution of
FODEs are covered. ODE solvers are covered in detail.10

10 Students undertaking the course Applied Computer Methods at The University of Kwazulu

Natal, SA will find this section particularly useful.

89
90
Chapter 8

Basic Integration and

Differentiation, and Data
Regression

This chapter will present

• some basic uses of M ATLAB in the numeric solution of problems. Most of the stuff
should have been covered in a Computer Fundamentals course. Some familiarity
with Mathematics I is required. However, if you have done additional maths in
secondary school then you should be OK.
• an introduction to loop vectorization
• rewriting program segments so that the code takes advantage of M ATLAB fea-
tures
• “matrix thinking”
• use of the polyfit and polyval and interp1 functions in data regression

8.1 Why Numerical Methods?

From the Mathematics covered in first year, we know that it is not always possible to
find an analytic solution to a problem e.g., one cannot evaluate the analytic solution of
Z
2
ex dx

by elementary methods of integration. However, we can evaluate the above integral

graphically by estimating the area under the graph of the integrand. We can, of course,
also find the derivative of a function by approximating the tangent to the graph of the
function.

91
8.2 Basic Numerical Integration
The integral of a function can be evaluated if we have points over a specified do-
main and the corresponding values in the range of the function. Say we have x =
[1, 2, 3, 4] and y = [2, 4, 6, 8]. It is clear in this case that the relationship between x and y
is y = 2x. Let us evaluate the integral of this function analytically and then numerically.
We shall at the end compare the two results. Although this example is quite trivial, it is
important that you understand the process involved in evaluating the numerical solu-
tion, and also how we link up the algorithm for the solution to its equivalent M ATLAB
coding, as we proceed. Before we proceed, we need to establish how to estimate the
integral of a function given pairs of points 1

TOOL 8.1 Let x and y be abscissa and ordinate vectors,respectively, each with length n. The
integral of the function y = f (x) may be approximated by the Trapezoidal Rule,
n−1
X (yi+1 − yi )
(xi+1 − xi ) (8.1)
2
i=1
or,
n
X (yi − yi−1 )
(xi − xi−1 ) (8.2)
2
i=2

Analytic Solution
• Clearly, x varies from 1 to 4. In M ATLAB, to access the first and last elements
of a 1 dimensional array, we could say x(1) and x(end), respectively. Note
that saying x(end) is interpreted as x(end,1) i.e. the last element in the first
dimension (rows). This gives us the limits of integration.
• We evaluate,
Z4
2x dx
1

which gives us 15.

Numeric Solution The equivalent M ATLAB code can be represented by

x=[1 2 3 4];
y=[2 4 6 8];
npoints=size(x,2) %number of points

A=0; %holds the value of the area under the curve

for i=1:npoints-1
deltax=x(i+1)-x(i) %delta x
1 It is possible to calculate integrals and derivatives using left-hand right-hand and central

differences. The central difference method requires an odd number of points and thus is restrictive.
For this reason it shall not be presented here.

92
yav=0.5.*(y(i+1)+y(i)) %average y value
dA=deltax.*yav;
A=A+dA;
end;

The numerical approximation gives 15.0 as an answer.

I think that I should mention at this point, that it is very easy to fall into the trap of
writing code such as that above, that does not harness the matrix capabilities of M AT-
LAB . Also, it is possible to vectorize loops. Vectorization of loops is considered an
artform by many as it is sometimes quite difficult to spot how one goes about vector-
izing a particular loop. Say, we wanted to find the sum of the squares of the first 50
counting numbers. We could write,

summ =0;
for num=1:50
sum =summ+numˆ2;
end;

A more elegant way of writing the same thing is using “matrix thinking”,

summ =[];
for num=1:50
summ =[summ numˆ2];
end;
sumsq=sum(summ)

Using loop vectorization, we get,

num=1:50
sumsq=sum(n.ˆ2)

There are obviously a few other ways of rewriting the above code. You are probably
wondering why you need to 3 different ways of writing the same bit of code. In my
opinion, the advantage(s) of each of the above versions is, respectively,
• if you have programmed in another language, then the first version is probably
the one that popped up in your mind. I think of this type of programming as be-
ing linear thinking. The advantage of this method is that the solution to a problem
is very apparent, albeit (usually) the version that will take the longest to run.
• matrix thinking is useful in iterative procedures. You shall see this later on. Also,
it sometimes simplifies the solution as M ATLAB has many functions for dealing
with matrices. It is also easier to visualize this solution, in some cases. Note that
this option can consume a large amount of memory.

93
 • if you want a program that runs as quickly as possible, then you will want to
vectorize your code. Note that sometimes it is impossible to vectorize code, and
that the increase in performance sometimes adds too much complexity actual
code. Beware, that for loops can only be vectorized if the order of evaluation of
the operations in the loop is irrelevant. Anyway, since the release of M ATLAB 6.5
the gain in performance by vectorization of code has become negligible, as the
looping structures have been optimized greatly.
Another popular method of quadrature is Simpson’s One Third Rule. It is a combination
of the Trapezoidal and Midpoint rules for quadrature. The integral over an interval is
defined as, (where a and b are the end-points of the interval)
N
P xi −xi−1 ¡ ¡ xi +xi−1 ¢¢
6
f (xi ) + f (xi−1 ) + 4f 2
i=2
(8.3)
N −1
P xi+1 −xi ¡ ¡ xi+1 +xi ¢¢
6
f (xi+1 ) + f (xi ) + 4f 2
i=1

The writing of the M ATLAB code shall be left as exercise for the reader.

8.3 Basic Numerical Differentiation

We base numerical differentiation on the definition of the derivative. Recall in Grade 2
you learned that the slope of a line was given by rise over run. Since this is one of the first
things we learned, we would expect the definition to behave erratically for non–linear
situations. This is indeed the case.
Once again, it is worthy to note that right–hand, left–hand and central differences
may be used. The former two shall be presented here.

TOOL 8.2 Given abscissa and ordinate vectors x and y both with n elements, we define the
derivative at xi using left–hand differences as,
yi − yi−1
f ′ (xi ) = (8.4)
xi − xi−1

and using right–hand differences as,

yi+1 − yi
f ′ (xi ) = (8.5)
xi+1 − xi

We are of course aware that the derivative is defined with regard to a particular
point. Sometimes we may wish to evaluate the average gradient over an interval. This
is achieved by extending the above definitions a little.

TOOL 8.3 Given abscissa and ordinate vectors x and y both with n elements, we define the
average derivative over x ∈ (x(1), x(n)) using left–hand differences as,
n
£ ¤ 1 X yi − yi−1
f ′ (x) = (8.6)
av n xi − xi−1
i=2

94
and using right–hand differences as,
n−1
£ ¤ 1 X yi+1 − yi
f ′ (x) = (8.7)
av n xi+1 − xi
i=1

It is important to note that for left–hand difference the summation starts at 2 and
ends at n. If the indexing started at 1, we would have x(0) and y(0), for i = 1. Unlike
a language like, say, Visual Basic, M ATLAB does not allow an Option Base 0 option.
For the right–hand difference, if the indexing ended at n, we would have x(n+1) and
y(n+1), which clearly do not exist, as x and y are of length n.
To end of this chapter, here is the M ATLAB code for finding the derivative vectors
using left–hand differences. The code for right–hand differences is similar and the onus
is on the you to write the code as an exercise.

%this code generates a vector that contains the derivative at

%each point of a vector x
%Assume that x and y are already defined i.e. have been input into Matlab
n=size(x,2);
dydx=[];
for i=2:n
dydx=[dydx (y(i)-y(i-1))/(x(i)-x(i-1))]
end;

Alternatively, using loop vectorization,

%this code generates a vector that contains the derivative at

%each point of a vector x
%Assume that x and y are already defined i.e. have been input into Matlab
n=size(x,2);
dydx=zeros(1,n); %this IS required-we need to "declare" variable
i=2:n
dydx(i)=(y(i)-y(i-1))/(x(i)-x(i-1))

There are just 2 other items that deserve mention in this chapter viz., fitting polyno-
mials to data and the use of the polyval function. The former will be presented as a
tool.

8.4 Data Regression

TOOL 8.4 Given x and y where the relationship between x and y can be represented by y =
axn + bxn−1 + cxn−1 + · · ·, we may generate a vector p where p = [a, b, c, · · ·] by using the
polyfit function i.e., polyfit(x,y,n)

Lastly, say we wanted to show graphically, the fit achieved by the use of the polyfit
function. We could do this by the following code (note the use of polyval)

95
%assuming that x and y are defined,and that we want a fit of order n

p=polyfit(x,y,n);
%generate a high resolution of x points within the domain of x values
xi=linspace(min(x,2),max(x,2),100);
yi=polyval(p,xi)
plot(x,y,’r’) %generate a red plot
hold(’on’); %do the next plot on the same figure
plot(x,y,’b’);

In the notes by Dr. Rawatlal, mention is made of the function interp1. interp1
performs one-dimensional interpolation and the general form of the function is,

interp1(known\_x,known\_y,’method’)

Note that the method is an optional argument. Linear interpolation is assumed if this
argument is omitted. I like to think of interp1 as being a special case of polyfit i.e.,
we are fitting a polynomial of order one to the given xy data. Thus, it is not a necessity
to use the interp function—it is worth your while to know of it, as it decreases the
number of lines that you have to type. Thus far, I think I have found use for using the
interp in the generation of a smooth curve that passes through a given set of data
points. An example would be the generation of a titration curve, where you do not
know the order of fit. The syntax I used in that case was,

%assume x and y are defined

xi=linspace(min(x,2),max(x,2),1000)
yi=interp1(x,y,xi,’splines’)
plot(x,y,xi,yi)

An example of using spline fitting will be presented in the exercises for this chapter.

8.5 When the Derivative Does Not Exist at a Point. . .

When the derivative fails to be defined at a point, the methods presented thus far also
fail. The teo simplest approaches to rectify this sort of situation is to use central differences
or the definition of the limit.

TOOL 8.5 To calculate an approximation of f ′ (xi ) where f ′ (xi ) does is not defined (but
f ′ (xi+1 ) and f ′ (xi−1 ) are defined), central differences may be applied i.e.,

f (xi+1 ) − f (xi−1 ) f (xi+1 ) − f (xi−1 )

f ′ (xi ) = =
xi+1 − xi−1 2∆x
It must be noted that in order to apply this technique, xi must not be a boundary point i.e.,
the technique is valid for interior points only. When the derivative over an interval exists at the
boundary points and its nearest neighbour, then left and right hand differences, together with
central differences may be used to approximate the derivative over the entire interval. Note also,
that ∆x represents a fixed step length.

96
 As you can see, the central difference technique is quite restrictive. If the central
difference technique of calculating the derivative at a point fails, then one may use the
definition of the limit. Since this requires choosing a successively smaller value for h, it
requires an iterative approach, where either the number of iterations is constrained or a
tolerance criterion is used. Due to the iterative nature of the process, this method shall
be discussed later.

97
98
Chapter 9

First Order Differential

Equations and Runge-Kutta
Methods

This chapter covers:

• Euler’s Method and the solution of First Order Differential Equations
• Transformation of mixed-order systems to FODE systems
• Runge-Kutta1 Methods in M ATLAB
• Solution of DEs in M ATLAB using ODE45

It is a certainty that you will encounter differential equations during your Engineer-
ing course. A famous example in Chemical Engineering is the stirred tank problem,
where one has to determine the mass fraction of salt in a continuously stirred tank as
a function of time, given an initial condition. In applied math, many systems are mod-
elled by differential equations e.g., the differential equation that models a simple har-
monic oscillator (ẍ(t) = ẋ(t)). Of course, one needs to be able to find solutions to such
systems. Sometimes, the analytic solutions can be exceedingly difficult to find, and in
some cases may not even exist. It is, then, imperative to develop numerical methods of
solving such problems.
Since this is not a book on numerical methods, the main focus of this chapter will
be on the use of M ATLAB in the finding numerical solutions to differential equations.

1 Named in honour of two German applied mathematicians, Runga-Kutta Methods refer to mul-

tistage, single-step methods for solving differential equations where y(tn+1 ) is based on a linear
combination of slopes calculated previously. Single-step methods only require y(tn−1 ) i.e. the so-
lution at the previous point to calculate y(tn ). Refer to a book on numerical methods if you want
to see how to use Runga-Kutta Methods to solve differential equations by hand.

99
9.1 Euler’s Method
Euler’s Method may be used to solve First Order Differential Equations(FODEs) where
the initial condition is known. Note that this can be regarded to be a very crude method
of solving FODEs.

dy
TOOL 9.1 Let dx
= f (x, y) where y(0) = y1 is known. The solution of the FODE may be
found by,

yi+1 = y1 + hf (xi , yi )

A few other points of note on Euler’s Method are:

1. x ∈ [a, b] i.e.x is defined over an interval with upper and lower bounds of b and a
respectively.
2. The step size is defined as h i.e., xi+1 − xi = h
3. The number of intervals(not the number of points!), n is defined as (b − a)/h
4. h must be chosen—there isn’t an method of selecting an optimal h. A tolerance
criterion can be applied to check if the solution obtained by Euler’s Method is
acceptable, if an (analytic) solution is available.
On the issue of number of points versus number of intervals, it is worthwhile to
think of it like this,
For 2 intervals, 3 points are required viz., the boundary points and the
partition point. In other words, for n intervals, there are 2 boundary points
and (n − 1) partition points yielding a total of 2 + n − 1 points i.e., n + 1
points.
Although, this may seem to be a trivial issue, it is important as one needs to keep
this in mind when using linspace and vector initialization methods. This shall be
seen in the example that follows. Based on this, one can either assume a value for h or a
value for n, as these 2 values are related by the definition of n. Anyway, if considering
that we defined n as being the number of intervals, it would not make sense if we chose
an h that resulted in a partial interval.

⊲ E XAMPLE
Given,
dy
= 4x2 , y(0) = 3
dx
Use Euler’s Method to find y(10). Show the solution graphically. Also find the analytic
solution for y(10) and compare the results obtained.

In order to apply Euler’s Method we require h. Let’s assume that n = 40. We also
know that the interval of interest is [0, 10]. The corresponding code is,

f=inline(’4.*x.ˆ2’,’x’,’y’);
nintervals=40;
a=0;
b=10;
h=(b-a)/nintervals;

100
y=zeros(1,nintervals+1);
x=a:h:b;
y0=3;
y(1)=y0;
for i=1:nintervals %perform over each interval
y(i+1)=y(i)+h.*f(x(i),y(i));
end

Comments on Code
1. Note that the y has n + 1 elements. All you have to do is remember that the
additional element in the array is the initial condition.
2. x=a:h:b could be replaced by x=linspace(a,b,n+1). It is a matter of pref-
erence. I prefer the latter method as the dimension of the arrays created can be
seen quite easily
3. Rather than using an inline function, you could have a function m-file that has
the definition for the function. I prefer using inline functions as the the resulting
program can be debugged more easily. Once again, it is a matter of preference.
4. Looking at the code, one gets a sense that the above program would be nicer if
you could replace y=zeros(1,n+1) by y=zeros(1,n). This can be done, but
I advise against it, as it just makes the code more difficult to follow and debug.
All you have to remember is that the number of points equals the number of in-
tervals plus one.

The analytic solution plotted on the same set of axes as the numeric solution is,

Solution of dy/dx=4*x2
1400

1200

1000

800
y

600

400

200

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

101
 The analytic solution is 1336 and the numeric solution is 1286. The accuracy of Eu-
ler’s Method can be improved by choosing a higher resolution of points i.e., increasing
the value of n. With n = 400, the numeric solution is 1333—much better than the pre-
vious result. However, it ca be seen that in order to get an acceptable numeric solution,
the value of n needs to be very high. The other downside of Euler’s Method is that we
do not have a way of assessing the accuracy of our result.

9.2 Vector-valued Functions and Systems of Differ-

ential Equations
M ATLAB only has solvers that can handle first-order differential equations. Thus, if one
has a system of DEs, of mixed order, then the system must be expressable in terms of
a system of FODEs, in order to be able to use the built-in solvers. Also, stiffness2 , chaos
need to be considered when finding solutions of DEs. Since these considerations can
become quite complex, it is something best left to a course on differential equations.

TOOL 9.2 3 Given a system of differential equations, with v distinct variables, where the
system of DEs is of mixed order. Let the highest order of differentiation of vi be ni . We may
generate a vector-valued function y(t) with v-by-n elements as follows,

(n−1) (n−2) (n−1) (n−2)

[y1 (t), y2 (t), · · ·] = [v1 , v1 , . . . , v 1 , v2 , v1 , v2 , · · ·] (9.1)

y(t) may then be used to express the mixed-order system as a system of FODEs by differen-
tiating y(t) and making the necessary substitutions from the original system, and substituting
for all vi from 9.1 i.e.,

(n) (n−1) (n) (n−1)

ẏ(t) = [v1 , v1 , . . . , v 1 , v2 , v1 , v2 , · · ·]
(9.2)
= f (y(t))

⊲ E XAMPLE
Given,
ü(t) = u(t)
v̈(t) = v(t)
Express the mixed-order system as a system of first-order differential equations.

Solution We note that there are 2 distinct variables, u and v. We expect y(t) to have
4 entries as both have a highest order of differentiation of two. The system may be
transformed as follows,
2 Stiffness refers to a type of differential equation or system of differential equations where the

solution of interest varies slightly, but at the same time, nearby solutions vary rapidly. In these
cases, small steps need to be taken so that the solution of interest is found.
3 It is customary to represent differential equations with the independent variable as being time,

hence it is so presented here. The presentation can of course be adapted to situations where the
independent variable is not time.

102
 y1 (t) u̇(t) ẏ1 (t) u(t) y2 (t)
         
 y2 (t)   u(t)   ẏ2 (t)   u̇(t)   y1 (t) 
 y (t)  =  v̇(t)  ⇒  ẏ (t)  =  v(t)  =  y (t) 
3 3 4
y4 (t) v(t) ẏ4 (t) v̇(t) y3 (t)

Hence we may write that,

y2 (t)
 
 y1 (t) 
ẏ(t) = 
y4 (t) 
y3 (t)

It is important to note that we eventually end up with a vector(-valued) function.

One gets a sense that this should be fairly easy to solve in M ATLAB, as M ATLAB has a
preference for vectors.

9.3 Using M ATLAB Built-in ODE Solvers: ODE45

In using the built-in ODE solvers of M ATLAB, the major part is expressing a system of
DEs of mixed-order as a system of FODEs. Once this has been accomplished, the coding
is actually quite easy.
M ATLAB has many ODE solvers. Your best bet is to first try ODE45 as it is fairly
versatile and quite accurate. It is also possible to set things such as tolerances using
odeset, but that shall not be dealt with here. If you for some reason require setting such
parameters, you can look it up in M ATLAB Help. Fear not, for setting such parameters
is an easy task!
Since the transformation of mixed-order systems has been dealt with already, the
example that follows will be based on a system of DEs that are have already been ex-
pressed as FODEs.

⊲ E XAMPLE
Given,
ẏ1 = y2 y3
ẏ2 = −y1 y3
ẏ3 = −0.51y1 y2
Solve the system of 1st order differential equations and show the solution graphically.

% watch the semicolons!

dydt= inline(’[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)]’,’t’,’y’);
[t y]=ode45(dydt,[0 12],[0 1 1]) % return t and y
plot(t,y);
xlabel(’time’);
ylabel(’y’);
title(’ODE45 solution’);

The plot of the solution is,

103
 ODE45 solution
1

0.8

0.6

0.4

0.2

0
y

−0.2

−0.4

−0.6

−0.8

−1
0 2 4 6 8 10 12
time

Comments
1. The order of the variables specified for the inline function (and the m-file, see be-
low) needs to be t then y. The ODE solver will not return the correct solution if
the order is swapped. Also, note that t needs to included as a variable (parame-
ter) even though it is not explicitly part of the right hand sides of the differential
equation system.
2. If you do not specify return arguments for the ODE solvers, the results will be
automatically plotted.
3. The return arguments [t y] could have been any other variable names e.g., [a
b]
4. Once again, the inline function need not be used if you have a preference for
creating a seperate m-file that has the definition of the function dydt. You could
have something like,

function f=dydx(t,y)
dydx=[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)];

The code would have to change accordingly i.e., you will have to use a function
handle. The change would be to remove the line that defines the inline function
and to change the following line to,
[t y]=ode45(@dydx,[0 12],[0 1 1]) % return t and y
Note that you would not have had to have define the function handle, as the
function handle is the same as the name of the function m-file. This is the default
behaviour of M ATLAB.

104
 This chapter has dealt with systems where the time interval is known. In many real-
life applications, one needs to find tf inal for a system. This can be done using events
in M ATLAB. It is an option of the ODE solvers. My initial opinion, was that this is too
advanced an issue to be included here, but I think it will come in useful in your final
year of Chemical Engineering. So get ready. . .

9.4 Using the events Option of ODE Solvers

We shall concern ourselves with the systems that are expressable as FODE systems
where y(t) = k(k is a constant). What we do require is the constraint function, which
will be given in the example that follows, and a function that takes care of the system of
differential equations. Note that the constraint equation needs to be in the form y(t) = 0
where y is the dependent variable. This is because the solvers solve for the zeros of the
event function. It is best explained with an example. There is one other consideration—
the constraint equation is best set up using a seperate m-file. Using a subfunction is
not advised as subfunctions are only have local scope, and thus cannot be called from
outside the function m-file4 .

⊲ E XAMPLE
Given,
ẍ = ẋ + 2, x(0) = 1 and ẋ(0) = 0
When is x(t) = 10?
Solution
The constraint equation is x(t) − 10 = 0. Transforming to a system of FODEs,
· ¸
y1 (t) + 2
ẏ(t) = , y2 (0) = 1; y1 (0) = 0
y1 (t)

The function m-file for the FODE system f.m is,

function ydot=f(t,y)
ydot=[y(1)+2;y(1)];

And the function m-file, c.m for the constraint equation is,

function [constrstop,isterminal,direction]=c(t,y)
constrstop=y(2)-10;
isterminal=1;
direction=[];

The m-file that contains the necessary commands to initialize the solution to the prob-
lem is,

opts = odeset(’events’,@c);
y0 = [0; 1];
[t,y,tfinal] = ode45(@f,[0 Inf],y0,opts);
tfinal %displays tfinal in cmd window
4 If you recall, I mentioned that subfunctions could be made global by the use of function han-

dles. Since I have not presented this use of function handles as yet, it is better to use a seperate
m-file for the constraint equation.

105
plot(t,y(:,2),’-’); %plots the solution NB. x(t)=y(2)
xlabel(’t’)
ylabel(’y’)
title(’Event Handling with ODE45’)
text(1.2, 2, [’tfinal = ’ num2str(tfinal)]) %print tfinal on graph

The result returned is 2.01716306153964 which is very close to the result that one gets
after evaluating the analytical solution. I used dsolve to check whether the answer
is acceptable. For t = 2.01716306153964, x(t) gives 9.99961282914699, which is fairly
close to 10. The plot is of the solution is,

Event Handling with ODE45

10

6
y

2 tfinal = 2.0172

1
0 0.5 1 1.5 2 2.5
t

Comments
1. For c.m, the output argument variable names isterminal and direction are
fixed i.e., they cannot be named differently. constrstop could be named dif-
ferently. If isterminal is 1 then, the ODE solver stops when the a (i.e., any)
zero of the constraint equation has been found. If isterminal is 0, then the
ODE solver does not stop when a zero has been reached. with direction=[]
or 0, all zeros are computed. If direction=1, then zeros where the constraint
function increases are computed and if If direction=-1, the zeros where the
constraint function decreases are computed.
2. Note that the time interval input for the ODE solver is from 0 to ∞.
The rest of the code is quite self-explanatory. As you can see using M ATLAB to solve
differential equations is not that difficult. . . it is setting up the system of equations from
the description of a situation that usually poses a problem!

106
Bibliography

[1] Hahn BD, Essential M ATLAB for Scientists and Engineers, 3rd ed.,Maskew
Miller Longman: 2002
[2] Rawatlal R, An Introduction to Programming in M ATLAB With Applications in
Chemical Engineering, University of Kwazulu-Natal: April 2003
[3] Online M ATLAB help for M ATLAB version 6.5 release 13

107
Index

, 62 fopen, 26, 27
’, 20 for, 17, 46, 50, 94
*, 16 fprintf, 23, 25, 26
-ascii, 29 gcbo, 78
-double, 29 get, 71, 77
-tab, 29 ginput, 37
.m, 13 guide, 70
/, 16 hold, 34
GUIDE, 69, 70, 76 if, 18, 19, 57
NaN, 57 importdata, 28
ODE45, 99, 103 inline, 59–61
&, 18, 58 inputargs, 14
˜, 18 input, 23, 24
and, 18 interp1, 91, 96
axis, 34 interp, 96
break, 17 intersect, 58
celldisp, 47 int, 61
cellplot, 47 isterminal, 106
cell, 44, 46 i, 17
char, 27 j, 17
children, 36 legend, 35, 36
cla, 85 length, 21
clear all, 15 linspace, 19, 20, 100
continue, 17 load, 26, 29
diff, 61 logical, 53
disp, 23, 24 max, 77
do, 17 meshgrid, 20, 21
dsolve, 106 min, 77
elseif, 19, 53, 58 not, 18
else, 18, 19 num2str, 85
end, 17–19 odeset, 103
eps, 56 ones, 44
evalin, 64, 79 or, 18
eval, 61 otherwise, 18, 19
events, 105 outputargs, 14
ezplot, 33, 61 plot3, 33
ezsurf, 33 plot, 31, 33
fclose, 26 polyfit(x,y,n), 95
feval, 59, 61, 62 polyfit, 91, 95, 96
find, 34, 53–55 polyval, 91, 95

108
repmat, 19, 20
reshape, 46
return, 17
rmfield, 50
save, 26, 29
set, 71, 77
size(m,1), 21
size(m,2), 21
size, 19, 21
sqrt, 62
str2double, 85
str2num, 27
struct, 48
subplot, 34
surf, 33
switch, 18, 19
textread, 28
title, 36
transpose, 19, 20
varargin, 43
varargout, 43
while, 17
xlabel, 36
ylabel, 36
zeros, 44
zlabel, 36

109