Unformatted text preview:

MATH2071: LAB #10: Iterative MethodsIntroduction Exercise 1Poisson equation matrices Exercise 2The conjugate gradient algorithm Exercise 3Compact storage Exercise 4Conjugate gradient by diagonals Exercise 5Exercise 6Exercise 7Exercise 8Exercise 9Exercise 101 IntroductionThe methods for solving matrix equations that we have seen already are called “direct” methods because theyinvolve an algorithm that is guaranteed to terminate in a predictable number of steps. There are alternativemethods, called “iterative” because they involve a repetitive algorithm that terminates when some specifiedtolerance is met. In this lab, we will focus on the conjugate gradient method applied to matrices arising fromelliptic partial differential equations. Atkinson presents the conjugate gradient method as a direct methodbecause in exact arithmetic it terminates after n steps for an n×n matrrix. In computer arithmetic, however,the termination property is lost and is irrelevant in many cases because the iterates often converge in farfewer than n steps.Iterative methods are most often applied to matrices that are “sparse” in the sense that their entries aremostly zeros. If most of the entries are zero, the matrices can be stored more efficiently than the usual formby omitting the zeros, and it is possible to take advantage of this efficient storage using iterative methodsbut not so easily with direct methods. Further, since partial differential equations often give rise to largesparse matrices, iterative methods are often used to solve systems arising from PDEs.You will see the conjugate gradient method (regarded as an iterative method), and will deploy it usingmatrices stored in their usual form and also in compact form. You will see examples of the rapid convergenceof the method and you will solve a matrix system stored in compact form that is so large that you could noteven construct and store the matrix in its usual form in central memory, let alone invert it.You may find it convenient to print the pdf version of this lab rather than the web page itself. This labwill take three sessions.2 Poisson equation matricesIn two dimensions, the Poisson equation can be written as∂2u∂x2+∂2u∂y2= ρ (1)where u is the unknown function and ρ is a given function. The simplest discretization of this equation isbased on the finite difference expression for a second derivatived2φdξ2=φ(ξ + ∆ξ) − 2φ(ξ) + φ(ξ − ∆ξ)∆ξ2+ O(∆ξ) (2)Suppose that the unit square [0, 1] × [0, 1] is broken into (N + 1)2smaller, nonoverlapping squares, eachof side h = 1/(N + 1). Each of these small squares will be called “elements” and their corner points willbe called “nodes.” The combination of all elements making up the unit square will be called a “mesh.” Anexample of a mesh is shown below.1uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuur r r r r r r r r r r rrrrrrrrrrrrrrrrrrrrrr r r r r r r r r r r rFigure 1: A sample mesh, with interior nodes indicated by larger dots and boundary nodes by smaller ones.The nodes in Figure 1 have coordinates (x, y) given byx = jh for j = 0, 1, . . . , 11y = kh for k = 0, 1, . . . , 11 (3)You will be generating the matrix associated with the Poisson equation during this lab. Employing (2)to (1) at a mesh node (x, y) yields the expressionu(x + h, y) + u(x, y + h) + u(x − h, y) + u(x, y − h) − 4u(x, y)h2= ρ(x, y).Denoting the point (x, y) as the “center” point, the point (x + h, y) as the “right” point, (x − h, y) as the“left” point, (x, y + h) as the “above” point, and (x, y − h) as the “below” point yields the expression1h2uright+1h2uleft+1h2uabove+1h2ubelow− 41h2ucenter= ρcenter. (4)The matrix equation (system of equations) can be generated from (4) by numbering the nodes in someway and forming a vector from all the nodal values of u. Then, (4) represents one row of a matrix equationΠu = ρ. It is immediate that there are at most five non-zero terms in each row of the matrix Π, no matter2what the size of N. The diagonal term is −4h2, the other four, if present, are each equal to1h2, and allremaining terms are zero. We will be constructing such a matrix in Matlab during this lab.The equation in the form (4) yields the “stencil” of the differential matrix, and is sometimes illustratedgraphically asuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuj jjjjFigure 2: A sample mesh showing interior points and indicating a five-point Poisson equation stencil.In order to construct matrices arising from the Poisson differential equation, it is necessary to give eachnode at which a solution is desired a number. We will be considering only the case of Dirichlet boundaryconditions (u = 0 at all boundary points), so we are only interested in interior nodes as illustrated in Figure1. These nodes will be counted, starting with 1 on the lower left and increasing up and to the right to 100at the upper right, as illustrated in Figure 3.Exercise 1:(a) Write a script m-file, named meshplot.m that, for N=10, generates two vectors, x(1:N^2) andy(1:N^2) so that the p oints with coordinates (x(n),y(n)) for n=1,...,N^2 are the interiornode points illustrated in Figures 1, 2 and 3 and defined in (3). You can use the following outlineN=10;h=1/(N+1);n=0;3±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯±°²¯r r r r r r r r r r r rrrrrrrrrrrrrrrrrrrrrr r r r r r r r r r r r123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100Figure 3: Node numbering% initialize x and y


View Full Document

Pitt MATH 2071 - Iterative Methods

Download Iterative Methods
Our administrator received your request to download this document. We will send you the file to your email shortly.
Loading Unlocking...
Login

Join to view Iterative Methods and access 3M+ class-specific study document.

or
We will never post anything without your permission.
Don't have an account?
Sign Up

Join to view Iterative Methods 2 2 and access 3M+ class-specific study document.

or

By creating an account you agree to our Privacy Policy and Terms Of Use

Already a member?