ELEC2103/ELE9103 Tutorial 3 – More Programming in
MATLAB and Introduction to Numerical Methods
The first part of this tutorial continues the basics of programming in MATLAB started in Tutorial 2. It is
designed to give you more practice in MATLAB programming. Mainly, the exercises go beyond those in the
previous tutorial in that you are now being asked to devise a solution and then translate this solution into
MATLAB code. Topics discussed include M-files, control flow, structures and cell arrays, and more on
graphics.
The second part of the tutorial presents an introduction to numerical methods. Matrix decompositions and
methods for solving systems of linear equations are introduced. The application of QR factorisation to find
the least-square solution to a set of over-determined linear equations is demonstrated.
3.1. Warm-up exercises
Exercise 1
a) Rewrite the following code using a while loop to avoid using the break command:
For k = 1:10
x = 50-k^2;
if x < 0
break
end
y = sqrt(x)
end
b) Write a script using a switch case block to input one angle, whose value can be 45°, -
45°, 135° or -135°, and display its quadrant (1, 2, 3 or 4).
3.2. Easy Plotting
MATLAB has a number of functions for use when data points do not have to be specified exactly.
The function fplot plots functions defined by a M-file name or function handle. Thus
fplot('sin', [0 2*pi])
or, equivalently,
fplot(@sin, [0 2*pi])
plot the sine function over the range shown. MATLAB does the job of working out numbers of data
points, etc. The function title, axis, etc can be used as with plot.
The function ezplot plot functions defined by symbolic or string expressions. Using:
ezplot('sin(x)/x', [-10 10])
it plots xx /)sin( over the range specified. It can also be used to plot an implicit function. For
example,
ezplot('x^2+(y/3)^2-1', [-2 2 -4 4])
plots the ellipse 1)3/(
22 =+ yx .
Tutorial 3 ELEC2103/ELEC9103

Page 2
Exercise 2
Plot the function ( , ) sin( )*cos( )f x y x y= for 2 2x −   and 2 2y −   .
3.3. Practice turning problems into code
Get into some good habits:
• Try to devise a solution on paper first, and then translate this solution into MATLAB code.
• Use comments to guide you as you develop your answers.
• Try out unit testing on individual code modules.
• Keep a log of your bugs and mistakes, and learn from them

Exercise 3
Honours are awarded on the basis of a “grand wam”, an average mark over all units of study,
weighted according to credit point value and according to the nominal year of each unit of study. The
formula is:
GWAM =  iiiii YCYCM
Where M is the mark, C the credit point value and Y the year (values of 1, 2, 3 or 4). Honours are
awarded as follows:
• H1 for GWAM 75,
• H2(1) for 75 >GWAM 70,
• H2(2) for 70 >GWAM 65.
Write a function which will accept a 3-column matrix containing a mark, credit point value and year
in each row and will return a value for the GWAM (it could consist of just one line).
Joe obtained the following results :

1st Year 2nd Year 3rd Year 4th Year
Mark CP Mark CP Mark CP Mark CP
66 6 68 4 76 4 74 12
54 6 77 4 70 4
70 4 69 4
What grade of Honours will he be awarded?
Answer: Second class, division 1.
Exercise 4
Legendre polynomials, ,1,0),( =nxPn , can be computed using Bonnet’s recursion formula, which
is defined as follows:
,3,2)()1()()12()( 21 =−−−= −− nxPnxxPnxnP nnn
with xxPxP == )(,1)( 10 .
Write a function that takes an integer n and returns the coefficients of )(xPn in descending order of
powers. Use only 1 for loop. Then use the function polyval to evaluate the Legendre polynomial
of degree 8 at x = 0.6.
Tutorial 3 ELEC2103/ELEC9103

Page 3
Answer: 8(0.6) = 0.2123
3.4. Numerical approximation
Series expansions are used to calculate the value of many trigonometric functions. For example, the
value of cos(x) can be expressed as a Maclaurin series:
cos = 1 −
2
2!
+
4
4!
−
6
6!
+ ⋯ = ∑ (−1)−1
2(−1)
(2( − 1))!
∞
=1

Likewise, the value of sin(x) is given by:
sin = −
3
3!
+
5
5!
−
7
7!
+ ⋯ = ∑ (−1)−1
2−1
(2 − 1)!
∞
=1

In practice, these functions’ values are approximated by truncating the corresponding series at an
appropriate level of accuracy. This approximation methods is the workhorse of many numerical
techniques, as demonstrated in the exercise below.
Exercise 5
The elliptic integral,
10,
)(sin1
)(
2
0
22
2 
−
=  k
tk
dt
kK


cannot be evaluated in terms of elementary functions. Carl Friedrick Gauss (1777-1855) devised an
algorithm to solve this integral which uses a sequence of arithmetic means }{ na and geometric
means }{ nb , where
( ) ,2,1,,2/
1,1
1111
2
00
==+=
−==
−−−− nbabbaa
kba
nnnnnn

Both sequences have the same limit g, and
g
kK
2
)( 2

= . Also, nn ba  for all n.
Write a MATLAB function that computes the elliptic integral
2
( )K k . Use a while loop to generate
the sequences and continue looping until epsba nn − , where eps is the relative floating point
accuracy value returned by the MATLAB function eps.
Write your code so that the function takes either a single number or an array of numbers as k2-values
(i.e. make it compatible with vector inputs). Your code should use only 1 loop, the while loop.
Return the answer in the same format (row or column vector) as the input.
What is the value of the elliptic integral for k2 = 0.3, 0.4, 0.5, 0.6?
Answer: 1. 7139, 1.7775, 1.8541, 1.9496

3.5. Circuit analysis applications
We will now move on to some circuit analysis applications.
Tutorial 3 ELEC2103/ELEC9103

Page 4
It may be useful to work through Example 10.5 in Moore – Solving Simultaneous Equations: An
Electrical Circuit – before moving on to Exercise 7.
Exercise 6
Electrical resistors are said to be connected “in series” if the same current passes through each and
“in parallel” if the same voltage is applied across each.
If in series, they are equivalent to a single resistor with resistance given by:
= 1 + 2 + 3 + ⋯ +
If in parallel, their equivalent resistance is given by:
=
1
1
+
1
2
+
1
3
+ ⋯ +
1


Write an m-file that prompts the user to input (i) a vector containing the individual resistors
resistances, and (ii) the type of connection (series or parallel) and then computes the equivalent
resistance.

Exercise 7

a) Derive the equations describing the circuit above, and encode them in a script for values R1 = 5,
R2 = 100, R3 = 200, R4 = 150 and R5 = 250 kΩ.
b) Suppose the resistors are rated to carry a current of no more than 1 mA. What is the allowable
range of positive values for the voltage V2?
c) Investigate how the resistance of R3 limits the allowable range for V2, by generating a plot of the
allowable limit on V2 as a function of R3 for the range 150 ≤ R3 ≤ 250 kΩ.

