翻译：插值与拟合

1、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 Newton polynomial. Given a function, the degree of the approximate polynomial, and the left/right boundary points

2、of the interval, the above MATLAB routine cheby()” uses this formula to make the Chebyshev polynomial approximation.The following example 川ustrates that this formula gives the same approximate polynomial function as could be obtained by applying the Newton polynomial with the Chebyshev nodes.Example

3、 3.1. Approximation by Chebyshev Polynomial. Consider the problem of finding the second-degree N = 2) polynomial to approximate the function2f (x) 1/(1 8x ). We make the following program do_cheby.m , which uses the MATLAB routine cheby()“ for this job and uses Lagrange/Newton polynomial with the Ch

4、ebyshev nodes to do the same job. Readers can run this program to check if the results are the same.do_cheby.nN ="?; fl = -2; b = 2;cBx1jy1 - chebyi f310 I 0°I 002 I iZ/0N d(x x ) d2(x x ) L +dz(x x )00(2) . 0(M N) . 0. rN,aRbi %Cheb/shev polynonial ftn%ft>r contparisori withpolynomial

5、ftnk - 0:N; m - cos( (2*N + 1 -个+ + 3 + :Chet)yshev nodesx - ( (b-a)-+ Eq. (3,3,1 b)y - f3i (x) ; n 二门stonpfC11 = lagrarp(x,y)chebyC = -0+3200-0.00001,00003.4 PADE APPROXIMATION BY RATIONAL FUNCTION/0、pm,n(x x )Qm(x x0)Dn(x x0)(3.4.1)Pade approximation tries to approximate a functiorf (x) around a p

6、ointxo by a rational functionq0 q1(x x0) q2(x x0)2+L +Qm(x x0)MHow do we find such a rational function? We write the Taylor series expansionof f (x) up to degreeM + N at x = xo as130 INTERPOLATION AND CURVE FITTINGf(x) Tm n(x x0)f(x0) f'(x°)(x x0)(2) / 0、f (x ) /0.2 .-(x x ) L2(M N) 0T (x )

7、 .0M N(x x ) (M N)!0020 M N出 a(x x ) 和(x x ) LaM n(x x) (3.4.2)Assuming x° = 0for simplicity, we get the coefficients of DN(x)andQM (x)such thatTM n(x)Qm(x)Dn(x)(a。ax L aM nxm N)(1 dx L dNxN) (q° qx L Qnxn) °1 d1x LdNxN,MN、“，.， N、,.N、(a0 aix LaM nx )(1 dx LdNx ) & qx LqNx ) (3.4.3

8、),d n and then substituted/ sby solving the following equations:a0q0ai加d1q1a2aa+ aod2q2LLLLLLaM+ aM 1d1+aM 2d2L+ aM N dNqM (3.4.4a ) (3.4.4b )aM 1+ aMd1+ aM 1d2L+ aM N 1dN0aM 2+ aM «+ aMd2LaM N 2dN0LLLLLLaM N+ aM N d1+ aM N d2L+ audN0Here, we must first solve Eq. (3.4.4b) for di,d2,into Eq. (3.

9、4.4a) to obtain q0 ,qi, , qMThe MATLAB routine padeap()" implements this scheme to find the coefficient vectors of the numerator/denominator polynomiaQM (x)/ Dn (x) of thePade approximation for a given functionf (x). Note the following things:? The derivativesf'(x0), f (x0),L , f(M N)(x0) u

10、p 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 x xo for x in pm,n(x) which has been obtained with the assumptionthat x0 = 0.PADE APPROXIMATION BY RATI

11、ONAL FUNCTION 131function numden =)%Input : f = function to be approximated around in xo, xfOutput: num = numerator coeffs of Pade approximation of degree M% den = denominator cceffs of Pade approximation of degree Na(i) = feval(ffxo);h = .01; tmp = 1;for i = 1: M + Ntmp = tmp*i*h; %ilhAidix = difap

