计算机科学与技术专业外文翻译--插值与拟合

外文原文：

PADE APPROXIMATION BY RATIONAL FUNCTION 129

We can apply this formula to get the polynomial approximation directly for a given function f (x), without having to resort to the Lagrange or Newton

polynomial. 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 program

to check if the results are the same.

3.4 PADE APPROXIMATION BY RATIONAL FUNCTION

Pade approximation tries to approximate a function f (x) around a point xo by a rational function

00

,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 as

130 INTERPOLATION AND CURVE FITTING

(2)0()000000200020012()()()()()'()()()()2()!()()()(3.4.2)

M N M N

M 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:

001

01121102211221

11212112212000

M 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 ’s

into Eq. (3.4.4a) to obtain M q q q ,,,10

The MATLAB routine “padeap()” implements this scheme to find the coefficient vectors of the numerator/denominator polynomial )(/)(x D x Q N M of the

Pade 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 ) are

computed 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 assumption

that 0x ＝0.

PADE APPROXIMATION BY RATIONAL FUNCTION 131

Example 3.2. Pade Approximation for ()x f x e = . Let’s find the Pade

approximation

3,232()()/()p x Q x D x = for ()x f x e = around 0x ＝0. We make the

MATLAB program “do_pade.m ”, which uses the routine “padeap()” for this

job 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 FITTING

Figure 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 degree

M + N = 5 for the given function ()x

as

f x e

()230()11()1!23!K M N

k y k f x T x x x x x k +===++++∑ (E3.2.1) whose coefficients are

0123111,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 get

2

212431223211211221()10,0

2/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) yields

0011012211023321121

11(2/5)3/5

1/21(2/5)1(1/20)3/20

1/6(1/2)(2/5)1(1/20)1/60

q 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 as

23

33,22

2322()1(3/5)(3/20)(1/60)()()1(2/5)(1/20)(1/3)31220

820Q x x x x p x D x x x x x x x x +++===+-++++-+ (E3.2.5)

3.5 INTERPOLATION BY CUBIC SPLINE

If we use the Lagrange/Newton polynomial to interpolate a given set of N + 1

data 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 data

points increases and so the degree of the polynomial gets higher, as illustrated

in Fig. 3.2. Then, how about a piecewise-linear approach, like assigning the

individual approximate polynomial to every subinterval between data points?

How about just a linear interpolation —that is, connecting the data points by

a straight line? It is so simple, but too short of smoothness. Even with the

second-degree polynomial, the piecewise-quadratic curve is not smooth enough

to 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 is

desirable. 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 spline

interpolation. (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 -≡==-for

These constraints (S1)–(S4) amount to a set of N + 1 + 3(N − 1) = 4N − 2

linear 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 FITTING

Table 3.4 Boundary Conditions for a Cubic Spline

as their unknowns. Two additional equations necessary for the equations to be

solvable 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 write

1,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()3

k k k k S h S S ---=- (3.5.1b) We substitute these equations into (S2) with k + 1 in place of k

321,31,21,11,01

32,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 as

11,2,2,2,1()3k k k k k k k k k k

h 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 135

to 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 write

1,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 write

1,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 as

listed 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 as

01,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 as

1,21,21,11(2)33N N N N N h S S S dy ----++=

and substitute (3.5.3)(k = N ) into this to write

01,20,20,10(2)33h S S S dy ++=

1,21,21,11,21,21

11,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 as

0,200''()/2S s x =,,2''()/2N N N S s x = (3.5.6) 136 INTERPOLATION AND CURVE FITTING

and, subsequently, we have just N − 1 unknowns. In the case where the secondorder

derivatives on the two boundary points are given as (iii) in Table 3.4

we 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 as

0,200,1001,2100

0111,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 them

into (S1), (3.5.2), and (3.5.1) to get the other coefficients of the cubic spline as

The 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 polynomial

expression, 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 of

the 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 where

the 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 points

subject to the boundary condition

(E3.3.2)

With the subinterval widths on the x

-axis and the first divided differences as

138

INTERPOLATION AND CURVE FITTING

we write Eq. (3.5.8) as

0,200,11,2102,2213,23,113()

210033()141063()0

14163()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 obtain

0,01,02,00

0,102,20,21

1,112,21,22

2,123,21,21,20,2

0,30

1,21,2

1,31

3,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 h

S 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,2

22( 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) as

We make and run the program “do_csplines.m”, which uses the routine

“cspline()” to compute the cubic spline coefficients

and obtain the value(s) of the cubic

spline function for i x(i.e., the

interpolation 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 condition

of type (iii) for the same lengths of x and y.

3.6 HERMITE INTERPOLATING POLYNOMIAL

In 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 FITTING

For simplicity, we consider a third-order polynomial

matching just two points and having the specified first derivatives

at the points. We can obtain the four coefficients by solving

As an alternative, we approximate the specified derivatives at the data points by their differences

and find the Lagrange/Newton polynomial matching the four points

The MATLAB routine “hermit()” constructs Eq. (3.6.2) and solves it to get

the Hermite interpolating polynomial coefficients for a single interval given the

two end points and the derivatives at them as the input arguments. The next routine “hermits()” uses “hermit()” to get the Hermite coefficients for a set

of multiple subintervals.

Example 3.4. Hermite Interpolating Polynomial. Consider the problem of finding the polynomial interpolation for the N + 1 = 4 data points

TWO-DI MENSIONAL INTERPOLATION 141

subject to the conditions

For this problem, we only have to type the following statements in the MATLAB command window.

3.7 TWO-DIMENSIONAL INTERPOLATION

In this section we deal with only the simplest way of two-dimensional interpolation—that is, a generalization of piecewise linear interpolation called

Figure 3.8 A two-dimensional interpolation using Zi = interp2() on the grid points [Xi,Yi] generated by the meshgrid() command.

142 INTERPOLATION AND CURVE FITTING

the bilinear interpolation. The bilinear interpolation for a point (x, y ) on the rectangular

sub-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 n

y 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 143

Example 3.5. Two-Dimensional Bilinear Interpolation. We consider interpolating the sample values of a function

for the 5 × 5 grid over the 21 × 21 grid on the domain

We make the MATLA B program “do_interp2.m”, which uses the routine

“intrp2()” to do this job, compares its function with that of the MATLAB

built-in routine “interp2()”, and computes a kind of relative error to estimate

how close the interpolated values are to the original values. The graphic results

of running this program are depicted in Fig. 3.9, which shows that we obtained

a 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 FITTING

When 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 FITTING

sample 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 the

cubic spline approach may be a roundabout toward the smoothness as explained

in 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 subinterval

needs 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, which

is 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 vertical

distance 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 Degree

If there is some theoretical basis on which we believe the relationship between

the two variables to be

we should set up the following system of equations from the collection of many experimental data:

Noting that this apparently corresponds to the overdetermined case mentioned

in Section 2.1.3, we resort to the least-squares (LS) solution (2.1.10)

which minimizes the objective function

CURVE FITTING 145

Sometimes 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 to

its accuracy/reliability. This policy can be implemented by the weighted leastsquares (WLS) solution

which minimizes the weighted objective function

If the weighting matrix is, then we can write the WLS solution (3.8.5) as

where

One may use the MATLAB built-in routine “lscov(A,y,V)” to obtain this

WLS solution.

3.8.2 Polynomial Curve Fit: A Polynomial Function of Higher Degree

If 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 θ as

The MATLAB routine “polyfits()” performs the WLS or LS scheme to

find 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 FITTING