Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 96
Текст из файла (страница 96)
For any Aeigenvalues as A (to be precise, compan(A) is normal if and only if a0 = 1 andai = 0 for i > 0).Companion matrices tend to have interesting pseudospectra, as illustratedin Figure 26.2. For more information on the pseudospectra of companionmatrices see Toh and Trefethen [1007, 1994].26.7. Notes and ReferencesThe earliest reference to the Hilbert matrix appears to be [570, 1894], whereinHilbert obtains the formula (26.2) for det(H n ).The formula (26.1) for Hn-l is taken from Choi [207, 1983], who describesvarious interesting properties of the Hilbert matrix and its infinite analogue.An excellent reference for the derivation of formulae for the inverse of theCauchy and Hilbert matrices is Knuth [667, 1973, pp.
35–37]. Another reference for the Cauchy matrix is Tyrtyshnikov [1034, 1991]. The formulae (26.3)and (26.4) for the Cholesky factor and its inverse are from Choi [207, 1983]and Todd [1004, 1954].Forsythe and Moler [396, 1967] have a chapter devoted to the Hilbert matrix, in which they describe the underlying least squares problem and discussnumerical computation of the inverse. There have been many papers on thePROBLEMS527Hilbert matrix; two of particular interest are by Todd [1004, 1954], [1005,1961].Other references on the eigenvalues and condition numbers of randommatrices include Edelman [336, 1991], [339, 1992] and Kostlan [671, 1992].Anderson, Olkin, and Underhill [20, 1987] suggest another way to construct random orthogonal matrices from the Haar distribution, based on products of random Givens rotations. Marsaglia and Olkin [728, 1984] discuss thegeneration of random correlation matrices (symmetric positive semidefinitematrices with ones on the diagonal).The involuntary triangular matrix pascal (n, 1) arises in the step-size changing mechanism in an ordinary differential equation code based on backwarddifferentiation formulae; see Shampine and Reichelt [912, 1995].We mention some other collections of test matrices.
The Harwell-Boeingcollection of sparse matrices, largely drawn from practical problems, is presented by Duff, Grimes, and Lewis [327, 1989], [328, 1992]. Bai [36, 1994] isbuilding a collection of test matrices for large-scale nonsymmetric eigenvalueproblems. Zielke [1128, 1986] gives various parametrized rectangular matricesof fixed dimension with known generalized inverses. Demmel and McKenney [299, 1989] present a suite of Fortran 77 codes for generating randomsquare and rectangular matrices with prescribed singular values, eigenvalues,band structure, and other properties.
This suite was the inspiration for therandsvd routine in the Test Matrix Toolbox and is part of the testing codefor LAPACK (see below).26.7.1. LAPACKThe LAPACK distribution contains a suite of routines for generating test matrices, located in the directory LAPACK/TESTING/MATGEN (in Unix notation).These routines (whose names are of the form xLAxxx) are used for testingwhen LAPACK is installed and are not described in the LAPACK Users’Guide [17, 1995].Problems26.1. Investigate the spectral and pseudospectral properties of pentadiagnal Toeplitz matrices. See Figure 26.3 (pentoep is a function in the TestMatrix Toolbox that generates matrices of this type).
References are Beamand Warming [85, 1993] and Reichel and Trefethen [866, 1992].26.2. (R ESEARCH PROBLEM ) Two methods for generating a real orthogonal matrix from the Haar distribution are Stewart’s method (Theorem 26.1),based on Householder transformations, and a method of Anderson, Olkin,and Underhill [20, 1987], based on Givens rotations. Compare the efficiencyA G ALLERY528OFT EST M ATRICESpentoep(32,0,1/2,0,0,1)inv(pentoep(32,0,1,1,0,.25))pentoep(32,0,1/2,1,1,1)pentoep(32,0,1,0,0,1/4)Figure 26.3.
Pseudospectm of 32 x 32 pentadiagonal Toeplitz matrices.of these two methods when used to generate randsvd matrices. Investigatethe use of the Givens rotations approach to construct random banded matrices with given singular values and random symmetric banded matrices withgiven eigenvalues, and compare with the technique of generating a full matrixand then using band reduction (as implemented in the routine randsvd in theTest Matrix Toolbox).26.3.
(RESEARCH PROBLEM) Develop an efficient algorithm for computing aunit upper triangular n x n matrix with prescribed singular values σ1, . . . . σ n ,wherePrevious HomeAppendix ASolutions to Problemswe have1.1 SinceHence if< 0.01, say, then there is no difference between Erel andpractical purposes.for1.2. Nothing can be concluded about the last digit before the decimal point. Evaluating y to higher precision yieldst3540y262537412640768743.99999999999925007262537412640768743.9999999999992500725972This shows that the required digit is, in fact, 3. The interesting fact that y is so closeto an integer was pointed out by Lehmer [697, 1943], who explains its connectionwith number theory.
It is known that y is irrational [959, 1991].1.3.1.2. 2sin3. ( x – y)(x + y). Cancellation has not been avoided, but it is now harmless if xand y are known exactly (see also Problem 3.9).4. sin x/(1 + cos x).5. c = ((a – b)2 + ab(2sinθ/ 2 ) 2 ) 1 / 2 .1.4. a + ib = (x + iy)2 = x2 – y2 + 2ixy, so b = 2xy and a = x2 – y2, givingx2 – b2 /(4x 2) = a, or 4x4 – 4ax2 – b2 = 0.
HenceIf a > 0 we use this formula with the plus sign, since x2 > 0. If a < 0 this formulais potentially unstable, so we use the rewritten form529NextS OLUTIONS530TOP ROBLEMSIn either case we get two values for x and recover y from y = b /(2x ). Note thatthere are other issues to consider here, such as scaling to avoid overflow.1.5. We need a way to compute f(x) = log(1 + x) accurately for all x > 0. Astraightforward evaluation of log(1 + x) is not sufficient, since the addition 1 + xloses significant digits of x when x is small. The following method is effective (foranother approach, see Hull, Fairgrieve, and Tang [592, 1994]): calculate w = 1 + xand then computeThe explanation of why this method works is similar to that for the method in§1.14.1.
We have = (1 + x)(1 + δ), |δ| < u, and if= 1 then |x| < u + u2+ u3+ · · ·,2so from the Taylor series f(x) = x(1 – x/2 + x /3 – · · ·) we see that f(x) = x is acorrectly rounded result. IfthenDefining=: 1 + z,so if(thus x0) thenThus, with 1 + θ :=showing thatis an accurate approximation to f(x) = log(1 + x).1.7. From (1.4) we haveS OLUTIONSTO531P ROBLEMSThis yields the following inequality, which is attainable to first order:Dividing byThe result forand taking the limit asfollows fromgives the expression forwhere the inequality is attainable (Cauchy–Schwarz), together with the relationThatis easily verified.1.8.
The general solution to the recurrence iswhere a, b, and c are arbitrary constants. The particular starting values chosenyield a = 0, b = c = 1, so thatRounding errors in the evaluation (and in the representation of x1 on a binarymachine) cause a nonzero “a” term to be introduced, and the computed values areapproximately of the formfor a constant η of order the unit roundoff. Hence the computed iterates rapidlyconverge to 100.
Note that resorting to higher precision merely delays the inevitableconvergence to 100. Priest [844, 1992, pp. 54–56] gives some interesting observationson the stability of the evaluation of the recurrence.1.9. WritingC:= adj(A) =we have532S OLUTIONSTOP ROBLEMSso– (1/d) A ∆Cb, so |r| < γ3 |A||A-l||b|, which implies the normwiseAlso, r = b –residual bound. Note that the residual bound shows thatwillbe small if x is a large-normed solution1.10.
LetThen, sinceFor any standard summation method we have (see §4.2)satisfiesHence, definingthat is,where2.1. There are 1 + 2(e max – emin + 1)(β – 1)β t-1 normalized numbers (the “1” is forzero), and 2(β t–1 – 1) subnormal numbers. For IEEE arithmetic we therefore haveS OLUTIONSTO533P ROBLEMSsingle precisiondouble precisionNormalized Subnormal1.7 x 1074.3 x 109191.8 x 109 x 10152.2.
Without loss of generality suppose x > 0. We can write x = m x β e-t, whereβ t-l < m < β t. The next larger floating point number is x + ∆x, where ∆x = β e - t ,andThe same upper bound clearly holds for the “next smaller” case, and the lowerbound in this case is also easy to see.2.3. The answer is the same for all adjacent nonzero pairs of single precision numbers. Suppose the numbers are 1 and 1 +(single) = 1 + 2-23.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
