Heath - Scientific Computing (523150), страница 24

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

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

Resolving the member forces into horizontal√ and vertical components and definingα = 2/2, we obtain the following system ofequations for the member forces fi :f2 = f6Joint 2 : f3 = 10αf1 = f4 + αf5Joint 3 :αf1 + f3 + αf5 = 0f4 = f8Joint 4 : f7 = 0αf5 + f6 = αf9 + f10Joint 5 : αf5 + f7 + αf9 = 15f10 = f13Joint 6 : f11 = 20f8 + αf9 = αf12Joint 7 : αf9 + f11 + αf12 = 0Joint 8 : f13 + αf12 = 0Use a library routine to solve this system oflinear equations for the vector f of memberforces. Note that the matrix of this systemis quite sparse, so you may wish to experiment with a band solver or more general sparsesolver, although this particular problem instance is too small for these to offer significantadvantage over a general solver.2.4 Write a routine for estimating the condition number of a matrix A.

You may use eitherthe 1-norm or the ∞-norm (or try both andcompare the results). You will need to compute kAk, which is easy, and estimate kA−1 k,which is more challenging. As discussed in Section 2.3.3, one way to estimate kA−1 k is topick a vector y such that the ratio kzk/kyk islarge, where z is the solution to Az = y. Trytwo different approaches to picking y:(a) Choose y as the solution to the systemAT y = c, where c is a vector each of whosecomponents is ±1, with the sign for each component chosen by the following heuristic. Using the factorization A = LU , the systemAT y = c is solved in two stages, successivelysolving the triangular systems U T v = c andLT y = v.

At each step of the first triangularsolution, choose the corresponding component79of c to be 1 or −1, depending on which willmake the resulting component of v larger inmagnitude. (You will need to write a customtriangular solution routine to implement this.)Then solve the second triangular system in theusual way for y. The idea here is that anyill-conditioning in A will be reflected in U , resulting in a relatively large v. The relativelywell-conditioned unit triangular matrix L willthen preserve this relationship, resulting in arelatively large y.(b) Choose some small number, say, five, different vectors y randomly and use the one producing the largest ratio kzk/kyk. (For this youcan use an ordinary triangular solution routine.)You may use a library routine to obtain thenecessary LU factorization of A.

Test yourprogram on the following matrices:0.641 0.242A=,0.321 0.12110 −7 0B =  −32 6.5 −1 5How do the results using these two methodscompare? To check the quality of your estimates, compute A−1 explicitly to determineits true norm (this computation can also makeuse of the LU factorization already computed).If you have access to linear equations softwarethat already includes a condition estimator,how do your results compare with its?2.5 (a) Use a single-precision routine forGaussian elimination to solve the systemAx = b, where21.0 67.0 88.0 73.07.0 20.0  76.0 63.0A=,0.0 85.0 56.0 54.019.3 43.0 30.2 29.4141.0 109.0 b=.218.093.7(b) Compute the residual r = b − Ax usingdouble-precision arithmetic, if available (butstoring the final result in a single-precision vector r).

Note that the solution routine may80CHAPTER 2. SYSTEMS OF LINEAR EQUATIONSdestroy the array containing A, so you mayneed to save a separate copy for computingthe residual. (If only one precision is availablein the computing environment you use, thendo all of this problem in that precision.)(c) Solve the linear system Az = r to obtainthe “improved” solution x + z. Note that Aneed not be refactored.(d ) Repeat steps b and c until no further improvement is observed.2.6 An n × n Hilbert matrix H has entrieshij = 1/(i + j − 1), so it has the form1 12 13 · · ·111 12 13 14 · · · .···3 4 5.. ..

.. . ... . .For n = 2, 3, . . ., generate the Hilbert matrix of order n, and also generate the n-vectorb = Hx, where x is the n-vector with all of itscomponents equal to 1. Use a library routinefor Gaussian elimination (or Cholesky factorization, since the Hilbert matrix is symmetricand positive definite) to solve the resulting linear system Hx = b, obtaining an approximatesolution x̂.

