MATLAB PROJECT 4
Please read the Instructions located on the Assignments page prior
to working on the Project.
BEGIN with creating a Live Script Project4.
Note: All exercises in this project have to be completed in the Live Script using the Live
Editor. Please refer to the MATLAB video that explains how to use the Live Script:
https://www.mathworks.com/videos/using-the-live-editor-117940.html?s_tid=srchtitle
The final script has to be generated by exporting the Live Script to PDF.
Each exercise has to begin with the line
Exercise#
You should also mark down the parts such as (a), (b), (c), and etc. This makes grading easier.
Important: we use the default format short for the numbers in all exercises unless it is
specified otherwise. We do not employ format rat since it may cause problems with
running the codes and displaying matrices in the Live Script. If format long had been used
in an exercise, please make sure to return to the default format in the next exercise.
Part I. Evaluating Eigenvalues
Exercise 1 (4 points): Difficulty: Moderate
In this exercise, you will calculate the eigenvalues of a matrix by using different methods and
different functions – one of the functions, quer, you created earlier, and the other functions
are MATLAB built-in functions.
NOTE: The function quer, which was created in Project 3, can be used for approximating
eigenvalues by successive iterations, but, unfortunately, not for all matrices the process
converges. In practice, for computing eigenvalues using MATLAB application, especially
complex eigenvalues, it is advisable to use a MATLAB built-in function eig.
**Create a function
function [] = eigenval(A)
**First, your function will calculate the eigenvalues by using QR factorization. You will
include into your code for igenval that part of the function quer from Project 3, Exercise 4
which deals with the case when m=n: it will be calculating consecutive iterations for a matrix
A, also named A, until the matrix closetozeroroundoff(A-triu(A),12)becomes the zero
matrix (notice that the parameter p=12 here). If the matrix B is your last iteration A, then your
only output for this part of the code, which uses QR factorization, will be the “sorted” vector
of the eigenvalues:
q=sort(diag(B))
supplied with the message: 'the eigenvalues of A from QR factorization are:'
(output q).
2
Note: The function sort will sort out the set of eigenvalues in the ascending order.
**Second, you will calculate the eigenvalues using a MATLAB built-in function eig:
e=sort(eig(A))
and output the sorted eigenvalues with a message 'the eigenvalues of A from a MATLAB
built-in function are:'
(output e).
**Third, there is a MATLAB built-in function poly(A) that outputs the vector of coefficients
of the characteristic polynomial of a matrix A written in the descending order according to the
degree of the polynomial. Calculate the vector C=poly(A) (do not display it) and, then, use
the MATLAB function roots(C) to output the zeros of the characteristic polynomial, which
are the eigenvalues of A. You will need to “sort” them as well, that is:
c=sort(roots(C))
Supply your output with the message: 'the eigenvalues from the MATLAB
characteristic polynomial are:'
and display c.
**Then, you will check whether the vectors q, e, and c are close to each other.
Output and display the three norms:
eq=norm(e-q)
cq=norm(c-q)
qc=norm(q-c)
**Finally, we will use the vector of the eigenvalues e to build a polynomial whose zeros are
the eigenvalues of a matrix A, i.e. to build the characteristic polynomial of A. To find the
coefficients of the characteristic polynomial, we will run here a MATLAB built-in function
poly(e), which returns the vector of the coefficients of the polynomial given its zeros e.
Then, we output the polynomial itself (in a symbolic form):
R=poly2sym(poly(e))
supplied with a message that R is the characteristic polynomial of A.
**Type the functions closetozeroroundoff, jord, and eigenval in your Live Script.
Note: function jord was created in Project 1.
**Run the function eigenval on the matrices below:
(a) A=[3 3;0 3]
(b) A=jord(5,4)
(c) A=ones(5)
(d) A=tril(magic(4))
(e) A=[4 0 0 0; 1 3 0 0; 0 -1 3 0; 0 -1 5 4]
% Analyze the outputs eq, cq, qc after running the function eigenval on the matrices
(a)-(e) and write a general comment for which of the methods the vectors of the eigenvalues
are matching most closely.
3
Part II. Eigenvectors Diagonalization
Exercise 2 (6 points) Difficulty: Hard
In this exercise we, first, find all eigenvalues of an n n× matrix A. Second, we will consider
distinct eigenvalues and find an orthonormal basis for each eigenspace. Then, we will decide
whether A is diagonalizable by applying the general theory. If a matrix A is diagonalizable,
the output has to contain an invertible matrix P and a diagonal matrix D, such that, 1A PDP−= ,
or, equivalently, AP PD= and P is invertible.
**Create a function in MATLAB that begins with:
function [P,D]=eigen(A)
Your function [P,D]=eigen(A)will have a set of commands which will generate the
following outputs for an n n× matrix A. Each output has to be supplied with the
corresponding message - you could use the commands disp and fprintf.
**Find a row vector L of all eigenvalues, where each eigenvalue repeats as many times as its
multiplicity. A basic MATLAB command for this part L=eig(A) returns a column vector of
all eigenvalues of A, and we will use the function transpose to get a row vector as an
output. You will also need to sort the entries of L – use the MATLAB command sort to
output the entries of L in ascending order.
Important: In your code, you have to ensure that the multiple eigenvalues show as the same
number. You will need to use closetozeroroundoff function with the parameter p = 7 to
compare the eigenvalues and assign the same value to the ones that are within 10^(-7)
tolerance to each other.
After performing all the tasks outlined above, you will get your final output L – a row vector
of the sorted eigenvalues of A where the multiple eigenvalues all equal to each other.
**Output a row vector M of all distinct eigenvalues (no repetitions are allowed). The
MATLAB command unique applied to L can be used here.
**For each distinct eigenvalue M(i), your function has to do the following:
(1) Find the multiplicity, m(i), of the eigenvalue M(i). Output it with the message:
fprintf('Eigenvalue %d has multiplicity %i\n',M(i),m(i))
(2) Produce an orthonormal basis W for the eigenspace for lambda = ( )iM . Output it
with a corresponding message. For example:
fprintf('A basis for the eigenspace for lambda=%d is:\n',M(i))
W
Hint: An appropriate command in MATLAB that creates an orthonormal basis for the Null
space of a matrix is null.
(3) Calculate the dimension d(i) of W. Output it with the message:
fprintf('The dimension of the eigenspace for lambda=%d is %i\n',M(i),d(i))
4
After completing the tasks (1)-(3), you will proceed to determining if A is diagonalizable:
**Compare the multiplicity of each eigenvalue, m(i), with the dimension of the
corresponding eigenspace, d(i), for every i. If, for at least one of the values of i, the
dimension of the eigenspace is less than the multiplicity of the eigenvalue, output the message
'A is not diagonalizable' and terminate the program. Your outputs will be the empty
matrices: P=[];D=[];
**If the multiplicity of each eigenvalue matches the dimension of the corresponding
eigenspace, your function will continue with constructing a diagonalization:
(1) Output a matrix P, whose columns are the vectors in the bases W for the eigenspaces.
(2) Output a diagonal matrix D with the corresponding eigenvalues on the main diagonal.
(3) Verify that AP=PD and P is invertible. To verify the first condition, you should use the
function closetozeroroundoff with the parameter p = 7. To verify the second condition,
please use the command rank. If both conditions hold, output a message: 'Great! I got a
diagonalization!' Otherwise, output a message: 'Oops! I got a bug in my code!'
and terminate the program.
**If diagonalization is confirmed, your code will continue.
There is a MATLAB built-in function
[U,V]=eig(A)
which, for a diagonalizable matrix A, generates a diagonal matrix V, with the eigenvalues of
A on its main diagonal, and an n n× invertible matrix U of the corresponding eigenvectors,
such that AU=UV. Place this function in your code and display the outputs U and V.
**Compare the output matrix U with the matrix P. Write a set of commands that would check
on the following case takes place for the matrices P and U:
• P and U have either the same set of columns or the columns match up to a scalar (-1)
(the order of the columns does not matter).
If it is the case, output the corresponding message.
**Verify that the matrices D and V have the same set of diagonal elements (the order does not
count either). If it is the case, output a message 'The diagonal elements of D and V
match'. If it is not the case, output something like: 'That cannot be true!'
**Type the functions eigen, closetozeroroundoff,jord (jord was created in Project 1).
**Run the function [P,D]=eigen(A) on the following matrices:
(a)A=[3 3; 0 3]
(b)A=[4 0 1 0; 0 4 0 1; 1 0 4 0; 0 1 0 4]
(c)A=jord(5,4)
(d)A=ones(4)
(e)A=[4 1 3 1;1 4 1 3;3 1 4 1;1 3 1 4]
(f)A=[3 1 1;1 3 1;1 1 3]
(g)A=magic(4)
(h)A=magic(5)
5
Exercise 3 (3 points) Difficulty: Moderate
In this exercise we will construct an orthogonal diagonalization of an n n× symmetric matrix.
Definition: A square matrix A is called symmetric if TA A= .
Theory: A matrix A is said to be orthogonally diagonalizable if there exists an orthogonal
matrix P ( 1 TP P− = ) and a diagonal matrix D, such that 1A PDP−= , or equivalently,
TA PDP= .
An n n× matrix A is orthogonally diagonalizable if and only if A is symmetric.
**Create a function in MATLAB
function []=symmetric(A)
**First, your function has to verify whether A is symmetric. If not, output a message that the
matrix A is not symmetric and terminate the program.
**If matrix A is symmetric, you will construct an orthogonal diagonalization. You should use
a built-in MATLAB function [P,D]=eig(A) in your code. The outputs will be an orthogonal
matrix P of the eigenvectors of A and the diagonal matrix D with the corresponding
eigenvalues on the main diagonal.
**You should verify in your code that you have got an orthogonal diagonalization, that is,
both conditions hold: AP=PD and P is an orthogonal matrix. If it is the case, output a
corresponding message. (Please make sure that you will have the message displayed after
running the function.)
Note: you will need to use the function closetozeroroundoff with p = 7 in your code for
verification of each of the two conditions above.
**type the functions symmetric and closetozeroroundoff in your Live Script.
**Run the function symmetric on the following matrices:
(a) A=[2 -1 1;-1 2 -2;1 -1 2]
(b) B=A*transpose(A)
where A is a matrix in (a)
(c) A=[3 1 1;1 3 1;1 1 3]
(d) A=[5 8 -4;8 5 -4;-4 -4 -1]
(e) A=[4 3 1 1; 3 4 1 1 ; 1 1 4 3; 1 1 3 4]
Part III. Orthogonal Projections Least-Squares Solutions
Exercise 4 (5 points) Difficulty: Moderate
In this exercise, you will create a function proj(A,b)which computes the projection p of a
vector b onto the Column Space of an m n× matrix A.
Theory: When the columns of A are linearly independent, a vector p is the orthogonal
projection of a vector ∉b Col A onto the Col A if and only if ˆA=p x , where xˆ is the unique
least-squares solution of A =x b , or equivalently, the unique solution of the normal equations
T TA A A=x b .
6
Also, if p is the orthogonal projection of a vector b onto Col A, then there exists a unique
vector z, such that, it is orthogonal to Col A and = +b p z .
**Your program should allow a possibility that the columns of A are not linearly independent.
In order for the algorithm to work, we will use the function shrink to create a new matrix,
also denoted A, whose columns form a basis for Col A. Thus, we assign in our code:
A=shrink(A);
The function shrink was created in Exercise 2 of Project 3 and you should have it in your
Current folder.
You will create a function proj as defined below. The inputs are an mxn matrix A and a
column vector b. The outputs will be the projection p of the vector b onto the Col A and the
vector z which is the component of b orthogonal to the Col A.
**Your function proj will begin with:
function [p,z]=proj(A,b)
A=shrink(A);
m=size(A,1);
**First, the program checks whether the input vector b has exactly m entries, where m is the
number of rows in A. If it doesn’t, the program breaks with a message: “No solution:
dimensions of A and b disagree”, and returns the empty vectors p and z.
If b has exactly m entries, proceed to the next step:
**Determine whether it is the case that ∈b Col A. You should use a MATLAB command
rank. If ∈b Col A, you should just assign the values to the vectors p and z based on the
general theory (no computations are needed!) and output a message “b is in the Col A”. After
that, the program terminates.
If ∉b Col A, proceed to the next step:
**Determine whether it is the case that b is orthogonal to Col A. Due to the round-off errors
in MATLAB computations, you will need to use closetozeroroundoff function with p = 7
to check orthogonality. If b is orthogonal to Col A, you should just assign the values to the
vectors p and z (no computations are needed!) and output a message “b is orthogonal to the
Col A”. After that, the program terminates.
**If b is not orthogonal to Col A, proceed to the next steps:
Find the solution of the normal equations (see the Theory above), a vector x (use the inverse
matrix or the backslash operator, \ ). Output the least squares solution, x, with a message:
'the least squares solution of the system is' (output x).
**Then, calculate and output the vector
x1=A\b
and check, by using the function closetozeroroundoff with p=12, whether the vectors x and
x1 match within the given precision. If yes, output a message: 'A\b returns the least-
squares solution of an inconsistent system Ax=b'
**Next, calculate the vectors p and z (see the Theory above).
7
**Last, check whether the vector p is, indeed, orthogonal to the vector z – the function
closetozeroroundoff with p=7 should be used here. If your code confirms that ⊥p z ,
display a message: 'p and z are orthogonal. Great job!' Otherwise, the message
could be something like: 'Oops! Is there a bug in my code?'
**Type the functions closetozeroroundoff, shrink, proj in your Live Script.
**Run the function [p,z]=proj(A,b)on the following sets of matrices A and vectors b.
Notice that some of the matrices are created in several steps and intermediate outputs have
been suppressed.
(a) A=magic(4); A=A(:,1:3),b=(1:4)'
(b) A=magic(6), E=eye(6); b=E(:,6)
(c) A=magic(6), b=(1:5)'
(d) A=magic(5), b = rand(5,1)
(e) A=ones(4); A(:)=1:16, b=[1;0;1;0]
(f) B=ones(4); B(:)=1:16; A=null(B,'r'), b=ones(4,1)
Bonus! (2 points)
% Analyze the outputs in parts (d) and write a comment that would explain a reason why a
random vector b belongs to the Col A.
% Based on the input matrix A and the outputs for part (f), state to which subspace the vector
b has to belong.
**You will need to run some commands in your Live Script to support your
comments/statements.
Part IV. Application to Polynomials
Applications to Linear Models - Regression Lines
In this part of the project, you will write a code that outputs the least-squares fit equation for
the given data. An equation is determined by the parameter vector c (that has to be
calculated), the design matrix X, and the observation vector y, such that, the 2-norm of the
residual vector e = y – Xc is the minimum.
Assume we are given the data points ( ),i ix y , where 1:i m= . A choice of an equation of a
curve that gives the least-squares fit to the data is, in general, arbitrary. First, we will take a
look at the equation of a straight line 1 0y c x c= + that minimizes the sum of squares of the
residuals – it is called the least-squares regression line. The parameter vector c can be
determined by computing the least-squares solution of X =c y , where
1
2
1
1
1m
x
x
X
x
=
, 1
0
c
c
=
c ,
1
2
m
y
y
y
=
y
.
8
Note: the MATLAB command c=lscov(X,y)calculates the least-squares solution of the
equation X =c y in the same way as we would have found it by solving the normal equations,
( )T TX X X=c y ,
by using either the backslash operator, that is, ( ) ( )\T TX X X=c y , or the inverse matrix.
The data vectors x and y will be represented by the row vectors – we use the transpose
function in the code below to convert them into the column vectors. The parameter vector c,
which is used for plotting with polyplot, has to be a row vector – in the code below it
appears as c'.
Exercise 5 (3 points) Difficulty: Easy
**Create the following two functions in MATLAB:
function []=polyplot(a,b,p)
x=(a:(b-a)/50:b)';
y=polyval(p,x);
plot(x,y);
end
and
function c=lstsqline(x,y)
x=x';
y=y';
a=x(1);
m=length(x);
b=x(m);
X=[x,ones(m,1)];
c=lscov(X,y);
% the same result will be obtained if c=inv(X'*X)*(X'*y) or c=(X'*X)\(X'*y)
% and we are verifying it by calculating:
c1=inv(X'*X)*(X'*y)
c2=(X'*X)\(X'*y)
% the next command calculates the 2-norm of the residual vector
N=norm(y-X*c)
% plot data points and the least-squares regression line:
plot(x,y,'*'),hold
polyplot(a,b,c');
% output the polynomial:
P=poly2sym(c)
end
Note: you may skip the comment lines in the body of lstsqline while creating the function –
they are placed there for your convenience.
**Type
format
**type the functions polyplot and lstsqline in your Live Script.
**Run the function c=lstsqline(x,y)
on the vectors
x = [1,3,4,5,6], y = [1,3,2,3,5]
9
The outputs have to be: the vectors c, c1, and c2; the 2-norm of the residual vector –
number N, the regression line in the form of a polynomial P in x of degree n = 1, and the plot
containing both the data points and the graph of the regression line.
Important: Please check the plot before submitting the exercise – you may need to Run
Section twice to have both the data points and the graph of the line displayed.
Exercise 6 (4 points) Difficulty: Moderate
**Create a new function in the file
function c=lstsqpoly(x,y,n)
which is a modification of the function lstsqline in Exercise 5 in the way that it outputs the
least-squares polynomial of degree n: you will need to generate a new design matrix X within
your code – it’s form depends on the degree n of the polynomial. For example, for 3n = , the
design matrix X will have a form:
3 2
1 1 1
3 2
2 2 2
3 2
1
1
1m m m
x x x
x x x
X
x x x
=
.
You will also need to have the function polyplot in your Current Folder in MATLAB - it
will be used inside the function lstsqpoly (see Exercise 5 of this project).
**type the functions polyplot and lstsqpoly in the Live Script.
**Run the function c=lstsqpoly(x,y,n) on the vectors from Exercise 5:
x = [1,3,4,5,6], y = [1,3,2,3,5]
for each n = 1, 2, 3, 4.
Your output for each n has to contain: the vectors c, c1 and c2 (see Exercise 5), the 2-
norm of the residual vector – number N, the least-squares polynomial P in x of degree n, and
the plot that contains both the data points and the graph of the polynomial.
Hint: to have both the polynomial and the data points displayed, you may need to insert a
Section Break before proceeding to the next n and Run each Section twice.
Please make sure that for each n the vectors of coefficients c, c1, and c2 are displayed and
match. If they don’t, make corrections in your code.
% Verify that your code for this exercise for n = 1 is consistent with the code for Exercise 5,
that is, the outputs of the functions lstsqline(x,y) and lstsqpoly(x,y,1) match. Write a
comment about it.
BONUS! (1 point)
% Analyze the outputs for n=4. Justify the fact that the least-squares polynomial of degree 4,
in practice, interpolates the 5 data points.
10
Part V. Applications to Differential Equations
In this part, you will work with the applications of eigenvalues and eigenvectors to the
continuous analogues of the difference equations, specifically, to the dynamical systems
involving differential equations.
Exercise 7 (5 points) Difficulty: Very Hard
Theory: Suppose 1 2( , ,..., )nx x x=x is a vector of differentiable functions of a single variable t
with the vector of its derivatives 1 2( , ,..., )nx x x′ ′ ′′ =x . Suppose also that A is an n n× matrix
and x0 is a (numerical) 1n× vector.
The initial value problem for the system of differential equations defined by the matrix A and
the initial vector x0 can be written in the form:
( ) ( )t A t′ =x x , x(0) = x0 (1)
From the theory it follows that, if a matrix A is diagonalizable, that is, A has n linearly
independent eigenvectors, then the general solution of the equation ( ) ( )t A t′ =x x has a form
( ) 1 21 1 2 2 ... ntt t n nt c e c e c eλλ λ= + + +x v v v (2)
where 1 2, ,..., nλ λ λ are real eigenvalues of A, the vectors 1 2, ,..., nv v v are the corresponding
eigenvectors (which form a basis for n ), and 1 21 2, ,..., n
tt t
ne e e
λλ λv v v are the eigenfunctions
for A, which form a basis for the solution set of ( ) ( )t A t′ =x x . Moreover, for any initial
condition x(0) = x0, there exists a unique set of weights ( )1 2., ,..., nC c c c= that defines the
unique solution (2) of the initial value problem (1).
For the case when n=2 and a diagonalaizable matrix A does not have a zero eigenvalue, the
origin can be either attractor, or repeller, or a saddle point for the dynamical system. From the
form of the general solution for n=2,
( ) 1 21 1 2 2t tt c e c eλ λ= +x v v ,
it follows that:
(1) the origin is a repeller if both eigenvalues are positive, and the direction of the greatest
repulsion is defined by an eigenvector corresponding to the largest eigenvalue;
(2) the origin is an attractor if both eigenvalues are negative, and the direction of the greatest
attraction is defined by an eigenvector corresponding the largest by absolute value eigenvalue;
(3) the origin is a saddle point if one eigenvalue is positive and the other eigenvalue is
negative; the direction of the greatest attraction is defined by an eigenvector corresponding to
the negative eigenvalue and the direction of the greatest repulsion is defined by an eigenvector
corresponding to the positive eigenvalue.
(For more details, please read Section 5.7 of the Textbook.)
**Create a function in MATLAB that begins with the following lines:
function []=invalprob(A,x0)
[~,n]=size(A);
It takes as inputs an n n× matrix A and an 1n× initial vector x0.
11
You will begin writing the code for invalprob by modify the code for the function eigen,
which is created in Exercise 2 of this project. The details on how the function eigen should be
modified are listed below – please remove or suppressed the tasks and outputs for the function
eigen that are not on the list. You will also need to replace the command null( ) with
null( ,'r') to create a “rational” eigenvector basis for n . After that, you will add some
more commands that are also listed below.
Your function invalprob will do the following:
**Check if A is diagonalizable.
If A is not diagonalaizable, output the message:
disp('A is not diagonalizable')
and terminate the program.
**If A is diagonalizable, output a message:
fprintf('There exists an eigenvector basis for R^%i\n',n)
and the function will continue with the following tasks:
(1) Output a row vector L of the sorted eigenvalues of A, where the multiple eigenvalues all
equal to each other (the entries of L have to be written in the ascending order and each
eigenvalue repeats as many times as its multiplicity).
(2) Output the matrix P of the corresponding eigenvectors 1 2[ ... ]nP = v v v with a message
that P forms an eigenvector basis for R^n (n has to be specified in your message).
(3) Output the matrix 1 21 2[ , ,..., ]n
tt t
nE e e e
λλ λ= v v v , whose columns form an eigenfunction basis
for the solution set – supply it with the corresponding message.
(Hint: type syms t in your code in order to define a symbolic variable t).
(4) Find and display the unique row vector of weights ( )1 2., ,..., nC c c c= , which you will
obtain by solving the linear system x(0) = x0. Output it with the message that the entries of the
vector C are the weights in a representation of the solution x through the eigenfunction basis.
**For the case when 2n = and A does not have a zero eigenvalue (please use the rank
command to check the second condition), your code has to continue with the following tasks:
(5) Determine whether the origin is an attractor, or a repeller, or a saddle point of the
dynamical system and display the corresponding message for each of the three cases. Also
indicate the direction of the greatest attraction (for an attractor), the direction of the greatest
repulsion (for a repeller), and both the direction of the greatest attraction and the greatest
repulsion (for a saddle point) - each of the outputs has to be supplied with a corresponding
message.
Important Note: please read the Theory above for the definitions of an attractor, a repeller,
and a saddle point – they are different from the ones for a discrete dynamical system.
**Type the function invalprob
**Run the function invalprob(A,x0) on the following pairs A,x0:
(a) A=ones(2), x0=[1;2]
(b) A=[-2 -5;1 4], x0=[2;3]
(c) A=[7 -1;3 3], x0=[-2;1]
(d) A=[-1.5 .5; 1 -1], x0=[5;4]
(e) A=[3 3;0 3], x0=[1;3]
(f) A=diag([1,2,2,3,3,3]), x0=ones(6,1)
YOU ARE DONE WITH THE PROJECTS!!!