MATLAB Manual 05
MATLAB Manual 05
MATLAB Manual 05
Multiplying a 1 x n (row) vector x by a matrix A, results in a 1 x k (row) vector:
1
n
j p
p
y x
=
=
pj
A
Secondly, an n x k matrix A can be multiplied by a matrix B, only if B has k rows, i.e. B is k x m
(with m an arbitrary value). As a result, you get an n x m matrix C, such that A B = C, where:
1
k
ij ip pj
p
C A
=
= B
>> b = [1 3 -2];
>> B = [1 -1 3; 4 0 7]
B =
1 -1 3
4 0 7
>> b * B % not possible: b is 1-by-3 and B is 2-by-3
??? Error using ==> *
Inner matrix dimensions must agree.
>> b * B' % this is possible: a row vector multiplied by a matrix
ans =
-8 -10
>> B' * ones(2,1)
ans =
5
-1
10
25-Jan-05
28
>> C = [3 1; 1 -3];
>> C * B
ans =
7 -3 16
-11 -1 -18
>> C.^3 % this is element-by-element power
ans =
27 1
1 -27
>> C^3 % this is equivalent to C*C*C
ans =
30 10
10 -30
>> ones(3,4)./4 * diag(1:4)
ans =
0.2500 0.5000 0.7500 1.0000
0.2500 0.5000 0.7500 1.0000
0.2500 0.5000 0.7500 1.0000
)
Exercise 17
Given: vectors x = [1 3 7] and y = [2 4 2], and matrices A = [3 1 6; 5 2 7] and B = [1 4;7 8;2 2],
determine which of the following statements will correctly execute (and if not, try to understand
why); also provide the result:
x + y
x + A
x + y
A [xy]
[x; y] + A
[x ; y]
[x; y]
A 3
A + B
B+ A
B * A
A. * B
A .* B
2 * B
2 .* B
B ./ x
B ./ [xx]
2 / A
ones(1,3) * A
ones(1,3) * B
)
Exercise 18
Let A be a random 3x3 matrix, and b a random 3x1 vector (use the command randn). Given that Ax
= b, then x = A
-1
b. Find x. Explain what the difference is between the operators \, / and the
command inv (look at Table 4). Having found x, check whether Ax - b is close to a zero vector.
5 Text and cell arrays
5.1 Character strings
Although MATLAB mainly works with numbers, it might occasionally be necessary to operate on
text. Text is saved as a character string, i.e., a vector of ASCII values which are displayed as their
character string representation. Since text is a vector of characters, it can be addressed in the same
way as any vector. For example:
>> t = 'This is a character string'
t =
This is a character string
>> size(t)
ans =
1 26
>> whos
Name Size Bytes Class
ans 1x2 16 double array
t 1x26 52 char array
Grand total is 28 elements using 68 bytes
25-Jan-05
29
>> t(10:19)
ans =
character
>> t([2,3,10,17]) % note brackets [ ] inside parentheses ( )
% to denote collection of indices
ans =
hi t
To see the underlying ASCII representation, the string has to be converted using the command
double or abs:
>> double(t(1:12))
ans =
84 104 105 115 32 105 115 32 97 32 99 104
The function char provides the reverse transformation:
>> t(16:17)
ans =
ct
>> t(16:17)+3
ans =
102 119
>> t(16:17)-3
ans =
96 113
>> char(t(16:17)-2)
ans =
ar
)
Exercise 19
Perform the commands described. Use string t to form a new string u that contains the word
character only. Convert string u to form the same word spelled backwards, i.e., u1 =
retcarahc.
With findstr, a string can be searched for a character or a combination of characters. Some
examples on the use of different string functions are given below:
>> findstr(t, 'c')
ans =
11 16
>> findstr(t, 'racter')
ans =
14
>> findstr(t, u)
ans =
11
>> strcat(u,u1)
ans =
characterretcarahc
>> strcmp(u,u1)
ans =
0
>> q = num2str(34.35)
q =
34.35
>> z = str2num(q)
z =
34.3500
25-Jan-05
30
>> whos q z
Name Size Bytes Class
q 1x5 10 char array
z 1x1 8 double array
Grand total is 6 elements using 18 bytes
>> t = str2num('1 -7 2')
t =
1 -7 2
>> t = str2num('1 - 7 2')
t =
-6 2
>> A = round(4*rand(3,3))+0.5;
>> ss = num2str(A)
ss =
4.5 2.5 2.5
1.5 4.5 0.5
2.5 3.5 3.5
>> whos ss
Name Size Bytes Class
ss 3x27 162 char array
Grand total is 81 elements using 162 bytes
>> ss(2,1), ss(2,1:27)
ans =
1
ans =
1.5 4.5 0.5
>> ss(1:2,1:3)
ans =
4.5
1.5
)
Exercise 20
Perform the commands shown in the example above. Define another string, e.g. s = Nothing
wastes more energy than worrying and exercise with findstr.
)
Exercise 21
Become familiar with num2str and str2num. Check, e.g., what happens with ss = num2str([1 2
3; 4 5 6]) or q = str2num(1 3.25; 5 5.5). Analyze also the commands str2double and
int2str (use help).
5.2 Text input and output
The input command can be used to prompt (ask) the user for numeric or string input:
>> myname = input(Enter your name: , s);
age = input(Enter your age: );
There are two common text output functions: disp and fprintf. The disp function only displays
the value of one parameter, either a numerical array or a string (the recognition is done
automatically). For example:
>> disp(This is a statement.) % a string
This is a statement
>> disp(rand(3))
0.9501 0.4860 0.4565
0.2311 0.8913 0.0185
0.6068 0.7621 0.8214
The fprintf function (familiar to C programmers) is useful for writing data to a file or screen when
precise formatting is important. Try for instance (an explanation will follow):
25-Jan-05
31
>> x = 2;
>> fprintf('Square root of %g is %8.6f.\n' , x, sqrt(x));
Square root of 2 is 1.414214.
Formally, fprintf converts, formats and prints its arguments to a file or displays them on screen
according to a format specification. For displaying on screen, the following syntax is used:
fprintf (format, a, ...)
The format string contains ordinary characters, which are copied to the output, and conversion
specifications, each of which converts and prints the next successive argument to fprintf. Each
conversion specification is introduced by the character % and ended by a conversion character.
Between the % and the conversion character there may appear:
a minus sign; controlling the alignment within the field defined.
a digit string specifying a minimum field length; The converted number will be printed in the
field at least this wide, and wider if necessary.
a period which separates the field width from the next digit string;
a digit string - the precision which specifies the maximum number of characters to be printed
from a string or the number of digits printed after the decimal point of a single or double type.
Format specification is given below in terms of the conversion characters and their meanings:
Command The argument
d is converted into decimal notation;
u is converted into unsigned decimal notation
c is taken to be a single character;
s is a string
e is taken to be a single or double and converted into decimal notation of the form:
[-]m.nnnnnnE[]xx, where the length of ns is specified by the precision;
f is taken to be a single or double and converted into decimal notation of the form:
[-]mmm.nnnnn, where the length of ns is specified by the precision;
g is specified by e or f, whichever is shorter; non-significant zeros are not printed;
Table 5, Format specification
The special formats \n, \r, \t, \b can be used to produce next line (i.e. makes sure that the next
command or output will be written on the next line), carriage return, tab and backspace respectively.
Use \\ to produce a backslash character and %% to produce the percent character.
Analyze the following example:
>> fprintf('Look at %20.6e!\nSee help fprintf', 1000*sqrt(2))
Look at 1.414214e+003!
See help fprintf
>> fprintf('Look at %-20.6e!See help fprintf\n', 1000*sqrt(2))
Look at 1.414214e+003 !See help fprintf
For both commands, the minimum field length is 20 and the number of digits after the decimal point
is 6. In the first case, the value of 1000*sqrt(2) is padded on the right, in the second case, because of
the -, it appears on the left.
)
Exercise 22
Try to understand how to use the input, disp and fprintf commands. For instance, try to read a
vector with real numbers using input. Then, try to display this vector, both by calling disp and
formatting an output by fprintf . Make a few variations. Try the same with a string.
)
Exercise 23
Study the following examples in order to become familiar with the fprintf possibilities and
exercise yourself to understand how to use these specifications:
>> str = 'life is beautiful';
>> fprintf('My sentence is: %s\n', str); % note the \n format
My sentence is: life is beautiful
25-Jan-05
32
>> fprintf('My sentence is: %30s\n', str);
My sentence is: life is beautiful
>> fprintf('My sentence is: %30.10s\n', str);
My sentence is: life is be
>> fprintf('My sentence is: %-20.10s\n', str);
My sentence is: life is be
>>
>> name = 'John';
>> age = 30;
>> salary = 6130.50;
>> fprintf('My name is %4s. I am %2d. My salary is $ %7.2f.\n', name, age,
salary);
My name is John. I am 30. My salary is $ 6130.50.
>>
>> x = [0, 0.5, 1]; y = [x; exp(x)];
>> fprintf('%6.2f %12.8f\n',y);
0.00 1.00000000
0.50 1.64872127
1.00 2.71828183
>> fprintf('%6.1e %12.4e\n',y);
0.0e+000 1.0000e+000
5.0e-001 1.6487e+000
1.0e+000 2.7183e+000
>>
>> x = 1:3:7; y = [x; sin(x)];
>> fprintf('%2d %10.4g\n',y);
1 0.8415
4 -0.7568
7 0.657
Warning: fprintf uses its first argument to decide how many arguments follow and what their
types are. If you provide a wrong type or there are not enough arguments, you will get nonsense for
an answer. So, be careful with formatting your output. Look, e.g. what happens in the following case
(age is not provided):
>> fprintf('My name is %4s. I am %2d. My salary is $ %7.2f.\n', name,
salary);
My name is John. I am 6.130500e+003. My salary is $
Remark: The function sprintf is related to fprintf, but writes to a string instead. Analyze the
example:
>> str = sprintf('My name is %4s. I am %2d. My salary is $ %7.2f.\n', name,
age, salary)
str =
My name is John. I am 30. My salary is $ 6130.50.
)
Exercise 24
Define a string s = A woodchuck couldnt chuck any amount of wood because a
woodchuck cant chuck wood, but if a woodchuck could chuck and would chuck
wood, how much wood would a woodchuck chuck?. Exercise with findstr, i.e., find all
appearances of the substring wood, o, uc or could. Try to build a new string ss by
concatenation of separate words in a few ways; use strcat, disp, sprintf or fprintf.
5.3 Cell arrays
Cell arrays are arrays whose elements are cells. Each cell can contain any data, including numeric
arrays, strings, cell arrays, etcetera. For instance, one cell can contain a string, another a numeric
array, etc. Below, there is a schematic representation of an exemplar cell array:
25-Jan-05
33
1 + 3i
cell 1,3
3 1 -7
7 2 4
0 -1 6
7 3 7
cell 1,1
5900 6000
35
John Smith
cell 1,2
this is a text
cell 2,1
implies effort
Living
cell 2,3
Bye -1 0
0 -1
[7,7] Hi
cell 2,2
Cell arrays can be built up by assigning data to each cell. The cell contents are accessed by brackets
{}, e.g.
>> A(1,1) = {[3 1 -7; 7 2 4; 0 -1 6; 7 3 7]};
>> A(2,1) = {'this is a text' };
>> B(1,1) = {'John Smith'}; B(2,1) = {35}; B(3,1) = {[5900 6000]}
B =
'John Smith'
[ 35]
[1x2 double]
>> A(1,2) = {B}; % cell array B becomes one of the cells of A
>> C(1,1) = {'Hi'}; C(2,1) = {[-1 0; 0 -1]}; C(1,2) = {[7,7]}; C(2,2) =
{'Bye'};
>> A(2,2) = {C}; %cell array C becomes one of the cells of A
>> A
A =
[4x3 double] {3x1 cell}
'this is a text' {2x2 cell}
>> A(2,2) %access the cell but not the contents
ans =
{2x2 cell}
>> A{2,2} %use {} to display the contents
ans =
'Hi' [1x2 double]
[2x2 double] 'Bye'
>> A{2,2}{2,1} %take the (2,1) element of the cell A{2,2}
ans =
-1 0
0 -1
There are also two useful functions with meaningful names: celldisp and cellplot. Use help to
learn more.
The common application of cell arrays is the creation of text arrays. Consider the following
example:
>> Month = {'January';'February';'March';'April';'May';'June';'July';...
'August';'September';'October';'November';'December'};
>> fprintf('It is %s.\n', Month{9});
It is September.
25-Jan-05
34
)
Exercise 25
Exercise with the concept of a cell array, first by typing the examples presented above. Next create a
cell array W with the names of week days, and, using the command fprintf, display on screen the
current date with the day of the week. The goal of this exercise is also to use fprintf with a format
for a day name and a date, in the spirit of the above example.
6 Control flow
Since variables are updated by commands (e.g., b=10; c=b+1), the order in which the commands are
executed is important. Control flow ensures that commands will be executed in a suitable order. This
chapter deals with ways of composing commands to achieve different control flows.
6.1 Expressions and logical expressions
To create different control flows logical expressions are a means to decide which group of
subcommands to execute. These expressions result in a logical value: TRUE or FALSE. In MATLAB
the result of a logical expression is 1 if it is true (e.g., 1 == 2-1) and 0 if it is false.
First check section 2.3 on MATLAB syntax.
6.2 Sequential commands
The most common control flow is the sequential composition of two (or more) commands: C
1
; C
2
;
denotes that command C
1
is executed before command C
2
.
6.3 Conditional commands
A conditional command has a number of subcommands (each consisting of a sequential group of
commands). Exactly one group of subcommands (starting with if, else, or elseif, and
ending with else, elseif, or end) is selected for execution. The choice among the various
candidate subcommands is based upon the truth-value yielded by a logical expression Elogical (e.g.
a>b). In the example below the disp function if frequently used.
if (a >= b)
max = a;
di
else
sp (a >= b);
max = b;
d
end
isp (b > a);
Example Syntax
if E
al
logoc
; C1
else
C
end
2;
if (height > 190)
max = a;
dis
elseif (height > 170)
p (a >= b);
dis
elseif (height < 150)
p(tall);
di
else
sp (small);
d
end
isp (average);
Example
if E1
C1;
elseif E2
% executed if E1 is True
C2
else
; % executed if E2 is True
C3; % executed if no other
end
% expression is True
Syntax
)
Exercise 26
Create a script that asks for a number and computes the drag coefficient C, depending on the
Reynolds number (make use of elseif). The command input might be useful here (use help if
needed).
25-Jan-05
35
0.7
0, 0
24/ , (0, 0.1]
(24/ )(1 0.14 ), (0.1, 1000]
0.43, (1000, 500000]
0.19 80000/ , 500000
N
N N
C N N N
N
N N
= +
>
Check whether your script works correctly. Compute the drag coefficient for e.g. N=-3e3, 0.01,
56, 1e3, 1e6 (remember that e.g. 3e3 is a MATLAB notation of 3*10
3
).
)
Exercise 27
Write a script that asks for an integer and checks whether it can be divided by 2 or 3. Consider all
possibilities, such as: divisible by both 2 and 3, divisible by 2 and not by 3 etc. (use the command
rem).
It is also possible to select a group of subcommands on the basis of values other than truth values
(e.g., scalars or a string). The value yielded by the expression E must equal the value of one of the
literals L
i
and then the corresponding subcommand C
i
is chosen.
method
swi ethod
= lin
tch m
case {lin,lin2}
d
case {near}
isp(Method is linear.);
disp(Met
otherwise
hod is nearest.);
end
disp(Unknown method.);
Example
swit ch E
case L 1
; % executed if L1 equals E C1
case {L , L , } 2 3
uted if L2 or L3 equals E C2; % exec
otherwise
C ; % executed if no literal Li 3
% equals E
end
Syntax
)
Exercise 28
Assume that the months are represented by the numbers from 1 to 12. Write a script that asks you to
provide a month and returns the number of days in that particular month. Alternatively, write a script
that asks you to provide a month name (e.g. June) instead of a number. Use the switch-
construction.
)
Exercise 29
Write a MATLAB function that implements :
if 0<x<7
y=4
elseif 7<x<55
x
y=
else
-10x
y
end
=333
This is not MATLAB syntax
Verify its functionality using
x=-1, x= 5, x=30, and x=56
6.4 Iterative commands
An iterative command (commonly known as a loop) has a subcommand (the loop body) that is to be
executed repeatedly, and some kind of phrase that determines when iteration will cease. We can
classify iterative commands according to whether the number of iterations is determined in advance
(definite) or not (indefinite).
25-Jan-05
36
6.4.1 Definite iteration
Definite iteration is characterized by the use of a control variable. The loop body is executed with
the control variable taking each of a predetermined sequence of values. This sequence of values is
called the control sequence.
note 1: the control sequence can be a vector.
note 2: the step can have a negative value.
for a = 0 : 0.5 : 4
d
end
isp(a^2)
Example
for a = first : step : last
C
end
1;
Syntax
v =
for a = v
[0 : 0.5 : 4];
d
end
isp(a^2);
Example Syntax
for a = vector
C
end
1;
In MATLAB, the for command is self-contained: control variables are created within the loop body.
Thus the control variable has no value (a in the example) outside the loop (it does not exist). The
loop body cannot cause an assignment to the control variable - it is in fact a constant!
6.4.2 Indefinite iteration
This concerns loops in which the number of iterations is not determined in advance. The while
command consists of a loop body preceded by a truth-valued expression E
logical
(the loop condition)
that controls whether iteration is to be (started or) continued.
note 1: The number of executions of a command can also be 0!
note 2: if Elogical is always true the loop will not have a (natural) end!
sum = 0;
delta = 1;
iter
while(delta > eps)
= 0;
sum = sum + delta;
delta = delta / 2;
end
iter = iter + 1;
disp(sum); disp(iter);
Example Syntax
while Elogical
C
end
1;
Given two vectors x and y, an example use of the loop construction is to create a matrix A whose
elements are defined, e.g. as A
ij
= x
i
y
j
. Enter the following commands in a script:
n = length(x); m = length(y);
for i = 1:n
for j = 1:m
A(i,j) = x(i)*y(j);
end
end
Now create A for x = [1 2 -1 5 -7 2 4] and y=[3 1 -5 7]. Note that A is of size n-by-m. The same
problem can be solved by using two while loops:
n = length(x); m = length(y); i = 1 % initialize i
while i<=n
j=1; % initialize j
while j <= m
A(i,j) = x(i)*y(j);
j=j+1 % increment j: it does not happen automatically!
end
i=i+1 % increment i
25-Jan-05
37
end
)
Exercise 30
Use a loop construction to carry out the computations. Write short scripts.
1. Given the vector x = [1 8 3 9 0 1], create a short set of commands that will:
add up the values of the elements (check with sum);
computes the running sum (for element j, the running sum is the sum of the elements from 1
to j; check with cumsum).
2. Given x = [4 1 6 -1 -2 2] and y = [6 2 -7 1 5 -1], compute matrices whose elements are created
according to the following formulas:
a
i
= x
i
y
i
and add up the elements of a;
b
ij
= x
i
/(2 + x
i
+ y
j
);
c
ij
= 1 / max (x
i
,y
j
);
3. Determine the sum of the first 50 squared numbers with a control loop.
4. Write a script to find the largest value n such that the sum:
3 3
1 2 ... n + + +
3
is less than
1000.
6.5 Evaluation of logical and relational expressions in the control flow structures
The relational and logical expressions may become quite complicated. It is not difficult to operate
on them if you understand how they are evaluated. To explain this in more detail, let us consider the
following example:
if (~isempty(data)) & (max(data) < 5)
...
end
This construction of the if-block is necessary to avoid comparison if data happens to be an empty
matrix. In such a case you cannot evaluate the right logical expression and MATLAB returns an error.
The & operator returns 1 only if both expressions: ~isempty(data) and max(data) < 5 are
TRUE, and 0 otherwise. When data is an empty matrix, the next expression is not evaluated since
the whole &-expression is already known to be false (lazy evaluation). The second expression is
checked only if data is a non-empty matrix. Dont forget to put logical expression units between
brackets to avoid unexpected evaluations!
)
Exercise 31
Consider the following example:.
max_iter = 50;
tolerance = 1e-4;
iter = 0;
xold = 0.1;
x = 1;
while (abs (x-xold) > tolerance) & (iter < max_iter)
xold = x;
x = cos(xold);
iter = iter + 1;
end
This short program tries to solve the equation cos(x) = x (x is the solution found). Make the script
solve_cos from the presented code. Note that the while-loop is repeated as long as both conditions
are true. If either the condition (|x-xold| <= tolerance) or (iter >= max_iter) is fulfilled then the loop
is quitted. Run the script solve_cos for different tolerance parameters, e.g. 1e-3, 1e-6, 1e-8,
1e-10, etc. Use format long to check more precisely how much the found x is really different
from cos(x). For each tolerance value check the number of performed iterations (iter value).
25-Jan-05
38
)
Exercise 32
Create the script solve_cos2, which is equal to the given, replacing the while-loop condition by:.
while (abs (x-xold) > tolerance) | (iter < max_iter)
Try to understand the difference and confirm your expectations by running solve_cos2. What
happens to iter?
7 Visualization
MATLAB can be used to visualize the results of an experiment. To that end, you should define your
variables such that each of them contains all values of one parameter to plot.
7.1 Simple plots
With the command plot, a graphical display can be made. For a vector y, plot(y) draws the
points [1,y(1)], [2,y(2), ..., [n, y(n)] and connects them with a straight line. plot(x,y) does the
same for the points [x(1),y(1)], [x(2),y(2)], ..., [x(n),y(n)]. Note that x and y have to be both either
row or column vectors of the same length (i.e., the same number of elements).
The commands loglog, semilog and semilogy are similar to plot, except that they use either one
or two logarithmic axes. Type the following commands, after predicting the result:
>> x = 0:10;
>> y = 2.^x; % the same as y = [1 2 4 8 16 32 64 128 256 512 1024]
>> plot(x,y); % to get a graphic representation
>> semilogy(x,y); % to make the y-axis logarithmic
As you can see, the same figure is used for both commands. The previous function is removed as
soon as the next is displayed. The command figure gives you an extra figure-window. Repeat the
previous commands, but generate a new figure before plotting the second function, so that you can
see both functions in separate windows. You can also switch back to a figure using figure(n),
where n is its number.
To plot a graph of a function, it is important to sample the function sufficiently well. Compare the
following examples:
>> n = 3;
>> x = 0:1/n:3; % coarse sampling
>> y = sin(5*x);
>> plot(x,y);
0 0.5 1 1.5 2 2.5 3
-1
-0.5
0
0.5
1
>> n = 30;
>> x = 0:1/n:3; % fine sampling
>> y = sin(5*x);
>> plot(x,y)
0 0.5 1 1.5 2 2.5 3
-1
-0.5
0
0.5
1
The plot command uses a black solid line by default. It is possible to change both the style and
colour, e.g.:
>> x = 0:0.4:3; y = sin(5*x);
>> plot(x,y,r--) % produces a dashed red line.
25-Jan-05
39
Symbol Colour Symbol Line style
r red .,o point, circle
g green * star
b blue x,+ x-mark, plus
y yellow - solid line
m magenta -- dash line
c cyan : dot line
k black -. dash-dot line
Table 6, Plot colours and styles
The third argument of plot specifies the colour (optional) and the line style. Table 6 shows a few
possibilities, help plot shows all. To add a title, grid and to label the axis, one uses:
>> title(Function y = sin(5*x));
>> xlabel(x-axis);
>> ylabel(y-axis);
>> grid % remove grid by calling grid off
)
Exercise 33
Make a plot connecting the coordinates: (2,6), (2.5, 18), (5,17.5), (4.2, 12.5) and (2,12) by a line.
)
Exercise 34
Plot the function y=sin(x) + x - x cos(x) in two separate figures for the intervals:
0<x<30 and -100<x<100. Add a title and axis description.
7.2 Several functions in one figure
There are different ways to draw several functions in the same figure-window. The first one uses the
command hold on. After this command, all functions will be plotted in the same figure until the
command hold off is used. When more than one function is plotted in a window, it is useful to use
different symbols and colours. An example is:
>> x1 = 1:.1:3.1; y1 = sin(x1);
>> plot(x1 , y1, md); % d for diamond symbol
>> x2 = 1:.3:3.1; y2 = sin(-x2+pi/3);
>> hold on
>> plot (x2,y2,k*-.)
>> plot(x1,y1,m-)
>> hold off
A second method to display a few functions in one figure is to extend the list of arguments handled
by the plot command. The next command will produce the same output as the commands in the
previous example.
>> x1 = 1:.1:3.1; y1 = sin(x1);
>> x2 = 1:.3:3.1; y2 = sin(-x2+pi/3);
>> plot (x1, y1,md, x2, y2,k*-., x1, y1,m-)
To make the axis better fitting the curves, perform
>> axis([1,3.1,-1m1])
The same can be achieved by axis tight. It might be also useful to exercise with axis options (see
help), e.g. axis on/of, axis equal, axis image or axis normal. A descriptive legend can be
included with the command legend, e.g.:
>> legend (sin(x), sin(-x+pi/3);
25-Jan-05
40
It is also possible to produce a few subplots in one figure window. With the command subplot(p,
r, figurenumber), the window can be horizontally and vertically divided into p x r subfigures,
which are counted from 1 to pr, row-wise, starting from the top left. Note that the commands: plot,
title, grid, etc. work only on the current subfigure.
>> x = 1:.1:4;
>> y1 = sin(3*x);
>> y2 = cos(5*x);
>> y3 = sin(3*x).*cos(5*x)
>> subplot(2,2,1); plot(x,y1,m-); title(sin(3*x))
>> subplot(2,2,2); plot(x,y2,g-); title(cos(5*x))
>> subplot(2,1,2); plot(x,y3,k-); title(sin(3*x) * cos(5*x))
0 1 2 3
-1
-0.5
0
0.5
1
sin(3*x)
0 1 2 3
-1
-0.5
0
0.5
1
cos(5*x)
0 0.5 1 1.5 2 2.5 3
-1
-0.5
0
0.5
1
sin(3*x) * cos(5*x)
)
Exercise 35
subplot(2,2,1): if the plot window is a figure matrix with dimensions n x k (row x
column) this would be the subfigure at the first position of a 2x2 matrix.
Plot the functions f(x) = x, g(x) = x^3, h(x) = e^x and z(x) = e^(x^2)) over the interval [0,4] using
both the normal and the log-log scales. Use an appropriate sampling to get smooth curves. Describe
your plots by using the functions: xlabel, ylabel, title and legend.
Command Result
grid on/off adds a grid to the plot at the tick marks or removes it
axis([xmin xmax ymin ymax]) sets the minimum and maximum values of the axes
axis([xmin xmax ymin ymax]) sets the minimum and maximum values of the axes
box off/on removes the axes box or shows it
xlabel(text) plots the label text on the x axis
ylabel(text) plots the label text on the y axis
title(text) plots a title above the graph
text(x,y,text) adds text at the point (x,y)
gtext(text) adds text at a manually (with the mouse) indicated point
legend(fun1,fun2) plots a legend box (move it with the mouse)
legend off deletes the legend box
clf clears the current figure
subplot creates a subplot in the current figure
Table 7, Plot/print commands
25-Jan-05
41
7.3 Printing
Before printing a figure, you might want to add some information, such as a title, or change
somewhat in the lay-out. Table 7 shows some of the commands that can be used.
)
Exercise 36
Plot the functions y1 = sin(4x), y2 = x cos(x), y3 = (x+1)^-1sqrt(x) for x = 1:0.25:10, and plot a
single point (x,y) = (4,5), all in one figure. Use different colours and styles. Add a legend, labels for
both axes and a title. Add also a text to the single point saying: single point. Change the minimum
and maximum values of the axes such that one can look at the function y3 in more detail.
When you like the displayed figure you can print it to paper. The easiest way is to click File in the
menu-bar and to choose Print. If you click OK in the print window, your figure will be sent to the
printer indicated there.
There also exists a print command, which can be used to send a figure to a printer or output it to a
file. You can optionally specify a print device (i.e., an output format such as tiff or postscript) and
options that control various characteristics of the printed file (i.e., which figure to print, etc.). You
can also print to a file if you specify the file name. If you do not provide an extension, print adds
one. Since there are many parameters they will not be explained here (click help print to learn
more) instead, try to understand the examples:
>> print -dwinc % print current figure to current printer in colour
>> print -f1 -deps myf.eps % print Figure no.1 to file myf.eps in black
>> print -f1 -depsc myg.eps % print Figure no.1 to file myg.eps in colour
>> print -dtiff myfile.tif % print current figure to file myfile.tif
>> print -f2 -djpeg myfile.jpg % print Figure no.2 to the file myfile.jpg
8 Numerical analysis
Numerical analysis can be used whenever it is very difficult or impossible to determine the
analytical solution. MATLAB can be used to find, for instance, the minimum, maximum or integral of
a particular function.
There are two ways to describe data with an analytical function:
curve fitting or regression; finding a smooth curve that best fits (approximates) the data
according to a criterion (e.g., best least square fit). The curve does not have to go through any of
the data points.
interpolation; the data are assumed to be error-free: desired is a way to describe what happens
between the data points. In other words, given a finite set of points, the goal of interpolation is
to find a function from a given class of functions (e.g. polynomials) which passes through the
specified points.
The figure below illustrates the difference between regression (curve fitting) and interpolation:
0 1 2 3 4 5 6 7 8 9 10
-2
0
2
4
6
8
10
linear interpolation
fitting a third-order polynomial
25-Jan-05
42
In this section we will briefly discuss some issues on curve fitting, interpolation, function
evaluation, integration and differentiation. It is also possible to solve differential equations (see
Matlab books mentioned in the Introduction).
8.1 Curve fitting
MATLAB interprets the best fit of a curve as least squares curve fitting. The curves used are
restricted to polynomials. Which the command polyfit any polynomial can be fitted to the data.
polyfit (x,y,n) finds the coefficients of a polynomial of degree n that fits the data (finds the
linear combination between x and y). Lets start with finding a linear regression of some data:
>> x = 0:10;
>> y = [-0.10 .24 1.02 1.58 2.84 2.76 2.99 4.05 4.83 5.22 7.51]
>> p = polyfit(x,y,1) % find the fitting polynomial of order 1
p =
0.6772 -0.3914
The output of polyfit is a row vector of the polynomial coefficients; the solution of this example is
therefore y = 0.6772x - 0.3914.
)
Exercise 37
Execute the above example. Plot the data and the fitted curve in a figure. Determine the 2nd and 9th
order polynomial of these data as well, and display the solution in a figure.
Hint: x1 = linspace(xmin, xmax, n); results in a vector with n elements evenly distributed
from xmin till xmax. z = polyval(p,x1); calculates the values of the polynomial for every
element of x1.
)
Exercise 38
Define x and y as follows:
>> x = -5:5;
>> y = 0.2 * x.^3 - x + 2; % a 3rd order polynomial
>> y = y + randn(11,1); % add some noise to y
Determine the 3rd order polynomial fitting y (note that because of added noise the coefficients are
different from those originally chosen). Try some other polynomials as well. Plot the results in a
figure.
8.2 Interpolation
The simplest way to examine various interpolations is by displaying the function plots with the
command plot. When the neighbouring data points are connected with straight lines, this is called
linear interpolation. When more data points are taken into account, the effect of interpolation
becomes less visible. Analyse:
>> x1 = linspace(0,2*pi, 2);
>> x2 = linspace(0,2*pi,4);
>> x3 = linspace(0,2*pi,16);
>> x4 = linspace(0,2*pi,256);
>> plot(x1, sin(x1), x2, sin(x2),x3,sin(x3),x4, sin(x4))
>> legend(2 points,4 points,16 points,256 points)
There also exist MATLAB functions that interpolate between points, e.g., interp1 or spline, as a
1D interpolation and interp2, as a 2D interpolation. Perform the following commands:
>> N = 50;
>> x = linspace(0,5,N);
>> y = sin(x) .* sin(6*x);
>> subplot(2,1,1); plot(x,y);
>> hold on
>> p = randperm(N);
25-Jan-05
43
>> pp = p(1:round(N/2)); % choose randomly N/2 indices of points from [0,5]
>> pp = sort(pp); % sort the indices
>> xx = x(pp); % select the points
>> yy = y(pp);
>> plot(xx,yy,ro-) % this is a coarse version
>> yn = interp1(xx,yy,x,nearest); % nearest neighbour interpolation on
all x points
>> plot(x,yn,g)
>> axis tight
>> legend(Original,Crude version); % Nearest neighbour interpolation
>> subplot(2,1,2); plot(xx,yy,ro-);
>> hold on
>> yc = interp1(xx,yy,x,linear); % spline interpolation
>> plot(x,yc,g);
>> ys = spline(xx,yy,x);
>> plot(x,ys,k)
>> axis tight
>> legend(Crude version, Linear interpolation, Spline interpolation);
8.3 Evaluation of a function
As said before, the command plot plots a function by connecting defined data points.
MATLAB Command Result
fplot(f,[min_x max_x] plots the function f on [min_x, max_x]
fminbnd(f, min_x, max_x, options) returns the x-value for which the function f
reaches a minimum on [min_x, max_x]
fzero(f, x0) returns the x-value for which the function f is
zero in the neighborhood of x0
Table 8, Evaluation of a function
When a function is constant and not very interesting over some range, and next shows an
unexpected wild behaviour over another, it might be misinterpreted by using plot. Then, the
command fplot is more useful. An example is:
>> x = 0: pi/8:2*pi; y=sin(8*x); plot(x,y,b)
>> hold on
>> fplot(sin(8*x),[0 2*pi],r)
>> title(sin(8*x))
>> hold on
Execution of these commands results in a figure with an almost straight blue line and a red sine.
fplot can also be used to plot a defined function:
>> f = 'sin(8*x)';
>> fplot(f,[0 2*pi],'r');title(f)
Let:
( )
( ) ( )
2 2
1 1
0.1 0.1 1 0.
f x
x x
= +
+ + 1
.
The point at which ( ) f x takes its minimum is found using fminbnd. By default, the relative error is
1e-4. It is possible however, to obtain a higher accuracy, by specifying an extra parameter while
calling fminbnd.
>> format long % you need this format to see the change in accuracy
>> f = '1./((x-0.1).^2 + 0.1) + 1. / ((x-1).^2 + 0.1)';
>> fplot(f,[0 2]);
>> options = optimset('TolX',1e-4);
>> [xm1,fm1] = fminbnd(f,0.3,1, options)
>> options = optimset('TolX',1e-12);
>> [xm2,fm2] = fminbnd(f,0.3,1, options)
25-Jan-05
44
>> [xm1,xm2] % compare the two answers
Important: Note that format does not change the accuracy of the numbers or calculations; it just
changes the numbers of digits displayed on the screen.
)
Exercise 39
Determine the maximum of the function f(x) = x cos(x) over the interval 10<x<15. Hint: there is no
function fmax. Find the minimum of -f(x), instead.
)
Exercise 40
Find the zero, i.e. x
0
, of the functions f(x) = cos(x) and g(x) = sin(2x) around the point x=2. Use the
command from table 8. Check the sign of the functions for x
0
, and that the values are approximately
zero.
8.4 Integration and differentiation
The integral, or the surface underneath a 2D function, can be determined with the command trapz.
trapz does this by measuring the surface between the x-axis and the data points, connected by
straight lines. The accuracy of this method depends on the distance between the data points:
>> x = 0:0.5:10; y=0.5*sqrt(x)+ x .* sin(x);
>> integral1 = trapz(x,y)
integral1 =
18.16550108689927
>> x = 0:0.05:10; y=0.5*sqrt(x)+ x .* sin(x);
>> integral2 = trapz(x,y)
integral2 =
18.38461257420812
A more accurate result can be obtained by using the command quad or quad1, which also
numerically evaluate an integral of the function f specified by a string.
Let again be:
( )
( ) ( )
2 2
1 1
0.1 0.1 1 0.
f x
x x
= +
+ + 1
.
>> f = '1./((x-0.1).^2 + 0.1) + 1./((x-1).^2 + 0.1)';
>> integral1 = quad(f,0,10)
>> integral1 = quadl(f,0,10)
You can also add an extra parameter to quad and quadl, specifying the accuracy.
)
Exercise 41
Find the integral of the function
2
2
( )
x
f x e
with x = 0.42:
s = 0; x = 1;
x0 =
while (x > 1e-6 * s) % relative
0.42;
s = s + x;
x
end
= x * x0;
Example
s = 0; x = 1;
x0 =
while (x > 1e-6) % absolute
0.42;
s = s + x;
x
end
= x * x0;
Example
The approximation of 1/(1-x) is returned in s (from sum). In the first case (on the left), the iteration
process is stopped as soon as the next term t is smaller than an absolute value of 1e-6. In the other
case, this value is relative to the outcome.
)
Exercise 43
Find the Taylor expansion of
1
1 x
and
3 5
sin( ) ...
3! 5!
x x
x x = + with an absolute precision of 1e-5
for x=0.4 and x=-0.4.
It is also possible to determine the Taylor expansion of
2
1
1 ...
1
x x
x
= + +