Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 52
Текст из файла (страница 52)
These can be found by the methods of §5.8, but not by economization ofseries. There is slightly less to economization of series than meets the eye.CITED REFERENCES AND FURTHER READING:Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 12.Arfken, G. 1970, Mathematical Methods for Physicists, 2nd ed. (New York: Academic Press),p. 631. [1]5.12 Padé ApproximantsA Padé approximant, so called, is that rational function (of a specified order) whosepower series expansion agrees with a given power series to the highest possible order. Ifthe rational function isMR(x) ≡ak xkk=01+Nk=1(5.12.1)bk xk5.12 Padé Approximants201then R(x) is said to be a Padé approximant to the seriesf (x) ≡∞c k xk(5.12.2)k=0ifR(0) = f (0)and alsodkdkR(x)=f(x),kkdxdxx=0x=0(5.12.3)k = 1, 2, .
. . , M + N(5.12.4)Equations (5.12.3) and (5.12.4) furnish M + N + 1 equations for the unknowns a0 , . . . , aMand b1 , . . . , bN . The easiest way to see what these equations are is to equate (5.12.1) and(5.12.2), multiply both by the denominator of equation (5.12.1), and equate all powers ofx that have either a’s or b’s in their coefficients. If we consider only the special case ofa diagonal rational approximation, M = N (cf. §3.2), then we have a0 = c0 , with theremaining a’s and b’s satisfyingNbm cN −m+k = −cN +k ,k = 1, . .
. , N(5.12.5)k = 1, . . . , N(5.12.6)m=1kbm ck−m = ak ,m=0(note, in equation 5.12.1, that b0 = 1). To solve these, start with equations (5.12.5), whichare a set of linear equations for all the unknown b’s. Although the set is in the form of aToeplitz matrix (compare equation 2.8.8), experience shows that the equations are frequentlyclose to singular, so that one should not solve them by the methods of §2.8, but rather byfull LU decomposition.
Additionally, it is a good idea to refine the solution by iterativeimprovement (routine mprove in §2.5) [1].Once the b’s are known, then equation (5.12.6) gives an explicit formula for the unknowna’s, completing the solution.Padé approximants are typically used when there is some unknown underlying functionf (x). We suppose that you are able somehow to compute, perhaps by laborious analyticexpansions, the values of f (x) and a few of its derivatives at x = 0: f (0), f (0), f (0),and so on. These are of course the first few coefficients in the power series expansion off (x); but they are not necessarily getting small, and you have no idea where (or whether)the power series is convergent.By contrast with techniques like Chebyshev approximation (§5.8) or economizationof power series (§5.11) that only condense the information that you already know about afunction, Padé approximants can give you genuinely new information about your function’svalues.
It is sometimes quite mysterious how well this can work. (Like other mysteries inmathematics, it relates to analyticity.) An example will illustrate.Imagine that, by extraordinary labors, you have ground out the first five terms in thepower series expansion of an unknown function f (x),f (x) ≈ 2 +11 249 3175 4x+x −x +x +···981874878732(5.12.7)(It is not really necessary that you know the coefficients in exact rational form — numericalvalues are just as good. We here write them as rationals to give you the impression thatthey derive from some side analytic calculation.) Equation (5.12.7) is plotted as the curvelabeled “power series” in Figure 5.12.1. One sees that for x >∼ 4 it is dominated by itslargest, quartic, term.We now take the five coefficients in equation (5.12.7) and run them through the routinepade listed below.
It returns five rational coefficients, three a’s and two b’s, for use in equation(5.12.1) with M = N = 2. The curve in the figure labeled “Padé” plots the resulting rationalfunction. Note that both solid curves derive from the same five original coefficient values.202Chapter 5.Evaluation of Functions10f(x) = [7 + (1 + x)4/3]1/386f (x)power series (5 terms)4Padé (5 coefficients)exact200246810xFigure 5.12.1. The five-term power series expansion and the derived five-coefficient Padé approximantfor a sample function f (x). The full power series converges only for x < 1. Note that the Padéapproximant maintains accuracy far outside the radius of convergence of the series.To evaluate the results, we need Deus ex machina (a useful fellow, when he is available)to tell us that equation (5.12.7) is in fact the power series expansion of the functionf (x) = [7 + (1 + x)4/3 ]1/3(5.12.8)which is plotted as the dotted curve in the figure.
This function has a branch point at x = −1,so its power series is convergent only in the range −1 < x < 1. In most of the rangeshown in the figure, the series is divergent, and the value of its truncation to five terms israther meaningless. Nevertheless, those five terms, converted to a Padé approximant, give aremarkably good representation of the function up to at least x ∼ 10.Why does this work? Are there not other functions with the same first five terms intheir power series, but completely different behavior in the range (say) 2 < x < 10? Indeedthere are.
Padé approximation has the uncanny knack of picking the function you had inmind from among all the possibilities. Except when it doesn’t! That is the downside ofPadé approximation: it is uncontrolled. There is, in general, no way to tell how accurateit is, or how far out in x it can usefully be extended.
It is a powerful, but in the end stillmysterious, technique.Here is the routine that gets a’s and b’s from your c’s. Note that the routine is specializedto the case M = N , and also that, on output, the rational coefficients are arranged in a formatfor use with the evaluation routine ratval (§5.3). (Also for consistency with that routine,the array of c’s is passed in double precision.)#include <math.h>#include "nrutil.h"#define BIG 1.0e30void pade(double cof[], int n, float *resid)Given cof[0..2*n], the leading terms in the power series expansion of a function, solve thelinear Padé equations to return the coefficients of a diagonal rational function approximation tothe same function, namely (cof[0] + cof[1]x + · · · + cof[n]xN )/(1 + cof[n+1]x + · · · +5.12 Padé Approximants203cof[2*n]xN ). The value resid is the norm of the residual vector; a small value indicates awell-converged solution.
Note that cof is double precision for consistency with ratval.{void lubksb(float **a, int n, int *indx, float b[]);void ludcmp(float **a, int n, int *indx, float *d);void mprove(float **a, float **alud, int n, int indx[], float b[],float x[]);int j,k,*indx;float d,rr,rrold,sum,**q,**qlu,*x,*y,*z;indx=ivector(1,n);q=matrix(1,n,1,n);qlu=matrix(1,n,1,n);x=vector(1,n);y=vector(1,n);z=vector(1,n);for (j=1;j<=n;j++) {Set up matrix for solving.y[j]=x[j]=cof[n+j];for (k=1;k<=n;k++) {q[j][k]=cof[j-k+n];qlu[j][k]=q[j][k];}}ludcmp(qlu,n,indx,&d);Solve by LU decomposition and backsubstitulubksb(qlu,n,indx,x);tion.rr=BIG;do {Important to use iterative improvement, sincerrold=rr;the Padé equations tend to be ill-conditioned.for (j=1;j<=n;j++) z[j]=x[j];mprove(q,qlu,n,indx,y,x);for (rr=0.0,j=1;j<=n;j++)Calculate residual.rr += SQR(z[j]-x[j]);} while (rr < rrold);If it is no longer improving, call it quits.*resid=sqrt(rr);for (k=1;k<=n;k++) {Calculate the remaining coefficients.for (sum=cof[k],j=1;j<=k;j++) sum -= x[j]*cof[k-j];y[k]=sum;}Copy answers to output.for (j=1;j<=n;j++) {cof[j]=y[j];cof[j+n] = -x[j];}free_vector(z,1,n);free_vector(y,1,n);free_vector(x,1,n);free_matrix(qlu,1,n,1,n);free_matrix(q,1,n,1,n);free_ivector(indx,1,n);}CITED REFERENCES AND FURTHER READING:Ralston, A.
and Wilf, H.S. 1960, Mathematical Methods for Digital Computers (New York: Wiley),p. 14.Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: NorthHolland), Chapter 2.Graves-Morris, P.R. 1979, in Padé Approximation and Its Applications, Lecture Notes in Mathematics, vol. 765, L. Wuytack, ed.
(Berlin: Springer-Verlag). [1]204Chapter 5.Evaluation of Functions5.13 Rational Chebyshev ApproximationIn §5.8 and §5.10 we learned how to find good polynomial approximations to a givenfunction f (x) in a given interval a ≤ x ≤ b. Here, we want to generalize the task to findgood approximations that are rational functions (see §5.3). The reason for doing so is that,for some functions and some intervals, the optimal rational function approximation is ableto achieve substantially higher accuracy than the optimal polynomial approximation with thesame number of coefficients.
This must be weighed against the fact that finding a rationalfunction approximation is not as straightforward as finding a polynomial approximation,which, as we saw, could be done elegantly via Chebyshev polynomials.Let the desired rational function R(x) have numerator of degree m and denominatorof degree k.
Then we haveR(x) ≡p0 + p1 x + · · · + pm xm≈ f (x)1 + q1 x + · · · + qk xkfor a ≤ x ≤ b(5.13.1)The unknown quantities that we need to find are p0 , . . . , pm and q1 , . . . , qk , that is, m + k + 1quantities in all. Let r(x) denote the deviation of R(x) from f (x), and let r denote itsmaximum absolute value,r(x) ≡ R(x) − f (x)r ≡ max |r(x)|a≤x≤b(5.13.2)The ideal minimax solution would be that choice of p’s and q’s that minimizes r. Obviouslythere is some minimax solution, since r is bounded below by zero. How can we find it, ora reasonable approximation to it?A first hint is furnished by the following fundamental theorem: If R(x) is nondegenerate(has no common polynomial factors in numerator and denominator), then there is a uniquechoice of p’s and q’s that minimizes r; for this choice, r(x) has m + k + 2 extrema ina ≤ x ≤ b, all of magnitude r and with alternating sign.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