12、x(i,-i l)*feval(fjXO+-i:i*h)1; %derivatlve(Ssction 5.3)a(i + 1) = dix/tmp； %Taylor series coefficient endfor m = 1 :Nn = 1 :N; A(nTn) = a(M + 1 + m - n);b(n) = -"M + 1 + n);endd = Ab' %Eq.(3.4.4b)for m = 1: M + 1mm = nin(in - 1,N);q(m) =d(1 :mm);*Eq, 3,4,4a)endnum = q(M + i:-1 :i)/d(N); den

13、 = d(N: 1:1), 1/d(N); %descending orderif nargout = 0 % plot the true ftn, Pade ftn and Taylor expansionif nargin <6, xo = xo - 1; xf = xo + 1； endx = x0+xf-x0/100*(0:100; yt = f9val(ftx);xl = x-xo; yp = polyval(nunfxl)./polyvalfden)xl);yT = polyval(a(M + N + 1:-1 ：1),x1);elf, plotfxt/k1,匹 yp J %

14、 KM 3) endExample 3.2. Pade Approximation for f(x) ex . Let ' s find the Pade approximationx0P3,2(x) Q3(x)/D2(x) for f(x) e around x =0. We make theMATLAB program do_pade.m' , which uses the routineoadeap。' " for this job and uses it again with no output argument to see the graphic

15、results asdepicted in Fig. 3.6.»do_pade %Pade approximationn = 0.33332.999611.999419.9908d = 1.0000 -7.999719.9988%d0_pade.il to get the Pade approximation for f(x) = eAxf1 = inline( 'exp(x)1 / x1);M = 3; N = 2; %the degrees of Numerator Q(x) and Denominator D(x) xo = 0; %the center of Tayl

16、or series expansionn,d = padeapff1jXOjMjN) %to get the coefficients of Q(x)/P(x)xO = -3,5; xf = 0.5; Ueft/right boundary of the interval padeaptfljXOjMjNfxOfXf) Mo see the graphic resultsTo confirm and support this result from the analytical point of view and to help the readers understand the inter

17、nal mechanism, we perform the hand-calculation whose coefficients are132 INTERPOLATION AND CURVE FITTING3.2.).procedure. First, we write the Taylor series expansion at = 0 up to degree 4 xx 一，、 xM + N = 5 for the given function f (x) e asTy(x)M N3xk1 x 1x3/23!k!(E3.2.1)/11 .a0 1,ai 1,a2,a3, L (E322)

18、26We put this into Eq. (3.4.4b) with M = 3,N = 2 and solve it for di' s to get2D2(x) 1 d1x d2xa4 a3dl a2d2 0,a3 a2dl a1d2 0(E3.2.3)1/6 1/2 d11/24 d12/51/24 1/6 d21/120 , d2a1/20Substituting this to Eq. (3.4.4a) yieldsqaa°d1q2a2 a1dlq。 a。1(E3.2.4)1 1 ( 2/5) 3/5a0d2 1/2 1 ( 2/5) 1 (1/20) 3/ 2

19、0 q3 a3 a2dl a1d2 1/6 (1/2) ( 2/5) 1 (1/20) 1/60INTERPOLATION BY CUBIC SPLINE 133With these coefficients, we write the Pade approximate function as(x) Qa1x)1 (3/5)x (3/20)x2 (1/60)x3自“D2(x)1( 2/5)x(1/20)x2(E3.2.5)(1/3)x33x2 12x20x28x 203.5 INTERPOLATION BY CUBIC SPLINEIf we use the Lagrange/Newton p

20、olynomial to interpolate a given set oN + 1 data points, the polynomial is usually of degreeN and so hasN - 1 local extrema (maxima/minima). Thus, it will show a wild swing/osc川ation (called' polynomialwiggle ' ), particularly near the ends of the whole interval as the number of data points

21、increases and so the degree of the polynomial gets higher, as 川ustrated 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

22、 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 adjacent subintervals can ' t be made to conform with eac

23、h other. In reals whlife, 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 rob

24、ot path/trajectory planning. That 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 po

