>> helpdeskThere you will find a whole section on MATLAB functions including the obvious ones like sin, cos, asin, acos, abs, exp, log, log10 and so on as well as matrix functions like inv, det, rref, norm, eig and many others. You may consult the help files to find out more about the built-in functions when you need them. Here we want to learn how to define our own MATLAB functions. As a simple example, we can define functions that correspond to the elementary row operations on a matrix.
Example 1: Given an m x n matrix A, a scalar k and an integer
i that represents a row index we want to define a function rowmultiple(A,k,i)
whose
output will be a matrix identical to A except that row i has been multiplied
by the scalar k.
To do this you must put the matlab commands in an m-file. (To
create such a file (assuming the path is set properly) on a PC go into
the file menu and select New, then M-file. Type in the commands for your
new function (given below), save it as "rowmultiple.m" and exit the text
editor):
The text you need in the file "rowmultiple.m" will be:
function B = rowmultiple(A,i,k)
B = A;
B(i,:) = k*A(i,:);It would be a good idea to add some explanatory comments at the beginning of the definition as in:
function B = rowmultiple(A,i,k)
% ROWMULTIPLE Given an m x n matrix A
% along with an integer i between 1 and m and a nonzero
% scalar k this function will produce an m x n matrix B
% that is the result of multiplying row i of A by k.
B = A;
B(i,:) = k*A(i,:);If this m-file is stored in your current directory, then you will be able to use it by typing commands like those given below. To review how to save your files in the correct directory see Intro To MATLAB
>> A = [ 1 2; 3 4]
>> B = [ 0 1 1; 2 1 -1; 2 3 4]
>> A = rowmultiple(A,1,5) % multiply row 1 of A by 5
>> B = rowmultiple(B,2,7) % multiply row 2 of B by 7
>> help rowmultiple % will display the help comments you put into
% the definition
>> type rowmultiple % will display the text of the file rowmultiple.mIn general it is a good idea to put some error checking and error messages into your function definitions. For example you might want to print an error message is the user tries to multiply row 5 by a scalar in a matrix that has fewer than 5 rows. However, we will not worry about such things now.
Example 2: Given a matrix A and two integers i and j we want a function that will interchange row i and row j of A. In a file called "rowswap.m" we put the following text:
function B = rowswap(A,i,j) % ROWSWAP Given an m x n matrix A and two integers % i and j between 1 and m, this function % produces an m x n matrix that results from % interchanging row i and row j of A. B = A; B(i,:) = A(j,:); B(j,:) = A(i,:);Example 3: Finally, we define a function rowcombine whose definition is stored in the file "rowcombine.m" and looks like
function B = rowcombine(A,i,j,k) % ROWCOMBINE % Makes a new matrix B by replacing % row i of A with % (row i of A) + k*(row j of A) B = A; B(i,:) = A(i,:) + k*A(j,:);Exercise 1: Define the functions rowcombine, rowswap and rowmultiple as above and use them to find the reduced row echelon form of the matrix M given below. Check your results using the built-in command rref(M).
M = [ 3 1 5 7 ]
[ 2 4 -1 0 ]
[ 5 7 3 1 ]
We could define a function that takes as its arguments an integer n,
representing the fineness of the grid and four vectors representing the
boundary temperatures. As its output the function will produce an
n x n matrix T that holds the equilibrium temperatures for each cell in
the grid. Such a function is given below:
function T = tempgrid(n,top,bottom,left,right) T = zeros(n); for i = 1:n T(1,i) = top(i); T(n,i) = bottom(i); T(i,1) = left(i); T(i,n) = right(i); end; variables = (n-2)^2; M = zeros(variables,variables + 1);
k = 0; for row = 2:n-1 for col = 2:n-1 k = k + 1; M(k,k) = 1; if row == 2 M(k,variables + 1) = M(k,variables + 1) + T(row - 1,col)/4; else M(k,k-n+2) = -1/4; end; if row == n-1 M(k,variables + 1) = M(k,variables + 1) + T(row +1,col)/4; else M(k,k+n-2) = -1/4; end; if col == 2 M(k,variables + 1) = M(k,variables + 1) + T(row,col-1)/4; else M(k,k-1) = -1/4; end; if col == n-1 M(k,variables + 1) = M(k,variables + 1) + T(row,col+1)/4; else M(k,k+1) = -1/4; end; end; end; R = rref(M); T(2:n-1,2:n-1) = reshape(R(:,variables+1),n-2,n-2)';Once you have the MATLAB commands above stored in a file called "tempgrid.m" you can use it to rework the problem we did in the first MATLAB homework as follows:
>> n = 5;
>> top = ones(1,n); bottom = zeros(1,n);
>> left = ones(1,n); right = zeros(1,n);
>> tempgrid(n,top,bottom,left,right)Exercise 2: As we let n increase we get a better and better approximation to the equilibrium temperature distribution. How large can you let n go before your computer crashes or it takes more than 5 minutes to complete the computation? Experiment.
________________________________ | | | | 1 | 1 | 1 | 1 | 1 a = (1 + 1 + b + d)/4 _____|______|______|______|_____ b = (a + 1 + c + e)/4 | | | | c = (b + 1 + 0 + f)/4 1 | a | b | c | 0 d = (1 + a + e + g)/4 _____|______|______|______|_____ e = (d + b + f + h)/4 | | | | f = (e + c + 0 + i)/4 1 | d | e | f | 0 g = (1 + d + h + 0)/4 _____|______|______|______|_____ h = (g + e + i + 0)/4 | | | | i = (h + f + 0 + 0)/4 1 | g | h | i | 0 _____|______|______|______|_____ | | | | 1 | 0 | 0 | 0 | 0 _____|______|______|______|_____And this gives rise to a matrix equation of the form
[ a ] [ 0 1/4 0 1/4 0 0 0 0 0 ] [ a ] [ 2/4 ]
[ b ] [ 1/4 0 1/4 0 1/4 0 0 0 0 ] [ b ] [ 1/4 ]
[ c ] [ 0 1/4 0 0 0 1/4 0 0 0 ] [ c ] [ 1/4 ]
[ d ] = [ 1/4 0 0 0 1/4 0 1/4 0 0 ]*[ d ] + [ 1/4 ]
[ e ] [ 0 1/4 0 1/4 0 1/4 0 1/4 0 ] [ e ] [ 0 ]
[ f ] [ 0 0 1/4 0 1/4 0 0 0 1/4] [ f ] [ 0 ]
[ g ] [ 0 0 0 1/4 0 0 0 1/4 0 ] [ g ] [ 1/4 ]
[ h ] [ 0 0 0 0 1/4 0 1/4 0 1/4] [ h ] [ 0 ]
[ i ] [ 0 0 0 0 0 1/4 0 1/4 0 ] [ i ] [ 0 ]Let t denote the vector of 9 interior temperatures and let M denote the 9 x 9 coefficient matrix of the equation above and let b denote the vector on the right that reflects how the boundary temperatures affect the interior temperatures. There is a standard method, called Jacobi iteration, that can be used to solve such a system. It works in the following way.
Step 1: Take a wild guess for the interior temperatures at equilibrium. If you don't feel too wild, just guess that all the interior temperatures are zero. Call this first guess t_0.
Step 2: Compute a new, perhaps better guess t_1 = M*t_0 + b.
Step 3: Continue in this fashion, computing a sequence of iterates t_2, t_3, ..., t_n, .... The claim is that this sequence will converge to a solution to the original system t = M*t + b.
With this point of view we are getting more information than we got simply by solving t = M*t + b. That solution was the equilibrium temperature distribution. With this approach we see how the temperature on the plate changes over time. Choosing an initial state t_0 = 0 means that we start with a plate whose temperature is everywhere zero and then we begin to heat along the boundary. As we look at the iterated temperatures t_0, t_1, t_3,... we are seeing how the heat gradually spreads through the plate. Eventually it stabilizes at the equilibrium temperature distribution we computed before. So this method is giving more information than the one we used before.
We can modify the function tempgrid defined above to do the Jacobi iteration problem. In the tempgrid function we formed the augmented matrix [M|b] for the system Mt = b. Now we want to form the (n-2)^2 x (n-2)^2 coefficient matrix M and a separate vector b. We will call the modified function iterate and store it in a file called "iterate.m". We give it one further argument, "t_0" that will be our initial "guess" for the interior temperatures.
Now we think of the unknown interior temperatures as x_1, x_2, ..., x_9 (in the 5 x 5 grid case) and for each temperature x_k we have an equation of the form
________________________________ | | | | 1 | 1 | 1 | 1 | 1 x_k = ((cell left)+(cell above)+(cell right)+(cell below))/4 _____|______|______|______|_____ | | | | where 1 | x_1 | x_2 | x_3 | 0 _____|______|______|______|_____ cell left corresponds to x_k-1 or is a boundary cell | | | | 1 | x_4 | x_5 | x_6 | 0 cell above <---> x_k-3 or is a boundary cell _____|______|______|______|_____ | | | | cell right <---> x_k+1 or is a boundary cell 1 | x_7 | x_8 | x_9 | 0 _____|______|______|______|_____ cell below <---> x_k+3 or is a boundary cell | | | | 1 | 0 | 0 | 0 | 0 _____|______|______|______|_____In general, if we have an n x n grid, then the cell above the x_k cell will correspond to the variable x_k-(n-2) and so on. Thus we should define our function iterate as below. Given the boundary temperatures and an initial guess t_0 it will return a new guess t_1 and the corresponding matrix of temperatures T:
function [t_1,T] = iterate(n,top,bottom,left,right,t_0) T = zeros(n); % starts the same way. Set up an n x n for i = 1:n % matrix to store the temperatures T(1,i) = top(i); % boundary temps are known. Below we T(n,i) = bottom(i); % compute the interior temps. T(i,1) = left(i); T(i,n) = right(i); end; variables = (n-2)^2; M = zeros(variables,variables); % Now M will be the coefficient matrix b = zeros(variables,1); % b is a column vector, % representing the contribution % of the boundary cells k = 0; for row = 2:n-1 for col = 2:n-1 k = k + 1; if row == 2 % cell above is a boundary cell b(k) = b(k) + T(row - 1, col)/4; else M(k,k-n+2) = 1/4; end; if row == n-1 % cell below is a boundary cell b(k) = b(k) + T(row +1,col)/4; else M(k,k+n-2) = 1/4; end; if col == 2 % cell to the left is a boundary cell b(k) = b(k) + T(row,col-1)/4; else M(k,k-1) = 1/4; end; if col == n-1 % cell to the right is a boundary cell b(k) = b(k) + T(row,col+1)/4; else M(k,k+1) = 1/4; end; end; end; t_1 = M*t_0 + b; T(2:n-1,2:n-1) = reshape(t_1,n-2,n-2)';
Exercise 3: Try Jacobi
iteration for the 5 x 5 temperature grid. Take the zero vector as
your initial "guess". Does your sequence converge to the same solution
you found using the function tempgrid? How many iterations
did you need to get the same solution to 4 decimal places? To use
the function iterate you could do the following:
>> n = 5;
>> top = ones(1,n); bottom = zeros(1,n);
>> left = ones(1,n); right = zeros(1,n);
>> t = zeros((n-2)^2,1); % this is t_0
>> [ t, T ] = iterate(n,top,bottom,left,right,t) % computes t_1 and first approx to the temp. array
>> [ t, T ] = iterate(n,top,bottom,left,right,t) % computes t_2 and second approximation to T
Exercise 4: Try at least two other choices of initial
guess
t_0. Did your sequence converge to the correct solution?
in how many steps?
Exercise 5: Use Jacobi iteration to solve the system with n = 15. Again use the zero vector as your initial guess. How does it compare to computing the exact solution with our function tempgrid? You will probably need to iterate many more times so you might want to use a loop as in:
>> t = zeros((n-2)^2,1);
>> for i = 1:100
[ t, T ] = iterate(n,top,bottom,left,right,t);
end
>> T % this will display temp array after 100 time intervals
>> for i = 1:100 % if we repeat the loop
[ t, T ] = iterate(n,top,bottom,left,right,t);
end
>> T % this will display temp array after 200 time intervalsExercise 6: To understand why Jacobi iteration works we can start by looking at a one variable version of the problem. For scalars m and b consider the equation
t = m*t + b whose solution is t = b/(1-m)provided that m is different from 1. What would the method of Jacobi iteration give for this problem? For which values of m and b would the sequence t_0, t_1, ... converge to the actual solution? Hint: Suppose your initial guess t_0 has error E. That is suppose that your initial guess is
t_0 = b/(1-m) + ECompute t_1. What is the error in t_1 as an estimate for the true solution? What about t_n? Under what conditions will the error tend to 0 as n increases?
Exercise 7: (continuation of exercise 6) We want to solve the matrix equation below
t = M*t + b.using Jacobi iteration and we want to know when this method will work. From an initial state t_0 we compute
t_1 = M*t_0 + b.The error in our first estimate compared to the error in our second estimate is obtained by subtracting these equations
(t - t_1) = M*(t - t_0)What will the error in t_n be? What conditions on the eigenvalues of M would guarantee that the error goes to zero as n increases?
Exercise 8: That Jacobi iteration works for this problem is a consequModify the function iterate to get a new function tempeig(n) that computes the eigenvalues and eigenvectors of the coefficient matrix M corresponding to an n x n grid. Note this does not depend on the boundary temperatures at all so you don't need to deal with them. You just need to set up the n x n coefficient matrix M and then find its eigenvalues and eigenvectors. Use the built-in MATLAB function eig that computes the eigenvalues and eigenvectors of a square matrix. An example of the use of the function eig is given below in which the eigenvalues are stored in a diagonal matrix D and the eigenvectors are stored as the columns of a matrix S. (For more information type "help eig" in the MATLAB window.)
>> a = [ 1 2; 3 4];
>> [ S, D ] = eig(a)