erdelyi addendum

Erdelyi Addendum PhL 1.20.04 These notes contain derivations and facts on orthogonal polynomials which Erdelyi did not provide. These things were used in Jim Ball's Gaussian Quadrature papers, but are really ortho poly facts. Contents: Theorem A1: ( The X Determinant) Theorem A2: (Angular Momentum Product Rule for Ortho Polys) 1. About the scaling freedom of the ortho polys. 2. How to get the ortho poly "A,B,C" recursion relation into "symmetric form" (via orthonormals) 3. How to get the "A,B,C" recursion relation into "Gautschi ',' form" (via monics) 4. How to compute the ortho polys from the symmetric recursion relation coefficients. 5. How to compute the symmetric recursion relation coefficients and the ortho polys at the same time, given knowledge of the w(x) and (a,b). [ This is the 1884 Stieltjes procedure ] 6a. Compute the ortho poly recursion coefficients A,B,C as integrals of the polys. 6b. Repeat previous section for the "symmetric" form. 6c. Repeat previous section for the "Gautschi ',' form". _____________________________________________________________________________________ On page 158 we find the following claim as part of (3): Theorem A1: ( The X Determinant) Consider this n+1 x n+1 matrix X, x00, x01, .... x0n x10, x11, .... x1n ..... = X xn0, xn1, .... xnn We claim that det(X) = r>s (xr - xs ). Example: = (x2 - x1)(x2- x0)(x1- x0) Proof: Let's do the proof for n+1 = 5 (n = 4), and it will be obvious what the general case looks like. Here is our matrix in question: x00, x01, x02, x03, x04 x10, x11, x12, x13, x14 x20, x21, x22, x23, x24 x30, x31, x32, x33, x34 x40, x41, x42, x43, x44 Our variables are x0, x1.....x4. Let's think first of x4 as our "variable". We know that detX is a polynomial in x4 of degree 4 ( just do the Laplace expansion along the bottom row). We know that such a poly can be written in the (x) form as a product of four factors (x4 - ai), so we have: detX = k(x0, x1, ..x3) (x4- a1)(x4- a2)(x4- a3)(x4- a4) where the roots are some a1 through a4, and where the leading constant is a function of the "other" variables. Now we know that detX = 0 if x4 = any of the other xi because then two rows are the same in the matrix. Thus, the roots ai must be "the other xi". So we now have detX = k(x0, x1, ..x3) (x4- x3)(x4- x2)(x4- x1)(x4- x0) Now consider this object k(x0, x1, ..x3) and let us ponder x3 as our variable now. We know that overall detX is a poly of degree 4 in x3. We have exactly one factor showing already which is (x4- x3), so it must be that k(x0, x1, ..x3) is a poly of degree 3 in x3. In other words, we already have one of the four factors, so it must be that k(x0, x1, ..x3) = k'(x0, x1, ..x2)(x3- x2)(x3- x1)(x3- x0) where we have now anticipated that the three roots for the x3 factors are the "other xi", just as above. Next, we recur again to say that k'(x0, x1, ..x2) = k"(x0, x1) (x2- x1)(x2- x0) by the same logic. In other words, for variable x2 we already extracted 2 of the 4 roots in our previous 2 steps [ namely, (x4- x2) in the first step, and (x3- x2) in the second step ] so we then write down the remaining two roots as above, (x2- x1)(x2- x0). Now we are close to the end with the next step: k"(x0, x1) = k'''(x0) (x1- x0) And we now assemble our total answer: detX = k'''(x0) [ (x4- x3)(x4- x2)(x4- x1)(x4- x0) ] [ (x3- x2)(x3- x1)(x3- x0) ] [ (x2- x1)(x2- x0) ] [(x1- x0) ] Now if we write det(X) in the form we get this det(X) = abcde (xa)0(xb)1(xc)2(xd)3(xe)4 = (x0)0(x1)1(x2)2(x3)3(x4)4 + other similar terms. The one term we have exposed is the product of the diagonal elements. We can identify this term in the above product of factors! The first four factors [ (x4- x3)(x4- x2)(x4- x1)(x4- x0) ] give us our (x4)4 and the next three [ (x3- x2)(x3- x1)(x3- x0) ] give us (x3)3, and so on. We know the coefficient of this one term must be 1, so we conclude that k'''(x0) = 1. So our final answer is: detX = [ (x4- x3)(x4- x2)(x4- x1)(x4- x0) ] [ (x3- x2)(x3- x1)(x3- x0) ] [ (x2- x1)(x2- x0) ] [(x1- x0) ] = r>s (xr - xs ) where r = 4,3,2,1 QED Comments: I don't think I have seen this theorem anywhere before. It is certainly not in any book I have (that I know of) except Erdelyi vol 2 page 158 top. It probably has a name associated with it, but I don't know how to find it because I don't know how to look this up on the web. Corollary: As long as all the xi are different, we have detX = 0 so X is invertible. Comment: There is a similar looking thing on page 53 of Scheid. _____________________________________________________________________________________ Theorem A2: (Angular Momentum Product Rule for Ortho Polys) If we expand pj1(x)pj2(x)as a lincom of ortho polys, we get a restricted range exactly as in the angular momentum rule of SO(3). pj1(x)pj2(x) = !Syntax Error, Iaj pj(x) Proof: We know that we can expand like so: pj1(x)pj2(x) = !Syntax Error, Iaj pj(x) just because we know that the product is a poly of degree j1 + j2. Now close with pk on both sides to get (pk, pj1pj2) = !Syntax Error, Iaj (pk,pj) = akhk // due to orthogonality Now assume that j1 j2 and write (pk, pj1pj2) = (pkpj2, pj1) We know that this will be zero if k+j2 < j1 because if we expand pkpj2 in pn, the highest will be n = k + j2. And if k+j2 < j1, then all of these are orthogonal to pj1. Therefore, all ak = 0 for k < j1 - j2. Thus, the actual sum of non-vanishing terms begins at j = j1 - j2. We can then use this symbolic notation for a product. n m = (n+m) + (n+m-1) + .....+ ( | n - m | ) where n is a shorthand for pn(x) and where the addition on the right ignores the coefficients ai . _____________________________________________________________________________________ 1. About the scaling freedom of the ortho polys. Consider the Erdelyi recursion formula pn+1 = (Anx + Bn)pn - Cnpn-1 where we show three parameters (An,Bn,Cn). These are based on (kn, kn', hn ) as shown page 159. Now, suppose we have this set { kn, kn', hn, pn(x) } and we define Pn(x) = sn pn(x) where sn is an arbitrary scale factor distinct for each n. Then we have a new set { Kn, Kn', Hn, Pn(x) }. But we know that Kn = snkn, Kn' = snkn' Hn = sn2 hn Pn(x) = sn pn(x) just from the definitions of these quantities. Now suppose we decide to scale things such that Hn = 1, so we are orthonormalized. We then have Hn = 1 = sn2 hn so that sn = 1/. We then get Kn = kn /, Kn' = kn' / Hn = 1 Pn(x) = pn(x) / At this point, we have to take what we get for Kn, Kn' and Pn. We cannot "force" things to be monic, for example. OK, suppose instead we scale to get monics. then Kn = snkn = 1, so sn = 1/kn and then we have Kn =1, Kn' = kn'/kn Hn = hn /kn2 Pn(x) = pn(x)/kn In this situation, we get what we get for Hn and we cannot control it. We cannot force these to be both monic and orthonormal, obviously. _________________________________________________________________________ 2. Show how to get the ortho poly "A,B,C" recursion relation into "symmetric form" The "genuine" recursion relation is this and it goes with (pn, pn) = hn , ( Erdelyi page 158) pn+1 = (Anx + Bn)pn - Cnpn-1 As we will show here, if you change the scale by doing n n(x) = pn(x) such that the n(x) are orthonormal, you arrive at what I have been calling "the symmetric form" of the recursion relation, x n =nn+1 + n n + n-1 n-1 Rewrite this as x pn = (1/An)pn+1 - (Bn/An) pn + (Cn/An) pn-1 Define n n(x) = pn(x). Then we have n x n = n+1 (1/An)n+1 - n (Bn/An) n + n-1 (Cn/An) n-1 Divide by n to get x n = (n+1/n) (1/An)n+1 - (Bn/An) n + (n-1 /n) (Cn/An) n-1 Now define n = (n+1/n) (1/An). Then we want to have n-1 = (n-1 /n) (Cn/An) which says n = (n /n+1) (Cn+1/An+1) So we then require that (n+1/n) (1/An) = (n /n+1) (Cn+1/An+1) which says that (n+1/n) 2 = Cn+1 (An/An+1) = (hn+1/hn) n+1 = n We can now iterate to get 1 = 0 = 0 2 = 0 = 0 3 = 0 = 0 ..... n = 0 (n, m) = (1/n)2(pn, pm) = (1/n)2 hnn,m = h0 (1/0)2 n,m Now, if we are handed the An, Bn, Cn we can think of that as determining kn, kn' and hn, so we cannot willy nilly go adjusting h0 . However, there is no reason not to choose 0 such that 0 = and then we end up with a nice orthonormal set of ortho polys. Meanwhile, we can look back at n = (n+1/n) (1/An) = (1/An) = (kn/kn+1) = (1/An) = So we can now summarize all the results: Start: pn+1 = (Anx + Bn)pn - Cnpn-1 // Erdelyi notation Define: n(x) = pn(x)/n Have: x n = (n+1/n) (1/An)n+1 - (Bn/An) n + (n-1 /n) (Cn/An) n-1 Require: x n = n n+1 - (Bn/An) n + n-1 n-1 Solve: n+1 = n = n Iterate: n 2 = (A0/An)C1C2....Cn => n = 0 Find: (n, m) = h0 (1/0)2 n,m Find: n = = (kn/kn+1) Define: n = - (Bn/An) = ( rn - rn+1) = ( kn'/kn - kn+1'/kn+1) An = kn+1/kn End up: x n =nn+1 + n n + n-1 n-1 Bn = An(rn+1 - rn) Select: 02 = h0 => (n, m) = n,m orthonormal Cn = (An/An-1)(hn/hn-1) The conclusion is that there is a simple way to define the recursion coefficients to get the symmetric Jacobi matrix and orthonormal polys. We have shown in the box above how you do this if you are given either (An,Bn,Cn) or (hn, kn, kn') as your independent variables. Wilf states this in another way. Suppose we selected pn(x) with hn = 1 from the get-go, which means that our pn are then orthonormal. Then in the above we would get n = 1 and n = pn and we would at once have the symmetric form with n = (kn/kn+1) = 1/An and n = ( rn - rn+1). So he regards the symmetric form just as a characteristic of having orthonormal pn . _________________________________________________________________________ 3. Show how to get the A,B,C recursion relation into Gautschi ',' form (via monics) Our starting point as before is this x pn = (1/An)pn+1 - (Bn/An) pn + (Cn/An) pn-1 where An = kn+1/kn and Bn = An (rn+1 - rn) and Cn = (An/An-1)(hn/hn-1) so we regard (An,Bn,Cn) as functions of (hn, kn, kn'). Since p-1 = 0, we don't care about the number C0. We can then just set C0 = any value we like, perhaps 0 or perhaps 1. So we only use the formula Cn = (An/An-1)(hn/hn-1) for n = 1,2,3... whereas we use the other two formulas for An and Bn we use for n = 0,1,2... Now, with monics we have kn = 1 so that An = 1 and then Cn = (hn/hn-1) and Bn = (k'n+1 - k'n). Note that by requiring all monics, we have to accept hn and kn' being whatever they are, we cannot "set them". For monics, then, we get a simplified recursion formula x pn = pn+1 - Bn pn + Cn pn-1 pn = monic x pn = pn+1 + n' pn + n' pn-1 pn = monic where we select the names n' and n' in order to match the expression (1.7) of Gautschi (except we have primes and he does not), so we have n' = - Bn = (rn - rn+1) = kn' - kn+1' n' = Cn = (hn/hn-1) Now suppose we apply the following scaling transformation pn(x) = qnfn(x) then the above becomes x qnfn = qn+1 fn+1 + n' qnfn + n' qn-1fn-1 x fn = (qn+1/qn) fn+1 + n' fn + n' (qn-1/qn)fn-1 Now select (qn+1/qn) = => (qn/qn-1) = to get x fn = fn+1 + n' fn + fn-1 The scaling factor qn can be determined by iterating this expression, (qn/qn-1) = = Cn = n = 1,2,3... as noted earlier Only the q ratio shows up in our work above, so we can have q0 be an arbitrary constant. We then have q12 = q02 (h1/h0) q22 = q12 (h2/h1) = q02 (h1/h0)(h2/h1) = q02 (h2/h0) q32 = q22 (h3/h2) = q02 (h2/h0) (h3/h2) = q02 (h3/h0) ..... qn2 = q02 (hn/h0) => (hn/qn2) = (h0/q02) The normalization for our monic pn is (pn, pm) = hn n.m . The scaled functions fn have this normalization (fn, fm) = (1/[qnqm]) (pn, pm) = n,m (hn/qn2) =(h0/q02) n,m How just for reference, we have p0(x) = 1 since monic, so h0 = (p0, p0) = (1,1) = a number determined by w(x) and (a,b). However, as noted above, q0 was an arbitrary number, so we can set q02 = h0 to get (fn, fm) = n,m But the orthonormalized ortho poly set is unique, so it seems to me that we must have fn = n of our previous section, up to a possible sign. That is, we could have f3 = - 3 and f2 = + 2 and things still work out OK (assuming everything real here). Let's assume for the moment no sign issues. Then we have x n = n+1 + n' n + n-1 x n = n n+1 + n n + n-1 n-1 We can then identify = n => = n-1 n' = n = n-1 So everything is consistent. All we have here is a stupid "name change". There is nothing special about the square roots appearing. The interesting fact is that when you deal with monics, it is convenient to use the n', n' type coefficients, and when the n n are expressed in terms of these, we get square roots. Comment: The monic form appears in Gautschi as pn+1 = (x-n')pn - n' pn-1. We can compare this to Ball/Beebe (27) to see that they use bn = n' and an =n' . Thus they would get x n = n+1 + bn n + n-1 with n = their Mn. This then explains the comment below eq (31) where I was wondering why the square roots appeared. _________________________________________________________________________ 4. Show how to compute the ortho polys from the recursion relation coefficients. Let's start with this simple form x n =nn+1 + n n + n-1 n-1 or n = n n+1 = (x - n)(1/n)n - (n-1/n) n-1 We know that 0 = constant, and we know -1 = 0. So away we go: n = 0 1 = (x - 0)(1/0)0 n = 1 2 = (x - 1)(1/1)1 - (0/1) 0 = [(x - 1)(1/1)(x - 0)(1/0)0 - (0/1) 0 ] n = 2 3 = (x - 2)(1/2)2 - (1/2) 1 = (x - 2)(1/2)[(x - 1)(1/1)(x - 0)(1/0)0 - (0/1) 0 ] - (1/2) [(x - 0)(1/0)0] = (x - 2)(1/2)(x - 1)(1/1)(x - 0)(1/0)0 - (x - 2)(1/2)(0/1) 0 - (1/2) [(x - 0)(1/0)0] We can see that 3 is cubic and so on, so we can see for sure how this works. Getting a closed form for the result may not be so easy. In practice, you would just compute them one at a time using the recursion and be happy with that. You can make a program compute the first 100 no problem. _________________________________________________________________________ 5. Show how to compute the recursion relation coefficients and the ortho polys at the same time, given knowledge of the w(x) and (a,b). [ This is the 1884 Stieltjes method ] We know from above that we can get things into this form: x n = nn+1 + n n + n-1 n-1 (n, m) = n,m where we can rewrite the first as n+1 = (x - n)(1/n)n - (n-1/n) n-1 First, close with (xn to get (xn, xn) = n2 + n2 + n-12 // This is Jim (3.13) which we can write as n2 = (n, x2n) - n2 - n-12 Next, close with (n to get n = (n, xn) // This is Jim (3.14) Here then is our process for computing everything at once. To get started, we always say -1 = 0 for an ortho poly system of any kind. In this case, the -1 in the recursion formula is a "don't care" so the simplest thing is to set -1 = 0. It has no effect on anything. So -1 = 0 -1 = 0 Next, 0 = constant, and we will just refer to this constant as 0 rather than make up a new symbol. We then have 0 = 0 0 = (0, x0) = 02 (1,x) 0 = (0, x20) - 02 - -12 = 02 (1,x2) - 02 (1,x2) Now we move up to the 1 level: 1(x) = (x - 0)(1/0)0 1 = (1, x1) 12 = (1, x21) - 12 - 02 And on up to the 2 level: 2(x) = (x - 1)(1/1)1 - (0/1) 0 2 = (2, x2) 22 = (2, x22) - 22 - 12 We can then continue this procedure up to any N. This method is called "The Stieltjes Procedure" and was first developed in 1884. 6a. Compute the ortho poly recursion coefficients A,B,C as integrals of the polys. Here is our recursion formula pn+1 = (Anx+Bn) pn - Cnpn-1 Consider: (pn, pn+1) = An(pn, xpn) + Bn (pn, pn) - Cn(pn,pn-1) From Lemma S3 of my Erdelyi notes we know that (pn, xpn) = kn[ rn - rn+1](pn, xn) = ( rn - rn+1) hn = - (Bn/An)hn. Thus we have 0 = An * (Bn/An)hn + Bnhn = - Bnhn + Bnhn and we have learned nothing at all! Hmmm. Well I guess you regard (xpn, pn) as an integral you can do knowing pn and w(x), then you can say - (Bn/An)hn = (xpn, pn) = w(x)dx x [pn(x)]2 // compare to Jim's (3.13) Now let's try another pk : (pn+1, pn+1) = An(pn+1, xpn) + Bn (pn+1, pn) - Cn(pn+1,pn-1) hn+1 = An (pn+1, xpn) So this tells us that hn+1(1/An)= (pn+1, xpn) And another pk: (pn-1, pn+1) = An(pn-1, xpn) + Bn (pn-1, pn) - Cn(pn-1,pn-1) 0 = An(pn-1, xpn) - Cn hn-1 Then we get hn-1 (Cn/An) = (pn-1, xpn) There is one other thing we can do of interest which is this: pn+1 = (Anx+Bn) pn - Cnpn-1 Anxpn = pn+1 - Bnpn + Cnpn+1 Then square both sides and integrate to get An2 (xpn,xpn) = hn+1 + Bn2 hn + Cn2 hn+1 so that hn+1 (1/An)2 + hn (Bn/An)2 + hn+1(Cn/An)2 = w(x)dx x2 [pn(x)]2 // compare to Jim's (3.14) So we end up with the following interesting set of relations: hn+1(1/An) = (pn+1, xpn) hn-1 (Cn/An) = (pn-1, xpn) hn (Bn/An) = - (pn, xpn) hn+1 (1/An)2 + hn (Bn/An)2 + hn+1(Cn/An)2 = (xpn,xpn) _______________________________________________________________________________ 6b. Repeat previous section for the "early Jim symmetrized recursion relation" polys. We know from notes elsewhere that we can get things into this form: x n = nn+1 + n n + n-1 n-1 (n, m) = n,m First, close with (xn to get (xn, xn) = n2 + n2 + n-12 // This is Jim (3.13) Next, close with (n to get (n, xn) = n // This is Jim (3.14) The second integral gives n, then we iterate the previous result to find the n . Therefore, this shows that we can "compute the recursion coefficients" if we know the orthonormalized polynomials n and the weight function w(x) and the interval (a,b) n = dx w(x)x [n(x)]2 n2 = dx w(x)x2 [n(x)]2 - n2 - n-12 6c. Repeat previous section for the monic normalized polys. We know from notes elsewhere that we can get things into this form for monics, x pn = pn+1 + n' pn + n' pn-1 (pn, pn) = hn First, close with (pn to get (pn, x pn) = n' hn Then close with (pn-1 to get (pn-1, x pn) = n' hn-1 Now we can manipulate the LHS as follows, due to monic-ness (pn-1, x pn) = (x pn-1, pn) = (xn + lower, pn) = (xn, pn) = ( pn + lower, pn) = (pn, pn) = hn causing this result: n' = (hn/hn-1) // this agrees with what we find below in section **** n' = (pn, x pn)/hn // this is the new result here. So we can write this as follows: hn (pn, pn) n' = (hn/hn-1) sn (pn, x pn) n' = (sn/hn) In Ball/Beebe this appears as tn (Mn, Mn) an = (tn/tn-1) sn (Mn, x Mn) bn = (sn/tn) where they use an = n' and bn = n' as indicated in their equation (27). Also, we know that pn = n since the n are orthonormal. Thus, B&B say Mn = n top of page 9.