Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 15

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 15)

Fortunately, many banded systems arising in practice arediagonally dominant.68CHAPTER 2SYSTEMS OF LINEAR EQUATIONSTRIDIAGONAL LINEAR SYSTEMSSome special cases ofand mu occur sufficiently often that special codes are writtenfor them. For example, whenthe coefficient matrix is called tridiagonal.For the numerical solution of partial differential equations by a number of importantmethods, it is necessary to solve a great many tridiagonal systems involving a greatmany unknowns, perhaps thousands of each. This would be impractical with Factor/Solve, but with a special-purpose algorithm it is not that difficult.

Let us assumethat the tridiagonal system is written asWhen no pivoting is done, elimination zeros out the bi to lead to a structure likeAs we shall see, the ci are unchanged from the original system. To show this and toderive formulas for fi and ei , first observe that f1 = a1 and e1 = d1 since Gaussianelimination without pivoting does not change the first equation. To eliminate b2 themultiplier is m2 = b 2 /f1 . Hence,(as we stated)Notice that x1 does not appear in any other equation, so this completes the first stageof the elimination. To complete the derivation we use induction. Assume that we havederived the fi and ei through row k. Then we have the patternin rows k and k + 1. Clearly the multiplier is mk+l = bk+l / fk; then eliminate on rowk + l to getwhich finishes the elimination of variable xk.

In algorithmic form this isf1 = a1for k = l,...,n - 1 begin2.5 MATRICES WITH SPECIAL STRUCTURE69end k.For the solution of many systems with the same coefficient matrix, we save the multipliers mk as well as the fk, ck and solve separately for the data d1 , .

. . , dn. Forwardelimination ise1 = d1for k = l , . . . , n - 1 beginek+l = dk+l - mk+lekend k.Back substitution is simplyx n = e n /f nfor k = n - l, n - 2, . . . , 1 beginxk = (ek - ckxk+l) / f kend k.Storage can be managed extremely efficiently. A general matrix requires storage forn2 entries, but a tridiagonal matrix requires storage for only 3n - 2 nonzero entries.A natural scheme is to store the three diagonal bands with ak, bk, and ck as threevectors of length n.

We can Write mk over bk and fk over ak as they are computed;during the forward and backward substitution stages ek and xk can overwrite dk so thatonly one additional n vector is needed. We leave as an exercise the task of countingarithmetic operations to see that they are dramatically less than for a general matrix(see Exercise 2.15). The above algorithm assumed no pivoting, but as was seen earlierin this chapter, it is not always possible to solve linear systems without pivoting, andeven when it is possible, the numerical results may be poor. For tridiagonal systemsthere is a simple condition that guarantees that all goes well, a condition often satisfiedby systems that arise in practice.

First note that if any ck or bk vanishes, the systemcan be broken into smaller systems that are also tridiagonal. Hence, we can assumethat ck0 and bk0 for all k. The key assumption is thatThis condition is a little stronger than being diagonally dominant, enough that we canprove that the matrix is nonsingular. The argument is by induction. By assumption|m2| = |b2/a1| < 1. Supposing now that |mj| < 1 for j = 2,. . .

, k,but70CHAPTER 2SYSTEMS OF LINEAR EQUATIONSThis implies that |mk+1| < 1 as desired. From the above inequality we have |fk| >|bk+l| > 0, soThus, under these assumptions all the quantities computed in the elimination are welldefined and they are nicely bounded in terms of the data. In particular, the matrix mustbe nonsingular.SYMMETRIC MATRICESSo far we have not considered how we might take advantage of a symmetric matrix A.It might seem possible to decompose A into the product of an upper triangular matrixU and its transpose UT, which is a lower triangular matrix. However, it is not alwayspossible to decompose a symmetric matrix in this way.

This follows from the factthat such a decomposition implies that the matrix must be positive definite and not allsymmetric matrices are positive definite. To see this, if A = UTU for a nonsingularmatrix U, thenwhere y = Uv. Because U is nonsingular, if v0, then y0 andshowing that A is positive definite. Although we shall not prove it, any positive definite,symmetric A can be decomposed as A = UTU for a nonsingular upper triangular matrixU. Symmetric, positive definite matrices are very important. In applications such asthe least squares fitting of data and the variational formulation of the finite elementsolution of partial differential equations, the quantity v T Av is a kind of “energy” andis naturally positive.