25、ints.)For a given set of data points ( xk , yk), k 0 : N , the cubic splines(x)consists of N cubic polynomial sk(x) s assigned to each subinterval satisfyingthe following constraints (S0) <S4).( 50)s(x)sk(x)Sk,3(xxk)Sk,2(xxk)2Sk,1(xxk)Sk,0,forxxk,xk1,k 0: N( 51)sk(xk) Sk,0 yk， for k 0 : N( 52)sk

26、1(xk) sk (xk) Sk,0 yk,for k 1: N 1( 53)s'k 1(xk) s'k(xk) Sk,1,for k 1: N 1( 54)s''k1(xk) s''k(xk) 2Sk,2,for k 1: N 1These constraints (S1)(S4) amount to a set ofN + 1 + 3(N - 1) = 4N - 2 linear equations having 4N coefficients of the N cubic polynomials Sk,0,Sk,1,Sk,2,Sk,3,k

27、0: N 1134 INTERPOLATION AND CURVE FITTINGTable 3.4 Boundary Conditions for a Cubic Spline<i Firsl-order dcrivalivcs中工八> =写"(支内)=5川s|>euilicddi) Second-ordr埼(回) = 25oj,瑞(，内)=25M2den waives spec i lied<cnLkcurv:Uure jdjusicd)(in) Second-order derivjiivcs 5：；(工力=s(x) + 搭(£；(h。一 5，

28、(h2) cxlrupalatcd”力*1北 WG =小一1)-琦_2&n2)(3.5.1b)fiN-2&,3hkSk,2hkSk,1hkyk yk ito eliminate Sk,3 s and rewrite it as4(0 1,2Sk,2)Sk,2hkSk,1” yk dyk3hk(3.5.2a )(3.5.2b )hk(Sk 1,22$" 3Sk,13dykhk 1 ( 1,3hk 1 (S<,2 Q 1,2)3We substitute these equations into (S2) with< + 1 in place of k32sk

29、 (xk 1)Sk,3 (xk1xk)Sk,2(xk1xk )Sk,1(xk1xk) Sk,0yk1k,22& 1,2) 3Sk 1,1 3dyk 1We also substitute Eq. (3.5.1b) into (S3)s'k 1(xk)3Q1,351Sk1,2hk1Sk1,1Sk,1INTERPOLATION BY CUBIC SPLINE 135to writeSk,1 Sk 1,1 hk 1(Sk,2 Sk 1,2 ) 2hk 1Sk 1,2 hk 1(Sk,2 Sk 1,2) (3.5.3)In order to use this for eliminati

30、ng Sk,1 from Eq. (3.5.2), we subtract (3.5.2b)from (3.5.2a) to writehk(Sk 1,22Sk,2)hk1(Sk,22Sk1,2)3(Sk,12Sk1,1) 3(dykdyk1)and then substitute Eq. (3.5.3) into this to writehk(Sk 1,22Sk,2)hk1(Sk,22Sk1,2) 3hk1(Sk,2Sk 1,2)3(dykdyk1)(3.5.4)hk 1Sk 1,2 (hk 1hk )Sk,2 hkSk 1,2 3(dyk dyk 1 )for k = 1 : N - 1

31、Since these areN - 1 equations with respect toN + 1 unknowns Sk,2 ,k 0 : 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 ar

32、e given as (i) in Table 3.4, we write Eq. (3.5.2a) for k = 0 ash0(S1,2 2S0,2) 3S0,1 3dy0,2h0S0,2 h0S1,23dy0 3S0,1 (3.5.5a)We also write Eq. (3.5.2b) for k = N ashN 1(SN,2 2SN 1,2) 3SN 1,1 3dyN 13dyN 1(3.5.5b)and substitute (3.5.3)k( = N) into this to write h0 (S1,2 2S0,2 ) 3S0,13dy0hN 1(SN,22SN 1,2)

33、 3SN 1,13hN 1(SN,2 2SN 1,2)hN 1SN 1,2 2hN 1SN ,23(SN ,1 dyN 1 )Equations (3.5.5a) and (3.5.5b) are two additional equations that we need to solveEq. (3.5.4) and that s it. In the case where th-oerdsercodnedrivatives on the two boundary points are given as (ii) in Table 3.4, S0,2 and SN,2 are directl