Optional Section – Optimisation
Engineers work to improve designs and operations. One way of doing this is called optimization,
which uses a mathematical description of a design to select the best values of certain variables. Many
sophisticated mathematical tools for large complicated optimization problems are available in the
MATLAB Optimisation Toolbox.
However, some problems can be solved in a more direct way, using exhaustive search. The next
problem is an example of distribution centre location optimization, which you might see is related to
planning problems that arise in power and telecommunications.
Tutorial 3 ELEC2103/ELEC9103

Page 5
Exercise 8 – Optional (not assessed)
A company wants to locate a distribution call centre that will serve six of its major customers in a
30x30km area. The customers’ locations are given in the following table, listing their (x,y)
coordinates, alongside the volume that must be delivered each week:

Customer x location (km) y location (km) Volume
(tonnes/week)
1 1 28 3
2 6 18 8
3 8 16 7
4 19 2 5
5 25 10 2
6 27 8 5
a) Plot the locations of each customer, and include a legend.
The weekly delivery cost ci for customer i depends on both the volume vi and distance di from the
distribution centre (assume the distances are straight-line distances).and is given by ci = 0.5divi.
b) Write a script that finds the location of the distribution centre that minimizes the total weekly
cost to service all six customers, to the nearest km, and returns the weekly cost.
c) Generate two plots of the costs for all locations: a contour plot and a mesh plot


Numerical Methods
We now consider efficient numerical methods for solving systems of linear equations and least squares
problems, using the LU and QR matrix factorisations.

3.6 Numerical Efficiency
Numerical methods are evaluated according to their computational efficiency. MATLAB provides
three methods to measure computation time: tic, cputime, clock. For example,
>> A = rand(1000,1000); B = rand(1000,1000);
>> tic; C = A*B; toc
>> t0 = cputime; C = A*B; t1 = cputime; t1-t0
>> t0 = clock; C = A*B; t1 = clock; t = etime(t0,t1)
If you are interested in benchmarking your computer, try typing bench at the MATLAB prompt.

Exercise 9
Determine the time required to multiply two 1000x1000 random matrices.

