accuracy of Gaussian quadrature

Accuracy of Gaussian Quadrature PhL 1.15.04 Contents: Theorem 1: The Gaussian Quadrature n+1 point formula is valid for an arbitrary polynomial of degree n where any distinct xi may be used. Theorem 2. If we assume that the n+1 point GQ formula is valid for polys of some degree M, then for such polys we have two equivalent ways to state the Gauss formula, and either implies the other: dx w(x)f(x) = i=0n Ai f(xi) (1) i=0n Ai pk(xi) = k,0 (h0/k0) k = 0,1,2....M (2) "sum rule" Theorem 3. If xi are the zeros of pn+1(xi), then the above "sum rule" is valid for a larger range of k: i=0n Ai pk(xi) = k,0 (h0/k0) k = 0,1,2....2n+1 (2) Theorem 4: If the xi are chosen as the zeros of the orthog poly pn+1(x), then the n+1 point Gaussian quadrature formula is true for any poly of degree 2n+1. [ Otherwise, it is only true up to degree n. ] this theorem follows at once from Theorem 3 and Theorem 2. Appendix A Theorem: There is a 1-to-1 correspondence between degree n polynomials written in the form i=0n ai xi and those written in the form i=0n Li(x)f(xi) for any set of n+1 distinct points { xi }. ____________________________________________________________________________________ Theorem 1: The Gaussian Quadrature n+1 point formula is valid for an arbitrary polynomial of degree n. Proof: Start with the formula where we use an arbitrary set of n+1 distinct points xi in our interval: dx w(x)f(x) = i=0n Ai f(xi) " Gaussian formula" with n+1 points where Ai = dx w(x)Li(x) We know from Appendix A that we can represent an arbitrary degree n poly as (it has n+1 params) f(x) = i=0n Li(x)f(xi). Then here is a proof of the Gauss formula for f(x) a arbitrary degree n poly: dx w(x) f(x) = dx w(x) {i=0n Li(x) f(xi)} = i=0n [ dx w(x) Li(x) ]f(xi) = i=0n Ai f(xi) Any distinct set {xi} of n+1 distinct points may be used. I don't think the xi have to be within the integration interval (a,b), but I might be wrong about that. Later when we decide to use the zeros of pn+1(x), they will be in the interval because that is where all n+1 zeros of pn+1 are located. ____________________________________________________________________________________ Theorem 2. If we assume that the Gauss formula is valid for polys of some degree M, then for such polys we have two equivalent ways to state the Gauss formula, and either implies the other: dx w(x)f(x) = i=0n Ai f(xi) (1) i=0n Ai pk(xi) = k,0 (h0/k0) k = 0,1,2....M (2) So far, we know only that (1) is true for M n, as shown in Theorem 1, but this fact does not interfere with our statement and proof of Theorem 2. In the proof that follows, pn(x) is the unique set of orthogonal polynomials (apart from normalization) that is defined by the interval (a,b) and weight function w(x). The integral in (1) is over this interval. We also use the scalar product notation (f,g). Proof: => First, we will show that (1) => (2): Using the Erdelyi notation and p0(x) = k0 and hn= (pn, pn) we know that dx w(x) p0(x) = (1,p0) = (1/p0)(p0, p0) = h0/k0 Since the pk are orthogonal, we know this is true for all k, dx w(x) pk(x) = (1,pk) = k,0 (h0/k0) = k = 0,1,2..... If we now assume that (1) is true, we can apply (1) in particular to f(x) = pk(x) for k = 0..M to get i=0n Ai pk(xi) = dx w(x)pk(x) = k,0 (h0/k0) k = 0,1,2....M and we have thus shown that (1) => (2) QED. => Second, we will show that (2) => (1): We expand an arbitrary degree M poly f(x) in the pk(x) as follows, f(x) = k=0M bk pk(x) // this is a poly of degree M We then integrate this to get an expression for the LHS of (1), dx w(x)f(x) = dx w(x){k=0M bk pk(x)} = k=0M bk (1,pk) = k=0M bk k,0 (h0/k0) = b0(h0/k0) and this is true for any positive integer M. So we have shown that LHS of (1) = b0(h0/k0). Now let's examine the RHS of (1): RHS = i=0n Ai f(xi) Insert the same expansion above for f(x) of degree M, RHS = i=0n Ai { k=0M bk pk(xi)} = k=0M bk [ i=0n Ai pk(xi) ] If we now assume (2) is valid for polys of degree M, we may make this replacement [i=0n Ai pk(xi)] = k,0 (h0/k0) k = 0,1,2....M and we find that RHS = k=0M bk [ i=0n Ai pk(xi) ] = k=0M bk k,0 (h0/k0) = b0(h0/k0) We have thus shown that RHS (1) = LHS (1) = b0(h0/k0), so (1) is true. QED. ____________________________________________________________________________________ Theorem 3. Consider the above version (2) of the n+1 point Gauss formula: i=0n Ai pk(xi) = k,0 (h0/k0) k = 0,1,2....M (2) We know from Theorem 1 that (2) is true for M n. The current theorem claims that, if we select the points xi to be the zeros of pn+1(xi), then (2) is true for the extended range M 2n+1. Proof: This proof uses the following ingredients: (a) the fact that pn+1(xi) = 0 since we are selecting the n+1 points xi to be the n+1 zeros of pn+1(x) (b) the recursion formula which schematically reads pk+1 = xpk + pk + pk-1 (c) the "angular momentum" theorem which says schematically that xpk = pk+1 + pk + pk-1 By "schematically" we mean we are not showing various constants like An, Bn, Cn in the recursion relation, or coefficients in the product expansion. We set all constants to 1 because we are only going to be concerned with whether terms are zero or non-zero; we won't be caring about exact values of non-zero terms. We start with a derivation of item (c) above. The "angular momentum" product rule (see my Erdelyi Addendum notes, this result does not appear in Erdelyi's book) states: pj1(x)pj2(x) = !Syntax Error, Iaj pj(x) and the proof of this fact you will see there is quite trivial and is based on nothing but the basic properties of the orthogonal polys, in particular, their orthogonality. We use integers j1, j2 and j for label names as a reminder that this is very similar to the "angular momentum addition theorem" in SO(3) group theory. Knowledge of group theory of course is irrelevant to what we are doing here, and the reader is welcome to ignore the phrase "angular momentum". The above expansion is just a simple fact. We now apply this fact to claim that x pk = (p0 + p1)pk = (p0pk) + (p1pk) = (pk) + (pk+1 + pk + pk-1) = pk+1 + pk + pk-1 where again we are setting all constants to 1. Certainly we can write x as a linear combination of p0 and p1with some specific coefficients. Clearly in this notation we have p0pk = pk by our theorem. And the second term is what gives the interesting part, p1pk = pk+1 + pk + pk-1. So the effect of applying x to pk is toe create a "spread" of unit in each direction. Now let's indicate by { pa pb } a descending range of pn . For example, pk+1+pk+pk-1 = { pk-1pk-1}. The range need not be "full" (as it is in our example), we just mean that a is the largest index and b is the smallest index in a sum of terms. Then rule (c) xpk = pk+1 + pk + pk-1 tells us that x{ pa pb } = { pa+1 pb-1 } (d) so application of x "stretches the range" by one unit in at each end. Now we begin our process of expanding the range for which (2) is true. As noted, we already know it is true for k = 0,1.2..n. i=0n Ai pk(xi) = k,0 (h0/k0) k = 0,1,2....n (2) It is obviously true for k = n+1 because by (a) we know pn+1(xi) = 0. For pn+2 we use the recursion (b) which says pn+2 = xipn+1 + pn+1 + pn where we have x = xi everywhere. This simplifies to pn+2 = pn since pn+1 = 0 and we already know that (2) is true for k=n, so now we know that it is true for k = n+2 as well. There is one exception, however. If we have n=0, meaning a 1-point Gaussian formula, then pn+2 = pn => p2 = p0, and we know that (2) 0 for p0, and therefore (2) is non-zero for p2 in this case. This agrees with the claim that the result is valid only for k 2n+1 =1. So we now assume that n > 0 and continue: [ By the way, the above case is the ONLY place we need to use the fact that pn+1(xi) = 0. It allows us to "get started". Accordingly, we will not delete pn+1 terms in the following, just to make this fact clear. ] For pn+3 the recursion now says pn+3 = xipn+2 + pn+2 + pn+1. But from above we know that pn+2 = pn, so we get pn+3 = xipn + pn+2 + pn+1 = { pn+1 pn-1 } + pn+2 + pn+1 = { pn+2 pn-1} where we have used rule (d) to expand xpn . The "danger" in this expression is the term pn-1. If we have n=1, then this is p0 and (2) is non zero for pn+3 = p4. So if n = 1, we would conclude that our theorem was not true for p4, and this agrees with our claimed validity limit of 2n+1 = 3. So we now assume that n > 1 and continue: For pn+4 the recursion now says pn+4 = xipn+3 + pn+3 + pn+2 . We insert pn+3 = { pn+2 pn-1} to get that pn+4 = xi{ pn+2 pn-1} + pn+3 + pn+2 = { pn+3 pn-2}. The higher terms all give 0 in (2), but now we are concerned about the pn-2 term. But let's now jump to the general case. Here is the general pattern of the recursion formulas. The first three in this list we have explicitly derived. We see that the "range" increases by 1 in each direction for each step, pn+2 = { pn+1 pn-0 } 1 = 2-1 0 = 2-2 pn+3 = { pn+2 pn-1 } 2 = 3-1 1 = 3-2 pn+4 = { pn+3 pn-2 } 3 = 4-1 2 = 4-2 ...... pn+m = { pn+(m-1) pn-(m-2)} (m-1) = m-1 (m-2) = m-2 and the last line will be "invalid" when m gets up to a value such that pn-(m-2) = p0 . Again, this is so because (2) 0 when k =0. By "invalid" we mean that (2) will be not true when we insert pn+m for pk. The last line is clear from the pattern, but we nevertheless give an induction proof below. Assuming for now it is correct, we find that our largest valid value of m occurs when pn-(m-2) = p1 which says m = n+1 so that n+m = 2n+1. Thus, the largest polynomial for which (2) = 0 is pn+m = p2n+1 . We have then shown that (2) is true up to this point, that is, i=0n Ai pk(xi) = k,0 (h0/k0) k = 0,1,2....(2n+1) (2) and our theorem is proved QED. Here is our induction proof of the pn+m expression above. Assume it is true for m: pn+m = { pn+(m-1) pn-(m-2)} (3) Then for m+1 we have pn+(m+1) = xipn+m + pn+m + pn+m-1 = xi { pn+(m-1) pn-(m-2)} + pn+m + pn+m-1 = { pn+m pn-(m+1-2)} + pn+m + pn+m-1 // where we used rule (d) = { pn+m pn-(m+1-2)} = { pn+(m+1-1) pn-(m+1-2)} so we have shown that if (3) is true for m, it is also true for m+1, and of course we have shown already that (3) is true explicitly for m = 2,3 and 4. ____________________________________________________________________________________ Theorem 4: If the xi are chosen as the zeros of the orthog poly pn+1(x), then the Gaussian quadrature formula is true for any poly of degree 2n+1. [ Otherwise, it is only true up to degree n. ] Proof: We have shown in Theorem 1 that the Gaussian formula is generally true for degree n. We have shown in Theorem 3 that, when the xi are chosen such that pn+1(xi) = 0, the alternate form (2) of the Gauss formula is true for k = 0,1,2 up to (2n+1), that is to say i=0n Ai pk(xi) = k,0 (h0/k0) k = 0,1,2....(2n+1) (2) Theorem 2 then tells us that, since (2) is true for M = 2n+1, then (1) is also true for M = 2n+1, and (1) is the official Gaussian quadrature formula. QED. ____________________________________________________________________________________ Comment #1: Sometimes one uses the n-point Gaussian formula (instead of our n+1 point one). In that case, the general Gauss formula is true for polys up to degree n-1 (instead of n), and the extended formula is then true up to 2(n-1)+1 = 2n-1 (instead of 2n+1 ). Also, the xi are then the n zeros of pn (not pn+1). Comment #2: The formula i=0n Ai pk(xi) = k,0 (h0/k0) for k = 0,1,2....M (with Ai = dx w(x)Li(x) ) is really a statement about the orthogonal polynomial set { pk}. It has nothing to do with any function we are integrating. Notice that the xi appear inside Li(x) (see appendix below). If we change from one general set of xi to another, the Ai all change, but then of course so do the pk(xi), and these changes occur in a way that the sum i=0n Ai pk(xi) is unchanged and still equals k,0 (h0/k0). For general xi this peculiar property of the {pk} is true for k = 0 n, but for the specific xi which are the zeros of pn+1 it is true for the much larger range k = 0 2n+1. Therefore, the accuracy of the Gaussian formula applied as an approximation to an arbitrary function f(x) is obviously going to be much improved when we use xi which are these zeros. You can "fit" an arbitrary f(x) much better with a 2n+1 poly compared to an n poly. Scheid gives detailed error formulas which show that n+1 point error is proportional to f(2n+2)() where is some point in (a,b). Clearly for a polynomial of degree 2n+1, the error formula then gives 0 error. Deriving that error formula is another way to get the result we have shown above, but we wanted here to show directly how it is that those xi work with the recursion formula to expand the range of the formula. ____________________________________________________________________________________ Appendix A: Theorem: There is a 1-to-1 correspondence between degree n polynomials written in the form i=0n ai xi and those written in the form i=0n Li(x)f(xi) for any set of n+1 distinct points { xi }. Proof: First, we show that g(x) = i=0n Li(x)f(xi) is in fact a polynomial of degree n. The reason is that each factor Li(x) is a poly of degree n, and we are then adding up polys of degree n to get an overall poly of degree n. To see that Li(x) is degree n, just write it out: (x) = (x - x0) ... (x - xn) = degree n+1 Li(x) = (x)/[(x-xi)'(xi)] = degree n Ignoring constants, Li(x) is just the product of n factors (x - xk) where the ith factor (x - xi) is missing. Second, it is pretty obvious that i=0n ai xi is also a poly of degree n. Now, assume f(x) = i=0n ai xi . Suppose we specify n+1 distinct xk points, so that yk f(xk) = i=0n ai (xk)i or i=0n (xk)i ai = yk . We then have this matrix equation where X is n+1 square Xa = y where Xki = (xk)i As I have shown in the Erdelyi Addendum notes (and see page 158 top), detX 0 as long as the n+1 points are all different, a fact that is not obvious. Therefore, we can invert to get a = X-1 y So there is a 1-1 correspondence between the two ways you can parameterize a n degree poly: { ai } { yk }. Each set of course contains n+1 numbers. If the points are not distinct, then detX = 0 and then X has a nullspace and we might get Xa = 0 for some choices of a, so then y cannot describe an arbitrary polynomial. Corollary: An "arbitrary" polynomial f(x) of degree n can be written as i=0n Li(xk)f(xi). The corollary just says this: pick an arbitrary deg n polynomial by picking coefficients { ai }. We can then compute y = Xa, and then we know how to write that arbitrary deg n poly in the form i=0n Li(xk)yi .