We shall discuss a very effective way of solving problems withA that are symmetric, positive definite. There are more complicated ways of dealingwith symmetric matrices that are not positive definite that approach the efficiency ofthe definite case for large systems, but they do not cope with the storage nearly sowell. Codes can be found in LINPACK [5] or LAPACK [l] for both cases.

Supposingthat A is symmetric, positive definite and using the fact stated that it can be factoredas A = UTU, we can obtain U by specializing the recipe given earlier. Now we are tohave LT = U. Thussoand as beforeu1j = alj/ull for j = 2, . . . , n.2.5 MATRICES WITH SPECIAL STRUCTURE71Now(2.28)from which we findThen, as before,This decomposition has excellent numerical properties. From (2.28) we see thathencewhich says the multipliers in Gaussian elimination cannot get large relative to A.

Thisdecomposition is called the Cholesky or square root decomposition. The square rootsin this algorithm can be avoided if we use a factorization of the form LDLT, where Dis diagonal and L is lower triangular with ones on the diagonal. As in the case of theLU decomposition, when A is a band matrix, so is U. The Cholesky decompositionpreserves more than the band structure of a matrix. By examination of its recipe it isseen that as one goes down a column of U, the first (possibly) nonzero element occursin the same place as the first nonzero element of A. This says that the “profile” or“envelope” or “skyline” of A is preserved.

Obviously it is more trouble to work with adata structure that takes advantage of this fact than with one suitable for a band, but itis not much more trouble and the storage can be reduced quite a lot. Renumbering theunknowns alters the band width and the envelope of a matrix. There are algorithms thatattempt to find the best numbering in the sense of minimizing the storage and cost ofcomputing the Cholesky factorization. Many techniques have been developed for thesolution of large, symmetric, positive definite systems when most of the components ofA are zero.

The monograph [8] explains the methods and presents codes; see [9] also. Itis possible to solve efficiently systems arising in many areas of scientific computationthat involve thousands of unknowns.72CHAPTER 2SYSTEMS OF LINEAR EQUATIONSEXERCISES2.19 Count the arithmetic operations required by the algorithm in Section 2.5.2 for a linear system of n equations with a tridiagonal coefficient matrix and m right-hand sides. Compare this with what is required for ageneral system.2.6 CASE STUDY 2All the standard methods of approximating the solution of elliptic partial differentialequations, PDEs, involve the solution of a system of linear equations. The numericalsolution of PDEs is a large subject and there are many books devoted to the solutionof particular classes of problems and to the use of particular kinds of methods.

TheGalerkin method is a very important one that we illustrate here with a simple examplefrom [4]. After some preparation, the velocity w(x,y) of the steady flow of a viscousfluid in a duct with square cross section can be found as the solution ofsubject to no-slip boundary conditionsw(x,y) = 0 on |x| = 1 and |y| = 1.The same mathematical problem arises in the torsion of a bar of square cross section.Galerkin methods approximate w(x,y) by an expansionThe choice of the basis functions φj is of fundamental importance.

For the examplewe shall useNotice that each φj satisfies the boundary conditions. Also, the problem is such thatw(x,y) has certain symmetry properties, properties shared by these basis functions.Because the φj reflect well the qualitative behavior of the solution, we can hope to geta reasonable approximation with just a few terms in the expansion. This is a globalGalerkin method because each φj approximates the solution over the entire domain.We shall see in Chapter 3 that approximating functions by piecing together polynomialapproximations over subdomains can be very effective. In the present context the subdomains are called elements, and a Galerkin method based on piecewise polynomialbasis functions is called a finite element method.

In any case, when wN is substitutedinto the differential equation, there is a residualThe idea is to find coefficients aj that make this residual small in some sense. Generally there is a residual arising from boundary conditions, too, but for simplicity we2.6 CASE STUDY 273discuss only the case of basis functions that satisfy the boundary conditions exactly. Toquantify the notion of being small, we make use of an inner product of two functionsf(X,Y) and g(x,y):The Galerkin method requires that the residual of the approximation be small in thesense that it is orthogonal to all the basis functions, that is, (φi ,R) = 0 for i = 1,.

Характеристики

Список файлов книги