34、y known from the boundary conditions asSo,2 s'0(Xo)/2, Sn,2 s''n(Xn)/2(3.5.6)136 INTERPOLATION AND CURVE FITTINGand, subsequently, we have jusN - 1 unknowns. In the case where the secondorder derivatives on the two boundary points are given as (iii) in Table 3.4$4 (上) (x 1 ) + y- ( $7(*

35、I ) 与，)fl IH Z _主'we can instantly convert these into two equations with respect tS0,2 and SN2 ash S()(2 (Aq + A| )5|j + h。= 0(3.5.7口)h拉 一 2sM2 gw-1 + An-2)Sn -1,2 + Hm-iSn -2,2 = 0(357b)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

36、form as2几ho0LLSd,23(dy。S0,1)ho2(ho h1)hLLS,23(dy1 dy°)0LLL0MM(3.5.8)LLhN 22(hN 2 hN 1)hN 1SN 1,23(dyN 1 dyN 2)LL0hN 12hN 1SN,23(Sn,1 dyN 1)After solving this system of equation for Sk2, k = 0 : N, we substitute them into (S1), (3.5.2), and (3.5.1) to get the other coefficients of the cubic spli

37、ne as竺X, Si/咤"dj% 一"(5从|,2+2Sg),""(3.5,9)33%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 expre

38、ssion, 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)

39、 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

40、 and fifth input arguments do not need to be fed.INTERPOLATION BY CUBIC SPLINE 137function yi, S = csplineThis function finds th* cubic splint* far th* input data point* lx.Y) 'Input: x = xO x1 i. xN, y = yQ y1 -，vN, xiint&riwlatior points %KC - 1/2 for 1st/2id derivatives cn t>c unci ary

41、 sp«cifi»d%KC = 3 for 2nd darWativ白 on bound口rp 白xtrap013tM零 dO = 5*(x0? - $G1: initial dorivati* dyN = S-xN)=钏1: final derivativeOutput; Stn.lc); n = 1：M K = 1,4 in descending orderif nargin y 6, dyFJ = 0; in% if rnroin J 6)d/0 = 0; «ndif nargin * 41KC = 0; #i)<iN = length(M) = 1;

