Blog

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 m

ultiply two 1000×1000 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 = Uy; 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 Mb, 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 5×3 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 = R0z0; 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= = = = = =