Thompson - Computing for Scientists and Engineers (523188), страница 30
Текст из файла (страница 30)
Within each subintervals(x) is a cubic polynomial.2. Each of the derivatives with respect to x, namely s, s (l), s (2), is continuousover the whole range3. At each knot, x = xj, the spline fit goes exactly through the data, so thats(xj) = yj, j = 1,2,..., n.You can see in Figure 5.1 the cubics that are used to go through each pair of points,and how the slopes of successive splines match at the knots. Note, however, thateach cubic is used for only two points. In the figure the extensions beyond theknots are drawn to indicate how use of a given cubic segment outside its rangewould often give a poor representation of the curve.With the above three conditions on the spline fit, we are ready to derive an appropriate algorithm.
For simplicity, and as often encountered in practice, we assumeequal spacing of the independent variable values xj by an amount h.Deriving the spline equationsSince the fitting conditions relate to derivatives, we make Taylor (polynomial) expansions (Sections 3.1, 3.2) in the variable x about the point xj for the j th interval.For a cubic polynomial, derivatives past the third are zero, so we can immediatelycalculate the third derivative at xj in terms of the second derivatives at the knots:(5.1)Thus, the cubic spline can be written in terms of its first two derivatives, using alsothe requirement that it pass through the knot value yj,(5.2)Our goal is to solve for the derivatives in terms of the yj by using the spline continuity conditions, item (2) in the above list.
By differentiating (5.2) we obtain5.1HOW TO FIT CURVES USING SPLINES157(5.3)and, on differentiating this slope equation,(5.4)We can relate the first derivatives at the knots by using (5.3) with x = xj -1, to findthat(5.5)Using property (3) enables us to equate the datum yj +1 with the spline fit at xj +l, sothat in (5.2) we have(5.6)Stepping back one point gives similarly(5.7)We can now obtain a relation between adjacent spline second derivatives by subtracting the second equation from the first, substituting (5.5), and rearranging, to obtainfor j = 2,3, . . ..n - 1(5.8)where the second differences of the data are(5.9)at the jth point.
This formula is arranged so that roundoff errors are minimizedwithout increasing the number of arithmetic operations, as discussed in Section 4.4.The equations in (5.8) form a set of n - 2 linear equations which can be solvedonce the derivatives at the endpoints have been chosen by the user. As we discussin more detail in Section 5.2, the most natural spline is one in which these derivatives are zero, that is, the fit sticks out straight at each end. Whatever our choice ofendpoints, the same method of solving the equations (5.8) can be used. We maythink of the relations (5.8) as defining the elements of matrices, with those on theleft forming a banded matrix, with elements only on the diagonal, j, and one placeabove and below the diagonal, j ± 1.
This property has arisen because we are using a spline of cubic order, whereas, one of higher order would have more off-diagonal elements. Because, also for simplicity, we are assuming equally spaced data, itis straightforward to solve the equations directly, rather than using the full power ofmethods useful for matrix inversion.158FITTING CURVES THROUGH DATAWe solve equations (5.8) iteratively by first eliminatingthese equations. Then we setbetween pairs of(5.10)(5.11)The spline equations (5.8) then become(5.12)From this, starting at j = n - 1 and working downward to smaller j, we find thespline second derivatives at the knots(5.13)The reader who is experienced in matrix operations may notice that we have just performed Gaussian elimination to find the inverse of a tridiagonal matrix. This connection is used in more general treatments of spline fitting, such as those in themonographs by De Boor and by Schumaker.Exercise 5.3 Verify all of the above steps in the spline-fitting derivation.(Spare the spline and spoil the child.) nWe have now completed the analysis required for deriving cubic-spline fits, andit is time to cast the formulas into an algorithm suitable for programming.The spline algorithmThe algorithm for cubic spline fitting of n data values y1,y2,...,yn can now be summarized as follows:(1) Compute the coefficients aj according to the iteration scheme(5.14)These spline coefficients depend only on the number of data points, n, and not onthe data values.
They could be computed and stored for reuse, provided that n didnot change.(2) With j increasing compute the first differences among the data values(5.15)and thence the second-difference quantities(5.16)5.2BOUNDARY CONDITIONS FOR SPLINE FITTING159(5.17)(5.18)and the coefficients(5.19)Note that the second derivatives at the boundaries must be supplied. They are zerofor the natural cubic spline.(3) With j decreasing, compute the second derivatives at the spline knots(5.20)(5.21)(4) With j increasing, compute the first derivatives at the knots(5.22)(5.23)and the third derivatives at the same points(5.24)(5.25)Steps (1) through (4) summarize the algorithm for generating the spline-fitting coefficients.
In the next two sections we discuss the boundary conditions for cubicspline fits, then the coding and testing of a program in the C language.5.2 BOUNDARY CONDITIONS FOR SPLINE FITTINGBecause a spline is a piecewise-continuous function, but with different coefficientsfor the cubic between each pair of knots, the derivatives at the boundaries x1 and xn(which are usually called the “endpoints”) have significant effects on the behavior ofthe fit even at the interior points. Therefore, we have to discuss the boundary conditions for spline fitting, for example (for a cubic spline) the second derivatives at eachendpoint.160FITTING CURVES THROUGH DATANatural splinesIf you use a flexible strip of wood and force it to pass through several points on acurve, it bends between the points but it sticks out straight at both ends, so that itssecond derivatives are zero.
This condition, when applied at both ends, is called thenatural-spline endpoint condition. Note that this naturalness comes from the draftinganalogy, rather than being related to the properties of the underlying function that weare trying to approximate.An interesting property of cubic splines with the natural endpoint conditions imposed is that they have a minimal value of the curvature, C, where C is defined by(5.26)For the natural-spline fit, sn (x), let us call the curvature CN, so that(5.27)Our claim is then that C CN for any cubic-spline fit. This condition is also called“minimum tension” because it corresponds to a the behavior of a flexible strip (themechanical drafting spline) when constrained to touch the knots but otherwise to beunconstrained.Exercise 5.4Derive the following result for the second derivatives (curvature) of spline fits:(5.28)To do this, compute the second integral on the right-hand side by expanding thesquare, then split the integral into a sum of integrals over segments of length h,integrate by parts, and resum.
Complete the proof by using the natural-splineendpoint conditions. nThus, since the integral of the squared difference of spline curvatures in (5.28) is necessarily nonnegative, the least curvature is obtained when the spline fit, s, coincideswith the natural-spline fit, sN, at all x in the range of integration.The effects of the endpoint conditions on the overall closeness of the spline fitcan be seen for the pathological case of our working function (4.1) with its closelyspaced roots and consequent very wiggly behavior over a small range of x values.With natural endpoint conditions the spline fit is as shown in Figure 5.1. By comparison, Figure 5.2 shows the working function with the same knots fitted by thecubic spline with the exact endpoint second derivatives, obtained by using in (4.3)the coefficients w2,k from Table 4.1.5.3PROGRAM FOR SPLINE FITTING161FIGURE 5.2 The working function (4.1) and its fit by a cubic spline through five knots(indicated by the circles) using exact endpoint conditions.
The solid curve is the working functionand the dashed curves are the cubics used between each pair of knots, with the cubic shown beyondthe region in which it is used, thus showing the continuity of the spline curvature.We make further comparisons of natural-spline and exact endpoint conditions,for both the badly behaved working function and the well-behaved cosine function,in Sections 5.4 and 5.5.5.3PROJECT 5: PROGRAM FOR SPLINE FITTINGWe now have the spline algorithm summarized at the end of Section 5.1, togetherwith some understanding of endpoint conditions (Section 5.2). Therefore, we areready to program, test, and use the program for spline fitting presented as Program 5.1.PROGRAM 5.1 Cubic-spline interpolation, differentiation, and integration with equally spacedknots and user-controlled endpoint conditions.#include <stdio.h>#include <math.h>#define MAX 101main(){/* Cubic Splines */doubleylist[MAX],sdl[MAX],sd2[MAX],sd3[MAX];doublexmin,xmax,h,x,sx,sdlx,sd2x,sd3x;162FITTING CURVES THROUGH DATAdoubleywl,yw2,yw3,yerr,ylerr,y2err,y3err;double xL,xU,Intx,yInterr;int n,j;char Interp,Integ;voidSplineFit(),SplineInterp();double yw(),SplineInt(),ywInt();printf("Cubic Splines");n = 3;while (n>O)printf("\nInput n <=%i (n=O to end): ",MAX-1);scanf ("%i", &n) ;if(n<l)printf("\nEnd Cubic Splines") ; exit(O);}if(n>MAX-1)printf("\n !! More than MAX = %i data points\n", MAX-l);else{if (n < 4){printf("\n !! Fewer than 4 data points\n");elseprintf("\nInput %i data: ",n);for ( j = 1; j <= n; j++ )scanf ("%lf",&ylist[j] ) ;printf("\nInput second derivatives at j=l & j=%i: ",n);scanf("%lf%lf",&sd2[1], &sd2[n] ) ;printf("\nInput minimum & maximum x: ");scanf ("%lf%lf",&xmin,&xrmx) ;h = (xmax-xmin)/(n-1); /* uniform step size */SplineFit(n,h,ylist,sdl,sd2,sd3);x = (xmin+max)/2;while ( x >= xmin && x <= xmax )printf("\nInterpolate? (y or n) ");if ( Inter-p == 'y' )scanf("%s",&Interp)5.3PROGRAM FOR SPLINE FITTING163printf("Input x (x<%g or x>%g to end): ",xmin,xmax);scanf("%lf",&x) ;if ( x >= xmin && x <= xmax ){SplineInterp(x,xmin,h,ylist,sd1,sd2,sd3,&sx, &sd1x, &sd2x, &sd3x) ;printf("Interpolated s=%g, s'=%g, s' '=%g, s' ' '=%g\n",sx, sd1x, sd2x, sd3x) ;ylerr = ywl-sd1x;yerr = yw(x,&yw1,&yw2,&yw3)-sx;y2err = yw2-sd2x;y3err = yw3-sd3x;printf("Error in y=%g, y'=%g, y''=%g, in y'''=%g\n",yerr,y1err,y2err,y3err);}}printf("\nIntegrate? (y or n) :"); scanf("%s",&Integ) ;if ( Integ == 'y' ){printf("Input xL (>=%g) xU (<=%g) : ",xmin,xmax);scanf("%lf%lf",&xL,&xU) ;if ( xL >= xmin && xU <= xmax ){if ( xU >= xL ){Intx = SplineInt(xL,xU,h,xmin,ylist,sd1,sd2,sd3);yInterr = (ywInt(xU)-ywInt(xL))-Intx;printf("Integra1 = %g & error = %g\n",Intx,yInterr);}elseprintf("\n !! Negative range (%g to %g)\n",xL,xU);}else printf("\n !! A limit is outside spline range\n");if ( (Interp != 'y') && (Integ != 'Y') ){x = 2*(fabs(xmin)+fabs(xmax)); /* forces x loop exit */}}/* end while x interpolation loop */}}} /* end while n loop */}voidSplineFit(n,h,ylist,sd1,sd2,sd3)/* Cubic spline fit; uniform step size (h) in x */doubleylist[],sdl[],sd2[],sd3[];double h;164FITTING CURVES THROUGH DATAint n;{doubleas[MAX],es[MAX],ds[MAX],bs[MAX];double hfact,hd2;int j;hfact = 6/(h*h); hd2 = h/2;/* Spline coefficients */as[2] = 4;for ( j = 3; j < n; j++ ){as[j] = 4-l/as[j-1];for ( j = 2; j <= n; j++ ){ /* First differences */es[j] = ylist[j]-ylist[j-1];for ( j = 2; j < n ; j++ ){ /* Second differences */ds[j] = hfact*(es[j+l]-es[j]);}ds[n-1] = ds[n-l]-sd2[n];/* Make b coefficients */bs[2] = ds[2]-sd2[1];for ( j = 3; j < n; j++ ){bs[j] = ds[j]-bs[j-l]/as[j-1];}/* Second derivatives of spline */sd2[n-1] = bs[n-l]/as[n-1];for ( j = n-2; j >= 2; j-- ) /* downward recurrence */{sa[j] = (bs[j]-sd2[j+l])/as[j];/* First derivatives of spline */sd1[l] = es[2]/h-(sd2[1]/3+sd2[2]/6)*h;for ( j = 2; j <= n; j++ ){sd1[j] = sd1[j-l]+(sd2[j-l]+sd2[j])*hd2;}/* Third derivatives of spline */for ( j = 1; j < n; j++ ){sd3[j] = (sd2[j+l]-sd2[j])/h;}sd3[n] = sd3[n-1];5.3PROGRAM FOR SPLINE FITTINGvoid SplineInterp(x,xmin,h,ylist,sdl1sd2,sd3,sx, sd1x, sd2x, sd3x)/* Cubic spline interpolation */doubleylist[],sdl[],sd2[],sd3[];double x,xmin,h;double *sx,*sdlx,*sd2x,*sd3x; /* interpolated values */{double e; /* distance above i-th point */int i; /* interpolation index */i = (x-xmin)/h+l.l; e = x-xmin-(i-l)*h;/* Interpolated value, sx, */*sx = ylist[i]+(sd1[i]+(sd2[i]/2+sd3[i]*e/6)*e)*e;*sd1x = sd1[i]+(sd2[i]+sd3[i]*e/2)*e; /* first derivative */*sd2x = sd2[i]+sd3[i]*e; /* second derivative */*sd3x = sd3[i]; /* third derivative */}double SplineInt(xL,xU,h,xmin,ylist,sd1,sd2,sd3)/* Cubic spline integral */doubleylist[l,sd1[],sd2[],sd3[];double xL, xU, h, xmin;double eL, eU,ysum, suml, sum2, sum3 ;doubleyU,s1U,s2U,s3U,yL,s1L,s2L,s3L,splint;int jL,jU,j;double Horner4Poly();jL = (xL-xmin)/h+l.Ol; eL = xL-xmin-(jL-l)*h;jU = (xU-xmin)/h+l.Ol;eU = xU-xmin-(jU-l)*h;ysum = sum1 = sum2 = sum3 = 0; /* sums over whole strips */for ( j = jL; j < jU; j++ )ysum = ysum + ylist[j];sum1 = sum1 + sd1[j];sum2 = sum2 + sd2[j];sum3 = sum3 + sd3[j];}splint = Horner4Poly(ysum,suml,sum2,sum3,h);/* Contributions from partial strips at ends */yU = ylist[jU]; s1U = sd1[jU]; s2U = sd2[jU]; s3U = sd3[jU];yL = ylist[jL]; s1L = sd1[jL]; s2L = sd2[jL]; s3L = sd3[jL];splint = splint + Horner4Poly(yU,s1U,s2U,s3U,eU)- Horner4Poly(yL,s1L,s2L,s3L,eL);165166FITTING CURVES THROUGH DATAreturn splint;}double Horner4Poly(y,s1,s2,s3,e)/* Homer 4th-order polynomial algorithm */double y,sl,s2,s3,e;{double poly;poly = (y+(sl/2+(s2/6+s3*e/24)*e)*e)*e;return poly;double yw(x,yw1,yw2,yw3)/* Working function and its first 3 derivatives */double x,*ywl,*yw2,*yw3;{double y;y = 120*(x+0.5)*(x+0.25) *x*(x-(1.0)/(3.0))*(x-0.2)*(x-l);*yw1 = -1+(6+(69+(-204+(-470+720*x)*x)*x)*x)*x;*yw2 = 6+(138+(-612+(-1880+3600*x)*x)*x)*x;*yw3 = 138+(-1224+(-5640+14400*x)*x)*x;return y;double ywInt(x)/* Working function integral */double x;{double yint;yint = (-0.5+(1+(5.75+(-10.2+((-47.0/3.0)+(120.0/7.0)*x)*x)*x)*x)*x)*x*x;return yint;We now describe the structure of the major parts of the spline program, whichare the main program and the fitting function SplineFit.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