42、为 constru cis 由 set of «t|u at lens u.r.t. (S( n3 2), n = 1 :N + 1)A = zeros(H + 1,N + 1i; b = zeros(N + 1,1);S = zti>s (H + 1,4);% Cubic spline co&fficlent natrixk - 1 ;Ni hk) = x(k + 1) * xk); dy(k1 - (y(k + 1 - y(k) )/h(k) i® Boundary conditionif KC <= 1 1 st derivatives 即&quo

43、t;IfiodA(191:2 = 2*11 CDb1 = 3*(dy(1J - dyO)； Eq. (3.5.5a)A|N + 1,N：N + 1) = h(Nj 2*h(N； b(N + 1) = 3«(dyN -明川；倏q. 8,5MKC = 2 2nd derivatives specifiedA(i ?1) = 2; bM) = dyO; A|N + 1 川+ 1) = 2; blN + 1) = dyN;号Eq,3,5.6 olta 与之nd dorivativa* oxtrapolatod111ftl = M0-h1) - h(2)A(N + 1|N-1；N 4 1 =

44、h(N - h(Nbl(N - 1> h(N - 1j endfor m = ! ；N SsEq. (3. 5,6)A(.再 1 :n + 1) = h(n - 1)- 1) + MW) h(m);h(fn) = 31dy (mJ - dy(m - 1);S(： J) = Ab；为 Cubic spline ccefficients+oni = 1: N和 *,4 =£) = $<%却幻 h(m); *Eq,修6 1S(flt2) dy(Ri) -h(n“史悻糖 + 1a>H2*S(ij3 i 8(1.1) y(MH色ndS = 8(1:N, 4:-1:1scend

45、ing orderpp = mkppfXjS); makfi ploc4wis0 poLynomiiiyi - ppval(ppjxi); 切守 1119s of pi0c9"is。polyncMiial ftn(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 splineinterpolation for the N + 1 = 4 dat

46、a points(0,0),(L )(2,4)，(3,5)(E3.3J)Go) = *。) = So j = 2. s"(川)=% (3) = ftjj = 2(E3.3.2)With the subinterval widths on the x-axis and the first divided differences as138 INTERPOLATION AND CURVE FITTINGwe write Eq. (3.5.8) as%q h=必=心=I芹Vtdy2 = = I (E3.3.3)hi2100s,23(d%So,i)3(E3.3.4) 3(E33.5)1410

47、s,23(dydy0)60141S,23(dy2dy)60 012s,23(S3,idy)3Then we solve this equation to get品上=一3, $.2 : 3,= - 3.and substitute this into Eq. (3.5.9) to obtainS0,oS0,iS,iS2,iS,3§,3S?,30,§,o1,S2,o4(E3.3.6)dyo h0 (S2,2 2So,2)2(E3.3.7a) 3dyI 9 s2,2 26,。2(E3.3.7b)dy2 |(S3,2 2sl,2) 2(E3.3.7c)Sl,2S0,23ho2(E

48、3.3.8a)6,2Sl,231TlS3,2S2,23h22(E3.3.8b)2(E3.3.8c)%do_csplines.niKC = 1; dyo - 2； dyn = 2； with specified 1st derivatives on boundaryx = O 1 2 3; y = O 1 4 5；Xi = x(1)+0:200*(x(erd)*x(1)/200; %intermediate points71,5 = csplinefx, KC,dyO,dyN); S %cubic spline interpolationelf pKo' iXipyi, 'k:&

49、#39;)yi - spline(x,dyO y dyNtxi); %for comparison with MATLAB built-in ftn hold on, pause, plot(x,yT * ro1 jXij/i,'r:1)yi = spline(x,yjXi); %for compar'ison with MATLAB built-in ftnpause, plot(x, bo1'b")KC = 3； yipSJ - csplinetx,ypxi,KC) ；%with the 2nd derivatives ext广即cdatecl pause

50、t plot(x,yJk。'JK')Figure 3.7 Cubic splines for Example 3.3.Finally, we can write the cubic spline equations collectively from (S0) as$0*) = $0,3(X Io)，+ So.2(X 1。广 + S。/。-v0)+ So.O213/+ 2 + 0l(-V)= S,3(K Xl)3 + S,2(工 >)"+ S1J(X - .X1) + S|、o=- 2(x 1)3 + 3(x I)2 + 2(* 1)+1$2。)= $33。-

51、*)3 +见2" -肛)2 + S2J(X M)+=2(x - 2)3 - 3a - 2)2 + 2(X - 1)+4We 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 cubicSa-,3, Sk,zy Sjm，S*小 k =0 : N 1spline function for xi (i.e., theinterpo

52、lation 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 §pline(x,y,xi) ” , which works with the boundary condition of type (i) for thesecond input argument given as dy0 y dyN, and with the boundary c

53、ondition of type (iii) for the same lengths of x and y.»do_csplin6s 兔cubic splineS= 2.0000-3.00002.00000-2.00003.00002.00001.00002.0000-3.00002.00004.00003.6 HERMITE INTERPOLATING POLYNOMIALIn some cases, we need to find the polynomial function that not only passes through the given points, but

54、 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 polynomialh(x) =+ Ho(361)matching just two points00,亢），（XQ1）and having the

55、 specified firstderivativesat the points. We can obtain the four coefficients : . """ bysolvingh(xo) = Ha反 + 8a丸 + */口 + Hi = y()= "%*： + 修工；+ "o ；(362)分'(大口)= 3H3飞 + ?-I- H =If (.t) = 373，+ 2/7*+ H = y jAs an alternative, we approximate the specified derivatives at the

56、data points by their differences, h(x(i +) - h(XQ) y2 To f /1U1) - Afjj -r) 甘一 先=, 月=££e8(3.63)and find the Lagrange/Newton polynomial matching the four points(m,比，(犯=十瓦 ys = >'o + 品£)，(刈=用 一 £, y3 = y -On, y) (3.64)The MATLAB routine hermit()" constructs Eq. (3.6.2) and solves it to get the Hermite interp