外文原文：PADE APPROXIMATION BY RATIONAL FUNCTION 129We can apply this formula to get the polynomial approximation directly for a given function f (x), without having to resort to the Lagrange or Newtonpolynomial. Given a function, the degree of the approximate polynomial, and the left/right boundary points of the interval, the above MATLAB routine “cheby()” uses this formula to make the Chebyshev polynomial approximation.The following example illustrates that this formula gives the same approximate polynomial function as could be obtained by applying the Newton polynomial with the Chebyshev nodes.Example 3.1. Approximation by Chebyshev Polynomial. Consider the problem of finding the second-degree (N = 2) polynomial to approximate the function 2()1/(18)f x x =+. We make the following program “do_cheby.m ”, which uses the MATLAB routine “cheby()” for this job and uses Lagrange/Newton polynomial with the Chebyshev nodes to do the same job. Readers can run this programto check if the results are the same.3.4 PADE APPROXIMATION BY RATIONAL FUNCTIONPade approximation tries to approximate a function f (x) around a point xo by a rational function00,00020012002012()()()()()()1()()()M M N N M M N N Q x x p x x D x x q q x x q x x q x x d x x d x x x x --=-+-+--+-+-+-＋＋＝＋d (3.4.1)where 00(2)0()0(),'(),(),,()M N f x f x f x f x + are known.How do we find such a rational function? We write the Taylor series expansion of f (x) up to degree M + N at x = xo as130 INTERPOLATION AND CURVE FITTING(2)0()000000200020012()()()()()'()()()()2()!()()()(3.4.2)M N M NM N M N M N f x f x f x T x x f x f x x x x x x x M N a a x x a x x a x x +++++≈-=+-+-++-+=+-+-++- Assuming 0x ＝0for simplicity, we get the coefficients of ()()N M D x andQ x such that()()0()M M N N Q x T x D x +-=011011()(1)()01M N N N M N N N N N a a x a x d x d x q q x q x d x d x ++++++++-+++=+++01101()(1)()M N N N M N N N a a x a x d x d x q q x q x ++++++++=+++ (3.4.3) by solving the following equations:0010112110221122111212112212000M M M M N N M M M M M N N M M M M N N M N M N M N M N a q a a d q a a d a d q a a d a d a d q a a d a d a d a a d a d a d a a d a d a d ---+--+++-++++=⎧⎪+=⎪⎪+=⎪⎪⎪⎨⎪=⎪+=⎪⎪⎪⎪=⎩＋＋＋＋＋＋＋＋＋＋＋＋（3.4.4a ）（3.4.4b ） Here, we must first solve Eq. (3.4.4b) for N d d d ,,,21 and then substitute di ’sinto Eq. (3.4.4a) to obtain M q q q ,,,10The MATLAB routine “padeap()” implements this scheme to find the coefficient vectors of the numerator/denominator polynomial )(/)(x D x Q N M of thePade approximation for a given function f (x). Note the following things:ž The derivatives (2)()000'(),(),,()M N f x f x f x + up to order (M + N ) arecomputed numerically by using the routine “difapx()”, that will be introduced in Section 5.3.ž In order to compute the values of the Pade approximate function, we substitute 0x x - for x in ,()M N p x which has been obtained with the assumptionthat 0x ＝0.PADE APPROXIMATION BY RATIONAL FUNCTION 131Example 3.2. Pade Approximation for ()x f x e = . Let’s find the Padeapproximation3,232()()/()p x Q x D x = for ()x f x e = around 0x ＝0. We make theMATLAB program “do_pade.m ”, which uses the routine “padeap()” for thisjob and uses it again with no output argument to see the graphic results as depicted in Fig. 3.6.To confirm and support this result from the analytical point of view and to help the readers understand the internal mechanism, we perform the hand-calculation 132 INTERPOLATION AND CURVE FITTINGFigure 3.6 Pade approximation and Taylor series expansion for f(x) = ex(Example 3.2.).procedure. First, we write the Taylor series expansion at x = 0 up to degreeM + N = 5 for the given function ()xasf x e()230()11()1!23!K M Nk y k f x T x x x x x k +===++++∑ (E3.2.1) whose coefficients are0123111,1,,,26a a a a ==== (E3.2.2)We put this into Eq. (3.4.4b) with M = 3,N = 2 and solve it for di ’s to get2212431223211211221()10,02/51/61/21/24,/201/241/61/120D x d x d x a a d a d a a d a d d d d d a =++++=++=--⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦ (E3.2.3)Substituting this to Eq. (3.4.4a) yields001101221102332112111(2/5)3/51/21(2/5)1(1/20)3/201/6(1/2)(2/5)1(1/20)1/60q a q a a d q a a d a d q a a d a d ===+=+⨯-==++=+⨯-+⨯==++=+⨯-+⨯= (E3.2.4) INTERPOLATION BY CUBIC SPLINE 133 With these coefficients, we write the Pade approximate function as2333,222322()1(3/5)(3/20)(1/60)()()1(2/5)(1/20)(1/3)31220820Q x x x x p x D x x x x x x x x +++===+-++++-+ (E3.2.5)3.5 INTERPOLATION BY CUBIC SPLINEIf we use the Lagrange/Newton polynomial to interpolate a given set of N + 1data points, the polynomial is usually of degree N and so has N − 1 local extrema (maxima/minima). Thus, it will show a wild swing/oscillation (called ‘polynomial wiggle’), particularly near the ends of the whole interval as the number of datapoints increases and so the degree of the polynomial gets higher, as illustratedin Fig. 3.2. Then, how about a piecewise-linear approach, like assigning theindividual approximate polynomial to every subinterval between data points?How about just a linear interpolation —that is, connecting the data points bya straight line? It is so simple, but too short of smoothness. Even with thesecond-degree polynomial, the piecewise-quadratic curve is not smooth enoughto please our eyes, since the second-order derivatives of quadratic polynomials for ad jacent subintervals can’t be made to conform with each other. In real life, there are many cases where the continuity of second-order derivatives isdesirable. For example, it is very important to ensure the smoothness up to order 2 for interpolation needed in CAD (computer-aided design)/CAM (computer-aided manufacturing), computer graphic, and robot path/trajectory planning. That’s why we often resort to the piecewise-cubic curve constructed by the individual thirddegree polynomials assigned to each subinterval, which is called the cubic splineinterpolation. (A spline is a kind of template that architects use to draw a smooth curve between two points.)For a given set of data points {(,),0:}k k x y k N =, the cubic spline s(x)consists of N cubic polynomial ()k s x ’s assigned to each subinterval satisfying the following constraints (S0)–(S4).（S0）32,3,2,1,01()()()()(),[,],0:k k k k k k k k k k s x s x S x x S x x S x x S x x x k N +==-+-+-+∈=for （S1）,0()0:k k k k s x S y k N ===，for（S2）1,0()(),1:1k k k k k k s x s x S y k N -≡===-for（S3）1,1'()'(),1:1k k k k k s x s x S k N -≡==-for（S4）1,2''()''()2,1:1k k k k k s x s x S k N -≡==-forThese constraints (S1)–(S4) amount to a set of N + 1 + 3(N − 1) = 4N − 2linear equations having 4N coefficients of the N cubic polynomials,0,1,2,3{,,,,0:1}k k k k S S S S k N =-134 INTERPOLATION AND CURVE FITTINGTable 3.4 Boundary Conditions for a Cubic Splineas their unknowns. Two additional equations necessary for the equations to besolvable are supposed to come from the boundary conditions for the first/secondorder derivatives at the end points (x 0, y 0) and (xN, yN ) as listed in Table 3.4.Now, noting from (S1) that ,0k k S y =, we will arrange the constraints (S2)–(S4)and eliminate Sk,1, Sk,3’s to set up a set of equations with respect to the N + 1 unknowns ,2{,0:}k S k N =. In order to do so, we denote each interval width by 1k k k h x x +=- and substitute (S0) into (S4) to write1,3,2111,2,31,2,2''()62''()21()3k k k k k k k k k k k k s x S h S s x S S h S S +++++=+≡==- (3.5.1a)1,31,21,21()3k k k k S h S S ---=- (3.5.1b) We substitute these equations into (S2) with k + 1 in place of k321,31,21,11,0132,3,2,11()()()()k k k k k k k k k k k k k k k k k k k k k s x S x x S x x S x x S y S h S h S h y y ++++++=-+-+-+≡+++≡to eliminate ,3k S ’s and rewrite it as11,2,2,2,1()3k k k k k k k k k kh y y S S S h S dy h ++--++== 1,2,2,1(2)33k k k k k h S S S dy +++= （3.5.2a ） 1,21,21,11(2)33k k k k k h S S S dy ----++= （3.5.2b ） We also substitute Eq. (3.5.1b) into (S3)211,311,211,1,1'()3k k k k k k k k s x S h S h S S ------=++≡INTERPOLATION BY CUBIC SPLINE 135to write,11,11,21,211,21,21,2()2()k k k k k k k k k k S S h S S h S h S S --------=-+=+ (3.5.3)In order to use this for eliminating Sk,1 from Eq. (3.5.2), we subtract (3.5.2b) from (3.5.2a) to write1,2,21,21,2,11,11(2)(2)3(2)3()k k k k k k k k k k h S S h S S S S dy dy +----+-++-=-and then substitute Eq. (3.5.3) into this to write1,2,21,21,21,21,2111,21,21,21(2)(2)3()3()()3()k k k k k k k k k k k k k k k k k k k k h S S h S S h S S dy dy h S h h S h S dy dy +--------+-+-+++=-+++=- (3.5.4)for k = 1 : N − 1 Since these are N − 1 equations with respect to N + 1 unknowns ,2{,0:}k S k N =, we need two more equations from the boundary conditions to be given aslisted in Table 3.4.How do we convert the boundary condition into equations? In the case where the first-order derivatives on the two boundary points are given as (i) in Table 3.4, we write Eq. (3.5.2a) for k = 0 as01,20,20,1000,201,200,1(2)33,233h S S S dy h S h S dy S ++=+=- (3.5.5a )We also write Eq. (3.5.2b) for k = N as1,21,21,11(2)33N N N N N h S S S dy ----++=and substitute (3.5.3)(k = N ) into this to write01,20,20,10(2)33h S S S dy ++=1,21,21,11,21,2111,21,2,11(2)33(2)323()N N N N N N N N N N N N N N h S S S h S S dy h S h S S dy ----------++-+=+=- (3.5.5b)Equations (3.5.5a) and (3.5.5b) are two additional equations that we need to solve Eq. (3.5.4) and that’s it. In the case where the second -order derivatives on the two boundary points are given as (ii) in Table 3.4, 0,2S and ,2N S are directly known from the boundary conditions as0,200''()/2S s x =,,2''()/2N N N S s x = (3.5.6) 136 INTERPOLATION AND CURVE FITTINGand, subsequently, we have just N − 1 unknowns. In the case where the secondorderderivatives on the two boundary points are given as (iii) in Table 3.4we can instantly convert these into two equations with respect to S 0,2 and SN,2 as Finally, we combine the two equations (3.5.5a) and (3.5.5b) with Eq. (3.5.4) to write it in the matrix –vector form as0,200,1001,21000111,2122211,2,11113()203()2()003()2()3()02N N N N N N N N N N N N S dy S h h S dy dy h h h h S dy dy h h h h S S dy h h -----------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦ (3.5.8) After solving this system of equation for {Sk,2, k = 0 : N }, we substitute theminto (S1), (3.5.2), and (3.5.1) to get the other coefficients of the cubic spline asThe MATLAB routine “cspline()” constructs Eq.(3.5.8), solves it to get the cubic spline coefficients for given x, y coordinates of the data points and the boundary conditions, uses the mkpp() routine to get the piecewise polynomialexpression, and then uses the ppval() routine to obtain the value(s) of the piecewise polynomial function for xi —that is, the interpolation over xi. The type ofthe boundary condition is supposed to be specified by the third input argument KC. In the case where the boundary condition is given as (i)/(ii) in Table 3.4,the input argument KC should be set to 1/2 and the fourth and fifth input arguments must be the first/second derivatives at the end points. In the case wherethe boundary condition is given as extrapolated like (iii) in Table 3.4, the input argument KC should be set to 3 and the fourth and fifth input arguments do not need to be fed.INTERPOLATION BY CUBIC SPLINE 137(cf) See Problem 1.11 for the usages of the MATLAB routines mkpp() and ppval(). Example 3.3. Cubic Spline. Consider the problem of finding the cubic spline interpolation for the N + 1 = 4 data pointssubject to the boundary condition(E3.3.2)With the subinterval widths on the x-axis and the first divided differences as138INTERPOLATION AND CURVE FITTINGwe write Eq. (3.5.8) as0,200,11,2102,2213,23,113()210033()141063()014163()00123S dy S S dy dy S dy dy S S dy --⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥-⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦⎣⎦(E3.3.4)Then we solve this equation to get and substitute this into Eq. (3.5.9) to obtain0,01,02,000,102,20,211,112,21,222,123,21,21,20,20,301,21,21,313,22,30,1,4( 3.3.6)(2)2( 3.3.7)3(2)2( 3.3.7)3(2)2( 3.3.7)32( 3.3.8)32( 3.3.8)3S S S E h S dy S S E a hS dy S S E b h S dy S S E c S S S E a h S S S E b h S S ====-+==-+==-+=-==-==--=2,222( 3.3.8)3S E c h =Figure 3.7 Cubic splines for Example 3.3.Finally, we can write the cubic spline equations collectively from (S0) asWe make and run the program “do_csplines.m”, which uses the routine“cspline()” to compute the cubic spline coefficientsand obtain the value(s) of the cubicspline function for i x(i.e., theinterpolation over xi) and then plots the result as depicted in Fig. 3.7. We also compare this result with that obtained by using the MATLAB built-in function “spline(x,y,xi)”, which works with the boundary condition of type (i) for the second input argument given as [dy0 y dyN], and with the boundary conditionof type (iii) for the same lengths of x and y.3.6 HERMITE INTERPOLATING POLYNOMIALIn some cases, we need to find the polynomial function that not only passes through the given points, but also has the specified derivatives at every data point. We call such a polynomial the Hermite interpolating polynomial or the osculating polynomial.140 INTERPOLATION AND CURVE FITTINGFor simplicity, we consider a third-order polynomialmatching just two points and having the specified first derivativesat the points. We can obtain the four coefficients by solvingAs an alternative, we approximate the specified derivatives at the data points by their differencesand find the Lagrange/Newton polynomial matching the four pointsThe MATLAB routine “hermit()” constructs Eq. (3.6.2) and solves it to getthe Hermite interpolating polynomial coefficients for a single interval given thetwo end points and the derivatives at them as the input arguments. The next routine “hermits()” uses “hermit()” to get the Hermite coefficients for a setof multiple subintervals.Example 3.4. Hermite Interpolating Polynomial. Consider the problem of finding the polynomial interpolation for the N + 1 = 4 data pointsTWO-DI MENSIONAL INTERPOLATION 141subject to the conditionsFor this problem, we only have to type the following statements in the MATLAB command window.3.7 TWO-DIMENSIONAL INTERPOLATIONIn this section we deal with only the simplest way of two-dimensional interpolation—that is, a generalization of piecewise linear interpolation calledFigure 3.8 A two-dimensional interpolation using Zi = interp2() on the grid points [Xi,Yi] generated by the meshgrid() command.142 INTERPOLATION AND CURVE FITTINGthe bilinear interpolation. The bilinear interpolation for a point (x, y ) on the rectangularsub-region having and as its left-upper/right-lower corner points is described by the following formula.11111,1111,11,11,11(,)(,)(,)1{()()()()()()()()()()},,n n n n n n n n m n m n m m n n m n m n m n m n m n m n m m n ny y y y z x y z x y z x y y y y y x x y y z x x y y x x y y z x x y y z x x y y z x x x y y y -----------------=+--=----+--+--+--≤≤≤≤for (3.7.2)This formula is cast into the MATLAB routine “intrp2()”, which is so named in order to distinguish it from the MATLAB built-in routine “interp2()”. Note that in reference to Fig. 3.8, the given values of data at grid points (x(m),y(n)) and the interpolated values for intermediate points (xi(m),yi(n)) are stored in Z(n,m) and Zi(n,m), respectively.CURVE FITTING 143Example 3.5. Two-Dimensional Bilinear Interpolation. We consider interpolating the sample values of a functionfor the 5 × 5 grid over the 21 × 21 grid on the domainWe make the MATLA B program “do_interp2.m”, which uses the routine“intrp2()” to do this job, compares its function with that of the MATLABbuilt-in routine “interp2()”, and computes a kind of relative error to estimatehow close the interpolated values are to the original values. The graphic resultsof running this program are depicted in Fig. 3.9, which shows that we obtaineda reasonable approximation with the error of 2.6% from less than 1/16 of the original data. It is implied that the sampling may be a simple data compression method, as long as the interpolated data are little impaired.3.8 CURVE FITTINGWhen many sample data pairs {(xk, yk), k = 0 : M} are available, we often need to grasp the relationship between the two variables or to describe the trend of the(a) True function (b) The function over (c) Bilinear interpolation Figure 3.9 Two-dimensional interpolation (Example 3.5).144 INTERPOLATION AND CURVE FITTINGsample grid data, hopefully in a form of function y = f (x). But, as mentioned in Remark 3.1,the polynomial approach meets with the polynomial wiggle and/or Runge phenomenon,which makes it not attractive for approximation purpose. Although thecubic spline approach may be a roundabout toward the smoothness as explainedin Section 3.5, it has too many parameters and so does not seem to be an efficient way of describing the relationship or the trend, since every subintervalneeds four coefficients. What other choices do we have? Noting that many data are susceptible to some error, we don’t ha ve to try to find a function passing exactly through every point. Instead of pursuing the exact matching at every data point, we look for an approximate function (not necessarily a polynomial) that describes the data points as a whole with the smallest error in some sense, whichis called the curve fitting.As a reasonable means, we consider the least-squares (LS) approach to minimizing the sum of squared errors, where the error is described by the verticaldistance to the curve from the data points. We will look over various types of fitting functions in this section.3.8.1 Straight Line Fit: A Polynomial Function of First DegreeIf there is some theoretical basis on which we believe the relationship betweenthe two variables to bewe should set up the following system of equations from the collection of many experimental data:Noting that this apparently corresponds to the overdetermined case mentionedin Section 2.1.3, we resort to the least-squares (LS) solution (2.1.10)which minimizes the objective functionCURVE FITTING 145Sometimes we have the information about the error bounds of the data, and it is reasonable to differentiate the data by weighing more/less each one according toits accuracy/reliability. This policy can be implemented by the weighted leastsquares (WLS) solutionwhich minimizes the weighted objective functionIf the weighting matrix is, then we can write the WLS solution (3.8.5) aswhereOne may use the MATLAB built-in routine “lscov(A,y,V)” to obtain thisWLS solution.3.8.2 Polynomial Curve Fit: A Polynomial Function of Higher DegreeIf there is no reason to limit the degree of fitting polynomial to one, then we may increase the degree of fitting polynomial to, say, N in expectation of decreasing the error. Still, we can use Eq. (3.8.4) or (3.8.6), but with different definitions of A and θ asThe MATLAB routine “polyfits()” performs the WLS or LS scheme tofind the coefficients of a polynomial fitting a given set of data points, depending on whether or not a vector (r) having the diagonal elements of the weighting matrix W is given as the fourth or fifth input argument. Note that in the case of a diagonal weighting matrix W, the WLS solution conforms to the LS solution with each row of the information matrix A and the data vector y multiplied by the corresponding element of the weighting matrix W. Let us see the following examples for its usage:146 INTERPOLATION AND CURVE FITTINGExample 3.6. Polynomial Curve Fit by LS (Least Squares). Suppose we havein two columns and we must fit these data into polynomials of degree 1, 3, 5, and 7.We make the MATLAB program “do_polyfit.m”, which uses the routine “polyfits()” to do this job and plot the results together with the given dataCURVE FITTING 147Figure 3.10 Polynomial curve fitting by the LS (Least-Squares) method.points as depicted in Fig. 3.10. We can observe the polynomial wiggle that the oscillation of the fitting curve between the data points becomes more pronounced with higher degree.Example 3.7. Curve Fitting by WLS (Weighted Least Squares). Most experimental data have some absolute and/or relative error bounds that are not uniform forall data. If we know the error bounds for each data, we may give each data a weight inversely proportional to the size of its error bound when extracting valuable information from the data. The WLS solution (3.8.7) enables us to reflect sucha weighting strategy on estimating data trends. Consider the following two cases. (a) Suppose there are two gauges A and B with the same function, but different absolute error bounds ±0.2 and ±1.0, respectively. We used themto get the input-output data pair (xm,ym) as{(1, 0.0831), (3, 0.9290), (5, 2.4932), (7, 4.9292), (9, 7.9605)} from gauge A {(2, 0.9536), (4, 2.4836), (6, 3.4173), (8, 6.3903), (10, 10.2443)} from gauge B Let the fitting function be a second-degree polynomial function148 INTERPOLATION AND CURVE FITTING(a) Fitting to a polynomial (b) Fitting toFigure 3.11 LS curve fitting and WLS curve fitting for Example 3.7.To find the parameters a2, a1, and a0, we write the MATLAB program“do_wlse1.m”, which uses the routine “polyfits()” twice, once without weighting coefficients and once with weighting coefficients. The resultsare depicted in Fig. 3.11a, which shows that the WLS curve fitting triesto be closer to the data points with smaller error bound, while the LScurve fitting weights all data points equally, which may result in larger deviations from data points with small error bounds.(b) Suppose we use one gauge that has relative error bound ±40[%] for measuring the output y for the input values x = [1, 3, 5, . . . , 19] and sothe size of error bound of each output data is proportional to the magnitudeof the output. We used it to get the input–output data pair (xm,ym) as{(1, 4.7334), (3, 2.1873), (5, 3.0067), (7, 1.4273), (9, 1.7787)(11, 1.2301), (13, 1.6052), (15, 1.5353), (17, 1.3985), (19, 2.0211)}Let the fitting function be an exponential functiony = axb (E3.7.2)To find the parameters a and b, we make the MATLAB program“do_wlse2.m”, which uses the routine “curve_fit()” without the weighting coefficients one time and with the weighting coefficients another time.The results depicted in Fig. 3.11b shows that the WLS curve fitting tries toget closer to the data points with smaller |y|, while the LS curve fitting pays equal respect to all data points, which may result in larger deviation fromdata points with small |y|. Note that the MATLAB routine “curve_fit()”appears in Problem 3.11, which implements all of the schemes listed inTable 3.5 with the LS/WLS solution.CURVE FITTING 149(cf) Note that the objective of the WLS scheme is to put greater emphasis on more reliable data.中文译文：插值与拟合我们能够运用这个公式直接得到给定函数f(x)的近似多项式估计，没有必要采取拉格朗日或牛顿多项式。
