MATH0033代写、Matlab语言编程代写

日期：2022-11-29 10:21

MATH0033 Numerical Methods,X

Computational homework 2

Nonlinear equations

Both exercises 1 and 2 (marked *) to be submitted.

At least one of these will be marked.

Deadline: 23:59hrs Sunday 27th November

It is recommended that you use some of the Matlab .m files that are provided on Moodle.

Please submit your solutions via Crowdmark. You should submit a single pdf file, created

using the Matlab publish command, formatted as per the template provided on the Moodle

page. Your file should show the code you used, its output (including figures) and your

comments. Please do not include the code of any of the mfiles I gave you. Note that before

submitting to Crowdmark you should “flatten” your pdf file by opening it in a pdf viewer

(like Adobe Reader) and printing it to a new pdf using the Print dialogue (see instructions

on Moodle). Otherwise Crowdmark may not display it properly.

EXERCISE 1(*) Consider the linear system Ax = b, where, for ∈ [0, 1] and n ≥ 5, the

pentadiagonal n× n matrix A and the n-vector b are defined by

.

For a given size n and a given value epsi of ?, you can construct the matrix A? and the vector

b? using the command [A,b] = matrix(n,epsi) (code available on the Moodle page).

(a) We know from lectures that if A? is strictly diagonal dominant then the Jacobi and

Gauss-Seidel methods both converge. For which ? does this hold? (Assume that n ≥ 5.)

(b) Compute the solution of the system in the case n=5 with ? = 0.3, using both the Jacobi

and Gauss-Seidel methods. To do this, use the command itermeth (code available on

the Moodle page) and choose first ’J’ and then ’G’ as the method. Use the tolerance

tol = 10?10, maximum number of iterations 103, and initial guess x(0) = (0, 0, 0, 0, 0)T .

Report the number of iterations necessary for convergence in both cases.

(c) Still using n=5, use the built-in Matlab command eig to compute the spectral radius of

the iteration matrices of the Jacobi method BJ and the Gauss-Seidel method BGS(?)

for  = 0, 0.01, 0.02, . . . , 1. Recall that

where D is the diagonal part of A, L is the lower triangular part of A? (excluding

the diagonal) and U is the upper triangular part of A? (excluding the diagonal). Plot

the spectral radius against the value of ? for both methods on the same axes (the hold

command might be useful).

1

For which (approximate) range of ? do you expect the two methods to converge? How

does your answer compare to the answer you obtained in part (a)?

When both methods converge, which do you expect to be faster?

Which method would you recommend using if n=5 and ? = 0.5? Run both methods for

this case and check that the outcomes are consistent with your recommendation.

EXERCISE 2(*) Consider the boundary value problem: find y : [0, 1] 7→ R such that

= sin(π x) in I := [0, 1], y(0) = y(1) = 0, (1)

the exact solution of which is y(x) = π?2 sin(π x).

In lectures we are soon going to discuss how the solution y can be approximated using a finite

difference method on a set of nodes {xi}Ni=0 defined by xn = nh, n = 0, . . . , N , where h = 1/N

is the step size. We seek a set of N + 1 numbers {ui}Ni=0 approximating the values of the

solution y at the nodes {xi}Ni=0. To impose the boundary conditions we set u0 = uN = 0.

To approximate the ODE we replace the second derivative operator by the finite difference

operator

D2un :=

un+1 2un + un1

h2

.

We find that the numerical approximation u = (u1, . . . , uN?1)T to the interior values satisfies

the linear equation Au = b, for a certain matrix A and vector b. For a given N ∈ N these

may be constructed in Matlab using the script

>> h=1/N;

>> A= (2/h^2)*diag(ones(N-1,1)) - (1/h^2)*diag(ones(N-2,1),1)...

- (1/h^2)*diag(ones(N-2,1),-1);

>> b = transpose(sin(pi*h*(1:N-1)));

(If you are interested in knowing more about the background theory, then the section of the

notes where this material is discussed is §5.9. But you don’t need to look at §5.9 to answer

this question - everything you need is given to you here.)

(a) Use the built-in Matlab command eig to compute the spectral radii of the iteration

matrices BJ and BGS of the Jacobi and Gauss-Seidel methods applied to the system

above, for each N ∈ {5, 10, 20, 40, 80}. Report the results and store them in two vectors

rhoBJ and rhoBGS.

Hint: the matrices D, L, U needed to construct the iteration matrices for the two methods

can be constructed using the commands

>> D = diag(diag(A));

>> L = tril(A)-D;

>> U = triu(A)-D;

Do you predict that the two iterative methods should be convergent for these values of

N? If so, which should converge faster?

What seems to be happening to ρ(BJ) and ρ(BGS) as the system size N increases? As a

result, what do you expect to happen to the performance of the Jacobi and Gauss-Seidel

methods?

To investigate this behaviour further, present the plots produced by the commands

2

>> loglog(Nvec,1-rhoBJ)

>> hold

>> loglog(Nvec,1-rhoBG)

where Nvec=[5,10,20,40,80] is the vector of matrix sizes.

Use your results to propose an approximate relationship between the spectral radii and the

size of the matrix N that you conjecture to hold in the limit as N →∞. Give estimated

values for any constants that appear in your relationships. (I.e. if you conjecture that

ρ(B) ≈ 1? CαN or ρ(B) ≈ 1? CNα then give numerical estimates for C and α.)

Assuming that the error in both methods is approximately proportional to ρ(B)k, and

using the relationship you derive above, and can you predict (roughly) how the number

of iterations ktol required to achieve a fixed solution accuracy tol will grow as N →∞, for

both methods? Comment on the implications for the practicality (or otherwise) of using

the Jacobi and Gauss-Seidel methods (compared to e.g. direct methods) for this problem

when N is large.

(b) Solve the system forN = 5, 10, 20, 40, 80 using the Gauss-Seidel method (use the command

itermeth as in Exercise 1, with a tolerance 10?10 and an increased maximum iteration

number 105), starting from an appropriate initial guess. Plot the resulting solutions on

the same axes (don’t forget u0 and uN ), along with a plot of the exact solution y(x)

(formula stated at the start of this exercise), evaluated on the finest mesh (N = 80). Do

the solutions appear to converge as N increases?

For each N compute the error

e(N) = max

n=0,...,N

|un ? y(xn)|, (2)

and store your results in a vector error vect.

In lectures we will show that e(N) ≤ CN?2 as N → ∞ (equivalently, e(N) ≤ Ch2 as

h = 1/N → 0), where the constant C is proportional to maxx∈[0,1] |y′′′′(x)|. (This comes

from the error in replacing d2/dx2 by D2 — recall theoretical exercise sheet 1.) Do your

numerical results support or contradict this theoretical estimate?

Hint: you may find it useful to plot your results using the command

>> loglog(hvect, error_vect)

where hvect=1./Nvec is a vector of step sizes. Recall that if E(h) = Chp then a loglog

plot of E against h will produce a straight line with slope p and offset logC. Once you

have determined appropriate values of p and C you can check them by adding an extra

plot to your figure using hold and the command

>> loglog(hvect, C*hvect.^p)

(c) Now consider the equation (1), but with the right hand side sin(πx) replaced by 1. The

exact solution in this case is y(x) = (1/2)x(1 ? x). Once again, solve the system for

N = 5, 10, 20, 40, 80 (you will need to modify the vector b in your code), plot the numerical

solutions and the exact solution on the same axes, and compute the corresponding errors

(2). You might initially be surprised by what you see when you plot

>> loglog(hvect, error_vect)

Can you explain what is going on here?