Tutorial 3 ELEC2103/ELEC9103

Page 6
Matrix Decompositions and Linear Systems of Equations
Linear systems of equations are a set of equations such as:

5 8 5
2 4 10 1
6 2 5 3 .
x y z
x y z
x y z
+ + =
− + =
+ − =

You will have solved several such equations in your Linear Algebra course using Gaussian
elimination. Although Gaussian elimination is useful for solving a matrix equation by hand, it is not
an efficient computational algorithm and is not used in numerical computations.
The best way to solve a system of linear equations such as = in MATLAB is to just type
\x b=M . This is much more efficient than computing inv(M)*b. The backslash operator solves
the equation using an LU matrix factorisation, which is the most common numerical method for
solving linear equations. Matrix decompositions refer to a factorization of matrices and are extremely
useful. We examine two matrix decompositions: (1) LU decomposition; and (2) QR decomposition.

LU Decomposition
A matrix, M , can be factored so that
 = P M L U
where L is a lower triangular matrix, U is an upper triangular matrix, and P is a permutation matrix
that permutes the rows of M (the permutation matrix results from pivoting to avoid round-off errors).
For example,
>> M = [1 5 8; 2 -4 10; 6 2 -5];
>> [L,U,P] = lu(M)
L =
1.00 0 0
0.33 1.00 0
0.1667 -1.00 1.00
U =
6.00 2.0000 -5.00
0 -4.6667 11.6667
0 0 20.5
P =
0 0 1.00
0 1.00 0
1.00 0 0
Using an LU decomposition, the matrix equation = can be re-written as () = and can
be solved easily in two steps:
(1) solve = ,
(2) solve = .


Tutorial 3 ELEC2103/ELEC9103

Page 7
For example, if Mand b are
>> M = [1 5 8; 2 -4 10; 6 2 -5];
>> b = [5 1 3] .';
then the solution to = can be determined as follows:
>> [L,U,P] = lu(M);
>> y = L\(P*b);
>> x = U\y;
The reason for using an LU decomposition is that equations (1) and (2) involve triangular matrices
which are easily solved by a sequence of substitutions. This is very efficient numerically (you are
asked to show this in Exercise 2). Importantly, L and U are reusable. Suppose we want to solve
= , = and = . In other words, suppose we have a number of linear systems with the
same matrix. Using the backslash operator three times is not the best method of solution because
MATLAB would perform an LU decomposition of the matrix M three times. Instead, it is better to
use a single LU decomposition to find the matrices L, U, P and then solve the three equations using
L, U, P as above.

Exercise 10
(a) Compare the computation time required for inv(M)*b and M\b, where
>> M = rand(1000,1000); b = [1:1000] .';
(b) From the LU decomposition example in the text or your own example, verify that L and U
are triangular matrices. Verify that P is a permutation matrix. Compare P M with M P .
What is the inverse of P? Explain why = and = are easily solved by a sequence
of substitutions.
(c) Compare the computation time required for solving = , = and = using the
backslash operator three times and using a single LU decomposition and the method outlined
above, where
>> M = rand(1000,1000);
>> a = rand(1000,1); b = rand(1000,1); c =rand(1000,1);

QR Decomposition
A matrix M can be factored so that
= M Q R
where Q is an orthogonal matrix and R is an upper triangular matrix. The fact that Q is an orthogonal
matrix means that the columns of Q are unit basis vectors. For example,
(:, ) * (:, ) 1 (for all n) and (:, ) * (:, ) 0 (m n)T Tn n n m= = Q Q Q Q . Furthermore, T Q Q I= . For
example,
>> A = [1 2 3;4 5 6; 7 8 10];
>> [Q,R] = qr(A)

Tutorial 3 ELEC2103/ELEC9103

Page 8
Q =
-0.1231 0.9045 0.4082
-0.4924 0.3015 -0.8165
-0.8616 -0.3015 0.4082
R =
-8.1240 -9.6011 -11.9399
0 0.9045 1.5076
0 0 0.4082
>> Q*Q.'
ans =
1.0000 0.0000 0.0000
0.0000 1.0000 0.0000
0.0000 0.0000 1.0000

