School of Mathematics and Statistics MAST30028 Numerical Methods & Scientific Computing 2023
Assignment 2: Simulation & errors
Due:15:15 p.m. Oct 5
Note: This assignment is worth 20% of the total assessment in MAST30028. When you submit the assignment, you should submit two files: one pdf file and one zip file. The pdf file contains the answer you write, the numerical results you generated including figures or tables (0 mark if numerical results are not included in the pdf file), the comment or the explanation of the results, and the printout of your code (you may use Matlab publish or Matlab live scripts). The zip file should only contains all the m-files.
To aid marking, please start each question on a new page in your write-up.
1 Root-finding [10 marks]
MATLAB has only one built-in function for finding roots : fzero. (The Optimization Toolbox has more). To find out about it, start up Matlab and type
help fzero
This hybrid algorithm is due to Richard Brent (1946-), emeritus Professor of Computer Science, ANU.
A simplified version called fzerotx can be found in the Asst2 folder. A description of the method as well as the code is given in Sections 4.6 & 4.8 of the textbook. You can use fzerogui from the Asst2 folder to see how the algorithm works.
a. Modify fzerotx to return, in addition to the computed root, the number of iterations taken.
b. Try your modified fzerotx on the following problems and see how it performs in terms of reliability and
number of iterations.
• Find the first 2 positive solutions of the equation
sin x = cos(2x2)
• try it on the function f(x) = 1/(x − π) on the interval [0,5] and comment.
• try it on the function f (x) = 1 − (1 + 3x) exp(−3x) on the interval [-1,1] and comment.
c. Exactly what kind of stopping criterion does fzerotx use?
d. (Tricky) Answer Q4.7 in the textbook (Numerical Computing with MATLAB by Moler).
2 Backward stability vs. accuracy
In this question, we compare 3 methods of solving a linear system Ax = b
a. LU factorization with partial pivoting, using \. [ Don’t use lu explicitly, just use \] b.
x = A−1b
using inv(A)
c. the QR factorization. Any real square nonsingular matrix can be factorized into two square matrices
A = QR
where Q is a real orthogonal matrix (satisfies QT Q = I) and R is an upper triangular matrix. Using this
factorization, we get
so that
QRx = b
QT QRx = Rx = QT b
a triangular system that can be solved with \. See doc qr for how to obtain the matrices Q and R.
To test these methods, use the following test matrices:
• gallery(’randsvd’, n, condno) for n = 25, 50 and κ = 1, 104 , 108 , 1012 • gfpp(n) for n = 25, 50
and construct the vector b so that you know the true value of x i.e. choose a true solution x (not 0) and let b = Ax. I suggest you use a (fixed) random vector for your true solution, not a vector of integers (which can give misleading results). gfpp.m can be found in the Asst2 folder.
Record both the normwise relative forward error
and the relative residual η
||xˆ − x||/||x|| η = | | b − A xˆ | |
||A||||xˆ||
[13 marks]
and output your results in a table.
Explain what your results reveal about the 3 methods for solving a linear system in terms of backward
stability and accuracy.
3 Application: Nonlinear systems of equations [6+2+9 =17 marks]
You have met Newton’s method as a method for solving nonlinear equations in the form f(x) = 0. xn+1 = xn − f(xn)
f′(xn)
which gives iterates {xn} that converge rapidly to the solution x∗ provided the initial guess x0 is sufficiently
close to x∗. It can be written in the alternative form: solve in turn f′(xn)sn = f(xn); xn+1 = xn − sn
For a system of two nonlinear equations
f1(x1, x2) = 0; f2(x1, x2) = 0
or in vector notation
Newton’s method generalizes in the following way: solve the 2 × 2 linear system
then update
where J is the Jacobian matrix formed from f
f(x) = 0 J(xn)sn = f(xn)
xn+1 = xn − sn
���∂f1 ∂f1��� ∂x1 ∂x2
J=
∂x1 ∂x2
∂f2 ∂f2
a. Write a MATLAB function to implement Newton’s method for a system of 2 nonlinear equations. It should have the calling sequence
[roots, count, resids, history] = ass2Q3(func,x0,tol)
where:
• the first input argument is a function handle to a function that returns the value of the function f and the Jacobian matrix J at the position x
e.g. function [ f J] = example2( xx )
The function example2 needs to be written by the user, including the expressions for f and J. There is no need for J to be found symbolically from f.
• x0 is the initial guess ( a column vector)
• tol is the absolute tolerance. The iterations should stop once the residual ||f ||, measured using the
vector ∞-norm, is below tol.
• the first output argument is a column vector containing the solution x
• count is the number of iterations required
• resids is a vector containing the norms of the residuals at each iteration
• history is a matrix whose columns are the iterates xn i.e. each column is a 2 × 1 array.
The function should:
• use a default tolerance of 10−10 if none is given
• use a while loop but loop no more than 50 times.
Although I only need it to work for systems of 2 equations, you should be able to write it so that it can
be used (in principle) on a system of any size. My code is about 20 lines long.
b. Write a driver function to test your code on for
x21 −x2 −1=0;−x1 +x2 −1=0
with initial condition x1 = 1 and x2 = 2. The driver should produce the following output
After 5 iterations, the roots are [1.6180340 1.6180340]
c.
Modify your driver to find the solution to
x21 +x1x32 −9=0;3x21x2 −x32 −4=0
Show results for 4 different initial guesses: x0 = (1.2, 2.5), (−2, 2.5), (−1.2, −2.5), (2, −2.5).
For each case, plot the convergence history - the residual norms at each iteration - on a suitable plot and
(on a different plot) the solution trajectory - the set of iterates xn. Comment on your results.
Least squares [10 marks]
The expression z = ax2 + bxy + cy2 + dx + ey + f is known as a quadratic form. The set of points (x,y) where z = 0 is a conic section. It can be an ellipse, a parabola, or a hyperbola, depending on the sign of the discriminant b2 − 4ac. Circles and lines are special cases. The equation z = 0 can be normalized by dividing the quadratic form by any nonzero coefficient. For example, by setting the coefficient of x2 to −1, we get the normalized equation
x2 =by2 +cxy+dx+ey+f
You can use the Matlab meshgrid and contour functions to plot conic sections:
• Use meshgrid to create arrays X and Y.
• Evaluate the quadratic form to produce Z.
• Then use contour to plot the set of points where Z is zero.
[X,Y] = meshgrid(xmin:deltax:xmax,ymin:deltay:ymax);
Z = a*X.^2 + b*X.*Y + c*Y.^2 + d*X + e*Y + f;
contour(X,Y,Z,[0 0])
A planet follows an elliptical orbit. Here are ten observations of its position in the (x, y) plane:
x = [1.02 .95 .87 .77 .67 .56 .44 .30 .16 .01]’;
y = [0.39 .32 .27 .22 .18 .15 .13 .12 .13 .15]’;
a. Determine the coefficients in the normalized equation that fits this data in the least squares sense by solving a 10-by-5 overdetermined system of linear equations for the five coefficients. Plot the orbit with x on the x-axis and y on the y-axis. Superimpose the ten data points on the plot.
b. Perturb the data slightly by adding to each coordinate of each data point a random number uniformly distributed in the interval [-.005,.005]. Compute the new coefficients resulting from the perturbed data. Plot the new orbit on the same plot with the old orbit. Comment on your comparison of the sets of coefficients and the orbits. Explain your findings by examining the singular values of your 10-by-5 matrix.
