Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 47
Текст из файла (страница 47)
This is the case, e.g., for (5.5.2) in the direction of increasingn, when x < n. You cannot generate Bessel functions of high n by forwardrecurrence on (5.5.2).To state things a bit more formally, iffn /gn → 0asn→∞(5.5.10)then fn is called the minimal solution of the recurrence relation (5.5.9).
Nonminimalsolutions like gn are called dominant solutions. The minimal solution is unique, if itexists, but dominant solutions are not — you can add an arbitrary multiple of fn toa given gn . You can evaluate any dominant solution by forward recurrence, but notthe minimal solution. (Unfortunately it is sometimes the one you want.)Abramowitz and Stegun (in their Introduction) [1] give a list of recurrences thatare stable in the increasing or decreasing directions. That list does not contain all180Chapter 5.Evaluation of Functionspossible formulas, of course.
Given a recurrence relation for some function fn (x)you can test it yourself with about five minutes of (human) labor: For a fixed xin your range of interest, start the recurrence not with true values of fj (x) andfj+1 (x), but (first) with the values 1 and 0, respectively, and then (second) with0 and 1, respectively. Generate 10 or 20 terms of the recursive sequences in thedirection that you want to go (increasing or decreasing from j), for each of the twostarting conditions. Look at the difference between the corresponding members ofthe two sequences. If the differences stay of order unity (absolute value less than10, say), then the recurrence is stable. If they increase slowly, then the recurrencemay be mildly unstable but quite tolerably so. If they increase catastrophically,then there is an exponentially growing solution of the recurrence.
If you knowthat the function that you want actually corresponds to the growing solution, thenyou can keep the recurrence formula anyway e.g., the case of the Bessel functionYn (x) for increasing n, see §6.5; if you don’t know which solution your functioncorresponds to, you must at this point reject the recurrence formula. Notice thatyou can do this test before you go to the trouble of finding a numerical method forcomputing the two starting functions fj (x) and fj+1 (x): stability is a property ofthe recurrence, not of the starting values.An alternative heuristic procedure for testing stability is to replace the recurrence relation by a similar one that is linear with constant coefficients. For example,the relation (5.5.2) becomesyn+1 − 2γyn + yn−1 = 0(5.5.11)where γ ≡ n/x is treated as a constant. You solve such recurrence relationsby trying solutions of the form yn = an .
Substituting into the above recurrence gives(ora = γ ± γ2 − 1(5.5.12)a2 − 2γa + 1 = 0The recurrence is stable if |a| ≤ 1 for all solutions a. This holds (as you can verify)if |γ| ≤ 1 or n ≤ x. The recurrence (5.5.2) thus cannot be used, starting with J0 (x)and J1 (x), to compute Jn (x) for large n.Possibly you would at this point like the security of some real theorems onthis subject (although we ourselves always follow one of the heuristic procedures).Here are two theorems, due to Perron [2]:Theorem A. If in (5.5.9) an ∼ anα , bn ∼ bnβ as n → ∞, and β < 2α, thengn+1 /gn ∼ −anα ,fn+1 /fn ∼ −(b/a)nβ−α(5.5.13)and fn is the minimal solution to (5.5.9).Theorem B. Under the same conditions as Theorem A, but with β = 2α,consider the characteristic polynomialt2 + at + b = 0(5.5.14)If the roots t1 and t2 of (5.5.14) have distinct moduli, |t1 | > |t2 | say, thengn+1 /gn ∼ t1 nα ,fn+1 /fn ∼ t2 nα(5.5.15)5.5 Recurrence Relations and Clenshaw’s Recurrence Formula181and fn is again the minimal solution to (5.5.9).
Cases other than those in thesetwo theorems are inconclusive for the existence of minimal solutions. (For moreon the stability of recurrences, see [3].)How do you proceed if the solution that you desire is the minimal solution?The answer lies in that old aphorism, that every cloud has a silver lining: If arecurrence relation is catastrophically unstable in one direction, then that (undesired)solution will decrease very rapidly in the reverse direction. This means that youcan start with any seed values for the consecutive fj and fj+1 and (when you havegone enough steps in the stable direction) you will converge to the sequence offunctions that you want, times an unknown normalization factor.
If there is someother way to normalize the sequence (e.g., by a formula for the sum of the fn ’s),then this can be a practical means of function evaluation. The method is calledMiller’s algorithm. An example often given [1,4] uses equation (5.5.2) in just thisway, along with the normalization formula1 = J0 (x) + 2J2 (x) + 2J4 (x) + 2J6 (x) + · · ·(5.5.16)Incidentally, there is an important relation between three-term recurrencerelations and continued fractions. Rewrite the recurrence relation (5.5.9) asbnyn=−yn−1an + yn+1 /yn(5.5.17)Iterating this equation, starting with n, givesynyn−1=−bnbn+1···an − an+1 −(5.5.18)Pincherle’s Theorem [2] tells us that (5.5.18) converges if and only if (5.5.9) has aminimal solution fn , in which case it converges to fn /fn−1 .
This result, usually forthe case n = 1 and combined with some way to determine f0 , underlies many of thepractical methods for computing special functions that we give in the next chapter.Clenshaw’s Recurrence FormulaClenshaw’s recurrence formula [5] is an elegant and efficient way to evaluate asum of coefficients times functions that obey a recurrence formula, e.g.,f(θ) =Nck cos kθorf(x) =k=0Nck Pk (x)k=0Here is how it works: Suppose that the desired sum isf(x) =Nck Fk (x)(5.5.19)k=0and that Fk obeys the recurrence relationFn+1 (x) = α(n, x)Fn (x) + β(n, x)Fn−1 (x)(5.5.20)182Chapter 5.Evaluation of Functionsfor some functions α(n, x) and β(n, x).
Now define the quantities yk (k =N, N − 1, . . . , 1) by the following recurrence:yN+2 = yN+1 = 0yk = α(k, x)yk+1 + β(k + 1, x)yk+2 + ck(k = N, N − 1, . . . , 1)(5.5.21)If you solve equation (5.5.21) for ck on the left, and then write out explicitly thesum (5.5.19), it will look (in part) like this:f(x) = · · ·+ [y8 − α(8, x)y9 − β(9, x)y10 ]F8 (x)+ [y7 − α(7, x)y8 − β(8, x)y9 ]F7 (x)+ [y6 − α(6, x)y7 − β(7, x)y8 ]F6 (x)+ [y5 − α(5, x)y6 − β(6, x)y7 ]F5 (x)(5.5.22)+···+ [y2 − α(2, x)y3 − β(3, x)y4 ]F2 (x)+ [y1 − α(1, x)y2 − β(2, x)y3 ]F1 (x)+ [c0 + β(1, x)y2 − β(1, x)y2 ]F0 (x)Notice that we have added and subtracted β(1, x)y2 in the last line. If you examinethe terms containing a factor of y8 in (5.5.22), you will find that they sum to zero asa consequence of the recurrence relation (5.5.20); similarly all the other yk ’s downthrough y2 .
The only surviving terms in (5.5.22) aref(x) = β(1, x)F0 (x)y2 + F1 (x)y1 + F0 (x)c0(5.5.23)Equations (5.5.21) and (5.5.23) are Clenshaw’s recurrence formula for doing the sum(5.5.19): You make one pass down through the yk ’s using (5.5.21); when you havereached y2 and y1 you apply (5.5.23) to get the desired answer.Clenshaw’s recurrence as written above incorporates the coefficients ck in adownward order, with k decreasing. At each stage, the effect of all previous ck ’sis “remembered” as two coefficients which multiply the functions Fk+1 and Fk(ultimately F0 and F1 ).
If the functions Fk are small when k is large, and if thecoefficients ck are small when k is small, then the sum can be dominated by smallFk ’s. In this case the remembered coefficients will involve a delicate cancellationand there can be a catastrophic loss of significance. An example would be to sumthe trivial seriesJ15 (1) = 0 × J0 (1) + 0 × J1 (1) + . . .
+ 0 × J14 (1) + 1 × J15 (1)(5.5.24)Here J15 , which is tiny, ends up represented as a canceling linear combination ofJ0 and J1 , which are of order unity.5.6 Quadratic and Cubic Equations183The solution in such cases is to use an alternative Clenshaw recurrence thatincorporates ck ’s in an upward direction. The relevant equations arey−2 = y−1 = 01[yk−2 − α(k, x)yk−1 − ck ],yk =β(k + 1, x)(k = 0, 1, . . . , N − 1)f(x) = cN FN (x) − β(N, x)FN−1 (x)yN−1 − FN (x)yN−2(5.5.25)(5.5.26)(5.5.27)The rare case where equations (5.5.25)–(5.5.27) should be used instead ofequations (5.5.21) and (5.5.23) can be detected automatically by testing whetherthe operands in the first sum in (5.5.23) are opposite in sign and nearly equal inmagnitude.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