Least-Square Solutions to Linear Problems
The QR decomposition is commonly used to find least-square solutions to linear problems such as
curve fitting. Suppose we would like to fit the curve
2/3 2
1 2 3 4 5cos(ln5 )
xC x C x C xe C x C−+ + + +
to a set of data. In this case, we cannot directly use MATLAB’s curve fitting utility via the GUI in
the plot window because it is not a polynomial function. Nonetheless, the least-square approach can
still solve this problem because we have a set of linear equations for the coefficients iC . For
example, suppose that the data is given in a set of column vectors x and y and that we want to use the
above function to map the data in x to the data in y. One first forms the matrix M
>> M = [x, x.^(2/3), x.*exp(-2*x), cos(log(5*x)), ones(size(x))].
The equation = represents the curve fitted to the data, where C is a vector of length 5 (the 5
coefficients). This is the important point: by forming the matrix M using the values for x, we obtain a
set of linear equations for the coefficients C describing the curve.
To perform the curve fitting, one finds the best least-square error solution to the equation = .
In other words, we find C that minimizes:
( ) ( )
1/ 2
T
y C y C y C −  = −   − 
  
M M M
r r rr r r
.
The expression with the double brackets is referred to as the norm of the vector. The norm is equal
to the square root of the dot product of the vector with itself. For example,
2
i
T
i
v v v v=  = 
r r r
.
In this case, the vector is = , which is the error vector between the true y and the curve
approximation MC. Thus, the C that minimizes || − || gives the coefficients for the fitted curve.

Tutorial 3 ELEC2103/ELEC9103

Page 9
QR Decomposition and Least-Square Minimization
The solution to the minimization problem using the QR decomposition method is based on the fact
that the norm of a vector is preserved when multiplied by an orthogonal matrix. That is to say, if Q is
an orthogonal matrix, then we have:
( ) ( )
1/ 2 1/ 2 1/ 2T T T Tx x x x x x x x     = = = =    
Q Q Q Q Q .
Using this property and the QR decomposition of M, we obtain:

( )


.
T
T T
T T
T
y C y C
y C
y C
y C
−  =  − 
=  −  
=  −   
=  − 
M Q M
Q Q M
Q Q Q R
Q R
r rr r
rr
rr
rr

In order to understand the advantage of the last form of the equation, we must consider the QR
decomposition of an over-determined system of equations (as occurs during curve fitting). If the
matrix M represents an over-determined set of equations then it is not a square matrix; it has more
rows than columns. For instance, the QR decomposition of a random 5x3 matrix gives:
>> M = rand(5,3)
M =
0.6813 0.4289 0.3028
0.3795 0.3046 0.5417
0.8318 0.1897 0.1509
0.5028 0.1934 0.6979
0.7095 0.6822 0.3784
>> [Q,R] = qr(M)

Q =
-0.4751 -0.1166 0.1676 0.3486 -0.7817
-0.2646 -0.2216 -0.4852 -0.7633 -0.2506
-0.5801 0.6539 0.3435 -0.2747 0.2061
-0.3506 0.2074 -0.7524 0.4647 0.2281
-0.4948 -0.6831 0.2291 0.0663 0.4813
R =
-1.4339 -0.7998 -0.8066
0 -0.4194 -0.1704
0 0 -0.5986
0 0 0
0 0 0
Note that the matrix R is upper triangular and has 0s in every row greater than three (the number of
columns of M). Therefore, when finding C that minimises || = ||, all rows greater than three
can be ignored! In other words, the over-determined system of equations can be reduced to a square
system of equations that solves the minimization problem. For exactly this reason, MATLAB has an
Tutorial 3 ELEC2103/ELEC9103

Page 10
“economy-size” QR factorization which only includes values for Q and R necessary to solve the
minimization problem:
>> [Q,R] = qr(M); % full-size QR factorization
>> [Q0,R0] = qr(M,0); % economy-size factorization.
In order to find C that minimizes || − || in a least-square sense, one uses the following
MATLAB commands to solve 0T y C −  =Q R
r rr
: =
>> [Q0,R0] = qr(M,0);
>> z0 = Q0.'*y;
>> C = R0\z0;
Try it and see.

Exercise 11
In this exercise, we try to fit the amplitude-modulated wave
cos(1* )*cos(10* / 3)y x x = +
by the sum of six sinusoidal functions:

1 2 3
4 5 6
*cos(9* ) *sin(9* ) *cos(10* )
*sin(10* ) *cos(11* ) *sin(11* ) .
C x C x C x
C x C x C x
+ + +  
+ +

Use the following commands to generate the sample data:
>> x = [-4*pi:0.01:4*pi]. ';
>> y = cos(1*x).*cos(10*x+ pi/3);
Use QR decomposition to find 1 2 6, , , C C C  that provide the best least-square fit to the data.
Answer: 1 2 3 4 5 60.25, 0.433, 0, 0, 0.25, 0.433C C C C C C= = − = = = = −