Compute the ∞-norm of the residual r = b − H x̂ and of the error ∆x = x̂ − x,where x is the vector of all ones. How largecan you take n before the error is 100 percent(i.e., there are no significant digits in the solution)? Also use a condition estimator to obtaincond(H) for each value of n. Try to characterize the condition number as a function ofn. As n varies, how does the number of correct digits in the components of the computedsolution relate to the condition number of thematrix?2.7 (a) What happens when Gaussian elimination with partial pivoting is used on a matrixof the following form?1000 1 −1100 1  −1 −1101  −1 −1 −11 1 −1 −1 −1 −1 1Do the entries of the transformed matrix grow?What happens if complete pivoting is used instead? (Note that part a does not require acomputer.)(b) Use a library routine for Gaussian elimination with partial pivoting to solve varioussizes of linear systems of this form, using righthand-side vectors chosen so that the solutionis known.

How do the error, residual, and condition number behave as the systems becomelarger? This artificially contrived system illustrates the worst-case growth factor cited inSection 2.4.1 and is not indicative of the usualbehavior of Gaussian elimination with partialpivoting.2.8 Multiplying both sides of a linear system Ax = b by a nonsingular diagonal matrixD to obtain a new system DAx = Db simply rescales the rows of the system and in theory does not change the solution. Such scalingdoes affect the condition number of the matrixand the choice of pivots in Gaussian elimination, however, so it may affect the accuracyof the solution in finite-precision arithmetic.Note that scaling can introduce some rounding error in the matrix unless the entries ofD are powers of the base of the floating-pointarithmetic system being used (why?).Using a linear system with randomly chosenmatrix A, and right-hand-side vector b chosen so that the solution is known, experimentwith various scaling matrices D to see whateffect they have on the condition number ofthe matrix DA and the solution given by alibrary routine for solving the linear systemDAx = Db.

Be sure to try some fairly skewedscalings, where the magnitudes of the diagonal entries of D vary widely (the purpose isto simulate a system with badly chosen units).Compare both the relative residuals and theerror given by the various scalings. Can youfind a scaling that gives very poor accuracy?Is the residual still small in this case?2.9 (a) Use Gaussian elimination withoutpivoting to solve the linear system 11 1x1x21+=2for = 10−2k , k = 1, .

. . , 10. The exact soluTtion is x = [ 1 1 ] , independent of the valueof . How does the accuracy of the computedsolution behave as the value of decreases?COMPUTER PROBLEMS(b) Repeat part a, still using Gaussian elimination without pivoting, but this time use oneiteration of iterative refinement to improve thesolution, computing the residual in the sameprecision as the rest of the computations. Nowhow does the accuracy of the computed solution behave as the value of decreases?2.10 Consider the linear system 11+x11 + (1 + )=,1−1x21where is a small parameter to be specified.The exact solution is obviously 1x=for any value of .Use a library routine based on Gaussian elimination to solve this system.

Experiment withvarious values for , especially values near√mach for your computer. For each value of you try, compute an estimate of the conditionnumber of the matrix and the relative errorin each component of the solution. How accurately is each component determined? Howdoes the accuracy attained for each component compare with expectations based on thecondition number of the matrix and the errorbounds given in Section 2.4.2? What conclusions can you draw from this experiment?2.11 (a) Write programs implementing Gaussian elimination with no pivoting, partial pivoting, and complete pivoting.(b) Generate several linear systems with random matrices (i.e., use a random number generator to obtain the matrix entries) and righthand sides chosen so that the solutions areknown, and compare the accuracy, residuals,and performance of the three implementations.(c) Can you devise a (nonrandom) matrix forwhich complete pivoting is significantly moreaccurate than partial pivoting?2.12 Write a routine for solving tridiagonalsystems of linear equations using the algorithmgiven in Section 2.5.3 and test it on some sample systems.

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

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