Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 17
Текст из файла (страница 17)
A little manipulationof the equations for x and x(k+1) shows that error e(k) = x - x(k) satisfiesThis then implies thatIf the numberthis inequality implies that the process converges for any starting guess x (0) . The error decreases by a factor of ρ at each iteration,so if M is close to A in the sense that ρ is small, the process converges quickly. Noticethat the crucial quantity ρ is a kind of relative error of the approximation of A by M.Sharp convergence results are known, but they are beyond the scope of this book; see,for example, [2] or [10].The situation with iterative refinement is special because we want a very accuratevalue for x.
For reasons taken up in Chapter 1, it would be better then to compute theThisiterate by means of the difference of successive iterates,leads to the form suggested in Section 2.3.2,Although we have not provided all the details, it is not surprising that even whenA is extraordinarily ill-conditioned, M is sufficiently close to A that the result justestablished guarantees convergence.The Jacobi iteration is more easily understood by looking at components.
Equation j isA little manipulation results inWe cannot expect this iteration to work unless the entries off the diagonal of A are smallcompared to those on the diagonal. To obtain a sufficient condition for convergence,let us suppose that A is strictly diagonally dominant by rows, meaning that there is anumber ρ < 1 such that for all j,78CHAPTER 2SYSTEMS OF LINEAR EQUATIONSA little calculation shows that the errorThis holds for each j, sotelling us that the worst error is reduced by a factor of ρ at each iteration. A (possibly)large system of linear equations arises in Chapter 3 when fitting data by smooth cubicsplines.
Equation j has the formand the first and last equations have the form2hlxl + hl = blhn-lxn-l + 2hnxn = bn.The quantities hj appearing in these equations are all positive. Evidently this system ofequations is strictly diagonally dominant by rows with ρ = l/2. The Jacobi iterationwould be a reasonable way to solve these systems, but A is a symmetric tridiagonalmatrix and the special version of elimination from Section 2.5.2 is so effective thatthere is little point to iterative methods. In components the Gauss-Seidel iteration isA modification of the proof given for the Jacobi iteration shows that this iteration alsoconverges whenever A is strictly diagonally dominant by rows. The finite differenceequations of the example do not satisfy the convergence criterion just developed, butthey do have a kind of diagonal dominance for which it is possible to prove convergence of both the Jacobi and Gauss-Seidel methods.The Jacobi iteration and the Gauss-Seidel iteration are too slow to be practical formost problems arising in practice.
However, they are still important as preconditionersfor more elaborate iterative procedures. For more information about iterative methods,see [2] or [10]. The latter provides some substantial applications to the numericalsolution of PDEs.REFERENCES1. E. Anderson, et al., LAPACK User’s Guide, SIAM, Philadelphia, 1992.2. O. Axelsson, Iterative Solution Methods, Cambridge University Press, New York, 1994.3. A. Cline, C. Moler, G. Stewart, and J.
Wilkinson, “An estimate for the condition number of amatrix,” SIAM J. Numer. Anal., 16 (1979), pp. 368-375.4. C. A. J. Fletcher, Computational Galerkin Methods, Springer-Verlag, New York, 1984.5. J, Dongarra, J. Bunch, C. Moler, and G. Stewart, LINPACK User’s Guide, SIAM, Philadelphia,1979.REFERENCES 796. G.
Forsythe and C. Moler, Computer Solutions of Linear Algebraic Systems, Prentice Hall, Englewood Cliffs, N.J., 1969.7. K. Gallivan and R. Plemmons, “Parallel algorithms for dense linear algebra computations,”SIAM Review, 32 (1990), pp. 54-135.8. A. George and J. Liu, Computer Solution of Large Sparse Positive Definite Systems, PrenticeHall, Englewood Cliffs, N.J., 1981.9.
G. Golub and C. Van Loan, Matrix Computations, 2nd ed., The Johns Hopkins University Press,Baltimore, M.D., 1989.10. L. A. Hageman and D. M. Young, Applied Iterative Methods, Academic Press, Orlando, Fla.,1981.11. N. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, Philadelphia, 1996.12. B. Noble, Applications of Undergraduate Mathematics in Engineering, Macmillan, New York,1967.13. B. Noble and J. Daniel, Applied Linear Algebra, 3rd ed., Prentice Hall, Englewood Cliffs, N.J.,1988.14. R. Skeel, “Iterative refinement implies numerical stability for Gaussian elimination,” Math.Comp., 35 (1980), pp. 817-832.15.
J. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press, Oxford, England,1988.MISCELLANEOUS EXERCISES FOR CHAPTER 22.20 A finite element analysis of a certain load bearing frame yields the stiffness equationswhere α = 482,317., β = 2,196.05, and γ = 6,708.43. Here x1, x2, x3 are the lateral and x4, x5, x6 the rotational (three-dimensional) displacements correspondingto the applied force (the right-hand side).(a) Solve for x.(b) How reliable is the computation? First assume exact data, then5 × 10-7.2.21 Wang (Matrix Methods of Structural Analysis, International Textbook Company,Scranton, Pa., 1966) considers a statically indeterminate pin-jointed truss. Withthis problem is associated a statics matrix A that defines the configuration of theframework, a member stiffness matrix S that relates the elastic properties of theconstituent members, and an external force vector p that describes the appliedforces at the joints.
A displacement vector bfx that accounts for the displacement80CHAPTER 2SYSTEMS OF LINEAR EQUATIONSat each degree of freedom and an internal force vector f acting on each membersatisfiesFor one exampleThe matrix S has all zero entries except along the diagonal where the entries are{4800: 10000, 4800, 10000, 10000, 10000, 3000, 4800, 4800, 3000}.Write a program to form matrix products and determine the elements of K.
Solvefor x using the three p vectorsFind the corresponding vectors f.2.22 This exercise assumes a familiarity with matrix multiplication. An appropriateorganization of Gaussian elimination as in Factor/Solve makes it efficient to solvesystems of equations with different right-hand sides but the same coefficient matrix A. It is more complicated to deal with changes in A, but a formula called theSherman-Morrison formula makes it possible to deal with certain modificationsefficiently. Assume that we have already factored A into LU and we want to solve( A + uvT)x = b for given column vectors u, v, and b.
Show that this can be doneby first solving AZ = u and Ay = b, then formingA proper choice of u and v handles the change of one row or one column of A.For example, if row i of A is to be changed by adding to it a given row vector vT,just take the column vector u to be zero in all entries except the ith, which is 1.(a) How do you change column j in A so as to add a given column vector u ?(b) How do you choose u and v in order to change aij into aij + δ?MISCELLANEOUS EXERCISES 81(c) A change in A may make it singular.
How would this be revealed when usingthe Sherman-Morrison formula? Note that this approach may not be an accurateway to solve for x when A is poorly conditioned because the elimination is doneon A. Still for small changes to A, the accuracy should be acceptable and this isan inexpensive way to study the effect of changes to the data of A.2.23 Occasionally it is desirable to compute the determinant of a matrix A with n rowsand columns.
Using the factorization PA = LU discussed in Section 2.2, it can beshown thatdetA = (-1) number of row interchanges × product of pivots.In terms of the output from Factor (FORTRAN version) this isdetA = PVTIDX(n) * A(l,1) *···* A( n,n).In the C or C++ versions this would bedetA = pivot_index[n - l] * a[0] [0] * ··· * a[n - l] [n - l].Use this formula to compute the determinant of(a) A in Exercise 2.13 and(b) the matrix in Exercise 2.15.Previous HomeCHAPTER 3INTERPOLATIONOften we need to approximate a function f(x) by another “more convenient” functionF(x). This arises in physical processes where f(x) is known only through its valuesat certain sample points x and F(x) is needed to approximate maximum or minimumvalues, to estimate integrals or derivatives, or merely to generate values for f(x) atpoints where experimental data are not available.
This also arises when f(x) is knownbut is difficult to evaluate, integrate, or differentiate. A familiar example is the functionf(x) = sinx that has to be approximated in a way that can be evaluated by calculatorsand computers. This is a fundamental principle of numerical analysis: if we cannotcarry out a basic computation with the function of interest, we approximate it by afunction for which we can do the computation.In this chapter a function f(x) is approximated by an interpolant F(x), a functionthat agrees with f(x) at certain points. A formal definition of the verb interpolate is asfollows.Definition.
interpolation:points {x l , . . . ,xN) ifA function F(x) is said to interpolate f(x) at theF(x j ) = f(x j ),j = 1,2 ,..., N.The process of constructing such a function F(x) is called interpolation.There are many types of approximating functions F(x) and which one to use depends to a large extent on the nature of the data and the intended use of the approximation. Perhaps the simplest approximating functions are polynomials. It can be shownthat any continuous function can be approximated arbitrarily well over a finite interval by a polynomial. More to the point, polynomials and their ratios (called rationalfunctions) are the only functions that can be evaluated directly on a computer. Forthis reason polynomials are used not only for interpolation but also as a foundation formost of the methods in the remaining chapters of the book.
Polynomial splines, that is,piecewise polynomial functions, are a very powerful tool for approximating functionsand are the main object of study in this chapter. In many applications the appearanceof the graph of F(x) is of great importance. For this reason it is very helpful to havea graphing package for visualization of the approximating functions derived in thischapter.82Next3.1 POLYNOMIAL INTERPOLATION83For a more thorough treatment of the theory of polynomial interpolation see [15,Chapters 5 and 6] and for more about approximation theory see [4]. The book [3] isan excellent introduction to polynomial splines and contains many FORTRAN codes.3.1 POLYNOMIAL INTERPOLATIONIn this section the approximation F(x) is a polynomial and it is traditional to use thenotation PN instead of F. The interpolation problem, formally stated, is as follows.Given the ordered pairs (x j ,f j ) for j = 1,2,.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
