Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 78
Текст из файла (страница 78)
The condition numberfor the harmonic points 1/(i + 1) grows faster than n!; by contrast, the condition numbers of the notoriously ill-conditioned Hilbert and Pascal matricesgrow only exponentially (see 26.1 and 26.4). For any choice of points therate of growth is at least exponential (V2), and this rate is achieved for pointsequally spaced on [0, 1]. For points equally spaced on [– 1, 1], the conditionnumber grows at a slower exponential rate than that for [0, 1], and the growthrate is slower still for the zeros of the nth degree Chebyshev polynomial (V6).For one set of points the Vandermonde matrix is perfectly conditioned: theroots of unity, for whichis unitary.21.2.
Primal and Dual SystemsThe standard Vandermonde matrix can be generalized in at least two ways:by allowing confluence of the points αi and by replacing the monomials by21.2 P RIMALANDD UAL S YSTEMS429other polynomials. An example of a confluent Vandermonde matrix is(21.4)The second, third, and fifth columns are obtained by “differentiating” theprevious column. The transpose of a confluent Vandermonde matrix arisesin Hermite interpolation; it is nonsingular if the points corresponding to the“nonconfluent columns” are distinct.A Vandermonde-like matrix is defined bywhere pi is apolynomial of degree i.
The case of practical interest is where the pi satisfy athree-term recurrence relation. In the rest of this chapter we will assume thatthe pi do satisfy a three-term recurrence relation. A particular application isthe solution of certain discrete Chebyshev approximation problems [11, 1984].Incorporating confluence, we obtain a confluent Vandermonde-like matrix,defined byP = P(α0, α1,. .
. ,αn ) = [q0 (α0), q1 (α1), . . . . qn (α n )]where the αi are ordered so that equal points are contiguous, that is,(21.5)and the vectors qj(x) are defined recursively byif j = 0 orotherwise.For all polynomials and points, P is nonsingular; this follows from the derivation of the algorithms below. One reason for the interest in Vandermonde-likematrices is that for certain polynomials they tend to be better conditionedthan Vandermonde matrices (see, for example, Problem 21.5). Gautschi [427,1983] derives bounds for the condition numbers of Vandermonde-like matrices.Fast algorithms for solving the confluent Vandermonde-like primal anddual systems Px = b and PTa = f can be derived under the assumption thatthe p j(x) satisfy the three-term recurrence relation(21.6a)(21.6b)wherefor all j.
Note that in this chapter γi denotes a constant in therecurrence relation and not iu/(1 – iu) as elsewhere in the book. The latternotation is not used in this chapter.V ANDERMONDE S YSTEMS430The algorithms exploit the connection with interpolation. Denote byr(i) > 0 the smallest integer for which αi = αi-1 = · · · = αr(i).
Considering first the dual system PTa = f, we note that(21.7)satisfiesi = 0:n.Thus ψ is a Hermite interpolating polynomial for the data {αi, fi }, and ourtask is to obtain its representation in terms of the basisAs a firststep we construct the divided difference form of ψ,(21.8)The (confluent) divided differences ci = f[α0, α1,. . . .
αi] may be generatedusing the recurrence relation(21.9)Now we need to generate the ai in (21.7) from the ci in (21.8). The idea isto expand (21.8) using nested multiplication and use the recurrence relations(21.6) to express the results as a linear combination of the pj. Defineq n (x) = cn ,(21.10)(21.11)from which q0(x) = ψ(x). Let(21.12)To obtain recurrences for the coefficientsof (21.11), givingwe expand the right-hand side21.2 P RIMALANDD UAL SYSTEMS431Using the relations, from (21.6),we obtain, for k = 0:n – 2,(21.13)in which the empty summation is defined to be zero. For the special casek = n – l we have(21.14)Recurrences for the coefficientsin terms ofj = k +1: n,follow immediately by comparing (21.12) with (21.13) and (21.14).In the following algorithm, stage I computes the confluent divided differences and stage II implements the recurrences derived above.Algorithm 21.2 (dual, PTa = f).
Given parametersa vecsatisfying (21 .5), this algorithm solves thetor f, and pointsdual confluent Vandermonde-like system PTa = f.% Stage I:Set c = fV ANDERMONDE S YSTEMS432for k = 0:n – 1clast = ckfor j = k + l:ni f αj αj-k-1thenc j = c j/(k + 1)elsetemp = cjc j = (cj – clast)/(αj – αj – k – 1 )clast = tempendendend% Stage II:Set a = ca n - l = a n - l + (β 0 - αn - 1 )a na n = a n /θ 0for k = n – 2: – 1:0for j = l:n – k –2enda n = a n /θ n - k - 1endAssuming that the valuesare given (note that γj appears only inthe termsthe computational cost of Algorithm 21.2 is at most 9n 2 / 2flops.
The vectors c and a have been used for clarity; in fact both can bereplaced by f, so that the right-hand side is transformed into the solutionwithout using any extra storage.Values offor some polynomials of interest are collected in Table 21.2.The key to deriving a corresponding algorithm for solving the primal system is to recognize that Algorithm 21.2 implicitly computes a factorization ofP - T into the product of 2n triangular matrices. In the rest of this chapterwe adopt the convention that the subscripts of all vectors and matrices runfrom 0 to n.
In stage I, letting c(k) denote the vector c at the start of the k t hiteration of the outer loop, we havec(0) = f,c(k+l) = Lkc(k),k = 0:n – 1.(21.15)The matrix Lk is lower triangular and agrees with the identity matrix in rows21.2 P RIMALAND433D UAL SYSTEMSTable 21.2. Parameters in the three-term recurrence (21.6).PolynomialMonomialsChebyshevLegendre*HermiteLaguerreθj12*2βj00002j + lγj01∗*θ0 = 1p j (1) = 12j0 to k. The remaining rows can redescribed, for k + l < j < n, byif αj = αj - k - 1 ,some s < j, otherwise,where ej is column j of the identity matrix.
Similarly, stage II can be expressedasa (n) = c (n) ,a ( k ) = U k a ( k + l ) , k = n – 1: – l:0.(21.16)The matrix Uk is upper triangular, it agrees with the identity matrix in rows0 to k – 1, and it has zeros everywhere above the first two superdiagonals.From (21.15) and (21.16) we see that the overall effect of Algorithm 21.2is to evaluate, step-by-step, the producta = U 0 . . .U n-l L n-l . . . L 0 f = P-T f.(21.17)Taking the transpose of this product we obtain a representation of P-1, fromwhich it is easy to write down an algorithm for computing x = P–1b.Algorithm 21.3 (primal, Px = b).
Given parametersa vecsatisfying (21.5), this algorithm solves thetor b, and pointsprimal confluent Vandermonde-like system Px = b.l% Stage I:Set d = bfor k = 0:n – 2for j = n – k: – 1:2endend434V ANDERMONDE S YSTEMS% Stage II:Set x = dfor k = n – 1: – 1:0xlast = 0for j = n: - l :k + li f αj = αj–k–l thenx j = x j/(k + 1)elset e m p = x j/(α j - αj - k - 1 )xj = temp – xlastxlast = tempendendxk = xk - xlastendAlgorithm 21.3 has, by construction, the same operation count as Algorithm 21.2.21.3. StabilityAlgorithms 21.2 and 21.3 have interesting stability properties. Depending onthe problem parameters, the algorithms can range from being very stable (ineither a backward or forward sense) to very unstable.When the pa are the monomials and the points αi are distinct, the algorithms reduce to those of Björck and Pereyra [121, 1970].
Björck and Pereyrafound that for the system Vx = b with αi = 1/(i + 3), bi = 2-i, n = 9, andon a computer with u 10 –16 ,Thus the computed solution has a tiny componentwise relative error, despitethe extreme ill condition of V. Björck and Pereyra comment “It seems as ifat least some problems connected with Vandermonde systems, which traditionally have been considered too ill-conditioned to be attacked, actually canbe solved with good precision.” This high accuracy can be explained with theaid of the error analysis below.The analysis can be kept quite short by exploiting the interpretation ofthe algorithms in terms of matrix–vector products.
Because of the inherentduality between Algorithms 21.2 and 21.3, any result for one has an analoguefor the other, so we will consider only Algorithm 21.2.The analysis that follows is based on the model (2.4), and so is valid onlyfor machines with a guard digit. With the no-guard-digit model (2.6) the21.3 S TABILITY435bounds become weaker and more complicated, because of the importance oftermsin the analysis.21.3.1. Forward ErrorTheorem 21.4. If no underflow or overflows are encountered then Algorithm 21.2 runs to completion and the computed solution satisfies(21.18)where c(n, u) := (1 + u)7n – 1 = 7nu + 0(u2).Proof. First, note that Algorithm 21.2 must succeed in the absence ofunderflow and overflow, because division by zero cannot occur.The analysis of the computation of the c(k) vectors is exactly the sameas that for the nonconfluent divided differences in 5.3 (see (5.9) and (5.10)).However, we obtain a slightly cleaner error bound by dropping the γ k notationand instead writing|∆Lk| < [(1 + u)3 - 1]|Lk|.(21.19)Turning to the equations (21.16), we can regard the multiplication a(k) =U ka (k+l) as comprising a sequence of three-term inner products.
Analysingthese in standard fashion we arrive at the equation|∆U k| < [(1 + u)4 - 1]|Uk|,(21.20)where we have taken into account the rounding errors in formingSince(21.19) and (21.20) imply that(21.21)Applying Lemma 3.7 to (21.21) and using (21.17), we obtain the desired boundfor the forward error.The product |U0| . . . |Un-1||Ln-1| . . . |L0| in (21.18) is an upper bound for|U0. . .Un-1Ln-1 . . . L0| = |P–T| and is equal to it when there is no subtractive cancellation in the latter product.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
