内插和外推方法.ppt

Chapter 3, Interpolation and Extrapolation,Interpolation & Extrapolation,(xi,yi),Find an analytic function f(x) that passes through given N points exactly.,低级多项式,高级多项式,低级多项式,高级多项式,Polynomial Interpolation,Use polynomial of degree N-1 to fit exactly with N data points (xi,yi), i =1, 2, , N. The coefficient ci is determined by a system of linear equations,Vandermonde Matrix Equation,But it is not advisable to solve this system numerically because of ill-conditioning.,Condition Number,cond(A) = |A| |A-1| For singular matrix, cond(A) = A linear system is ill-conditioned if cond(A) is very large.,Norms,Vector p-norm Matrix norm,sup: supremum,Commonly Used Norms,Vector norm Matrix norm,Where is the maximum eigenvalue of matrix ATA.,Lagranges Formula,It can be verified that the solution to the Vandermonde equation is given by the formula below: li(x) has the property li(xi)=1; li(xk)=0, ki.,Joseph-Louis Lagrange (1736-1813),Italian-French mathematician associated with many classic mathematics and physics Lagrange multipliers in minimization of a function, Lagranges interpolation formula, Lagranges theorem in group theory and number theory, and the Lagrangian (L=T-V) in mechanics and Lagrange equations.,Nevilles Algorithm,Evaluate Lagranges interpolation formula f(x) at x, given the data points (xi,yi). Interpolation tableau P(x),x1: y1 = P1 P12 x2: y2 = P2 P123 P23 P1234 x3: y3 = P3 P234 P34 x4: y4 = P4,Pi,i+1,i+2,i+n is a polynomial of degree n in x that passes through the points (xi,yi), (xi+1,yi+1), , (xi+n,yi+n) exactly.,Determine P12 from P1 & P2,Given the value P1 and P2 at x=x1 and x2, we find linear interpolation P12(x) = (x) P1+1-(x) P2 Since P12(x1) = P1 and P12(x2) = P2, we must have (x1) = 1, (x2) = 0 so (x) = 12= (x-x2)/(x1-x2),Determine P123 from P12 & P23,We write P123(x) = (x) P12(x) + 1-(x) P23(x) P123(x2) = P2 already for any choice of (x). We require that P123(x1)=P12(x1) =P1 and P123(x3)=P23(x3) =P3, thus (x1) = 1, (x3) = 0 Or,Recursion Relation for P,Given two m-point interpolated value P constructed from point i,i+1,i+2,i+m-1, and i+1,i+2,i+m, the next level m+1 point interpolation from i to i+m is a convex combination:,Use Small Difference C & D,P1,P2,P3,P4,P12,P23,P34,P123,P234,P1234,C1,1 = P12-P1,D1,1,C1,2,C1,3,D1,2,D1,3,C2,1=P123-P12,C2,2,D2,1,C3,1=P1234-P123,D3,1=P1234-P234,D2,2=P234-P34,Deriving the Relation among C & D,P = PA + (1-)PB,PA,PB,P0,C2=P-PA,C1=PB P0,D2=P-PB,D1=PA P0,Evaluate P(3) given 4-points (0,1), (1,2), (2,3),(4,0).,P1234= 11/4,C1,1 = 3,D1,1= 2,C1,2= 2,C1,3= -3/2,D1,2= 1,D1,3 = 3/2,C2,1= 0,C2,2= -5/3,D2,1= 0,C3,1= - 5/4,D3,1= 5/12,D2,2=5/6,x1=0, y1=1 =P1,x2=1, y2=2 =P2,x3=2, y3=3 =P3,x4=4, y4=0 =P4,C0,1=1,D0,1=1,C0,2=2,D0,2=2,C0,3=3,D0,3=3,D0,4=0,C0,4=0,polint( ) Program,polint(), continued,Subroutine polint(xa,ya,n,x,y,dy) Integer n, nmax Parameter (nmax=10) Integer i,m,ns Real den,dif,dift,ho,hp,w,c(nmax),d(nmax) Ns=1 Dif=abs(x-xa(1) Do 11 i=1,n Dift=abs(x-xa(i) If (dift.lt.dif) then Ns=i Idf=dift Endif,C(i)=ya(i) D(i)=ya(i) Enddo 11 Y=ya(ns) Ns=ns-1 Do 13 m=1,n-1 Do 12 i=1,n-m Ho=xa(i)-x Hp=xa(i+m)-x W=c(i+1)-d(i) Den=ho-hp If(den.eq.0) pause failure in polint Den=w/den,D(i)=hp*den C(i)=ho*den Enddo 12 If(2*ns.lt.n-m) then Dy=c(ns+1) Else Dy=d(ns) Ns=ns-1 Endif Y=y+dy Enddo 13 Return end,Piecewise Linear Interpolation,(x1,y1),(x2,y2),(x3,y3),(x4,y4),(x5,y5),P12(x),P23(x),P34(x),P45(x),x,y,Piecewise Polynomial Interpolation,(x1,y1),(x2,y2),(x3,y3),(x4,y4),(x5,y5),P1234(x),P1234(x),P2345(x),P2345(x),x,y,Discontinuous derivatives across segment,Cubic Spline,Given N points (xi,yi), i=1,2,N, for each interval between points i to i+1, fit to cubic polynomials such that Pi(xi)=yi and Pi(xi+1)=yi+1. Make 1st and 2nd derivatives continuous across intervals, i.e., Pi(n)(xi+1) = P(n)i+1(xi+1), n = 1 and 2. Fix boundary condition to P(x1 or N)=0, or P(x1 or N ) = const, to completely specify.,Curve Fitting by Interpolation,2. We now discuss Interpolation & Extrapolation The function passes through all (or at least most) points.,Curve Fitting: 1. We have discussed Least-Squares Regression where the function is “best fit“ to points but does not necessarily pass through the points.,Curve Fitting by Interpolation: C&C covers four approaches: Polynomials (C&C 18.1-18.5 “skim” only) n + 1 equations & n + 1 unknowns Lagrange Polynomials Newtons Divided Difference (NDD) Polynomials Splines (C&C 18.6 assigned reading) The first 3 approaches find the same polynomial. We will only cover superficially, and concentrate on Splines.,Curve Fitting by Interpolation: General Scheme Given: Set of points (xi,yi), not necessarily evenly spaced or in ascending order. Assume: x = independent variable; y = dependent variable. Find: y = f(x) at some value of x not in the set (xi). Method: Determine the function f(x) which passes through all (or most) points.,n+1 Equations and n+1 Unknowns (C&C 18.3) Given n+1 data points (xi,yi), find an nth order polynomial: y = a0 + a1x + a2x2 + + anxn that will pass through all the points.,Too much work! Equations are notoriously ill-conditioned for large n. Equations are not diagonally dominant. Method is rarely used.,Lagrange Interpolating Polynomials (C&C 18.2) Given n+1 data points (xi,yi), find the nth order polynomial: y = pn(x) = a0 + a1x + a2x2 + + anxn that passes through all of the points. The Lagrangian polynomials approach employs a set of nth order polynomials, Li(x), such that:,where Li(x) satisfies the condition:,Newtons Divided Difference (NDD) Polynomial (C&C 18.1),Gives the same polynomial as the Lagrange method but is computationally easier. General form for n+1 data points: pn(x) = b0 + b1(xx0) + b2(xx0)(xx1) + + bn(xx0)(xx1)(xx2)(xxn-1) with b0, b1, . , bn all unknown. Note that ith term is zero at xj for j i. Each term insures that the polynomial correctly interpolates at one new point. The algorithm is recursive and readily suited for spreadsheet or other programmed calculation.,Newtons Divided Differences (NDD) versus Lagrange Polynomials: 1. Both methods give the same results. 2. Comparison based on a count of the FLOPS: Evaluate coefficients: Interpolate for one x: Lagrange: (n+1) (n+1) n2 (n+1) (n+1) n2 NDD: (n+1) (n+1)/2 n2/2 n 3. Easy to add a node with NDD. Need to start over with Lagrange. 4. Both methods share a major problem: as the number of points increases, so does the order of polynomial. This may cause excessive “wiggles“ or “waves“ between points.,Splines (C&C 18.6) Issue: Need to overcome the “wiggle“ or “wave“ problem Idea: Use a piecewise polynomial approximation Simplest idea: Straight line on each segment. The problem is that g(x) and g“(x) are discontinuous,Splines (cont.),Most frequently used: Cubic Splines = Separate Cubic polynomial on each interval. This is the analytical/numerical analog of a flexible edge (spline) which is used by draftsmen. With this tool the first and second derivatives (curvature) are continuous, and the function appears “smooth“.,Cubic Splines,Objective: Define a 3rd-order polynomial for each interval: fi(x) = aix3 + bix2 + cix + di For n+1 data points, (x0,y0), (x1,y1), , (xn,yn), there are n intervals with 4 unknowns per interval (ai, bi, ci, and di). = Total 4n unknowns.,n=3, 3 segments, n+1 points,We need 4n equations to compute all 4 unknowns for each interval,Cubic Splines,How do we obtain the required 4n equations? functions must pass through fi(x) at knots (points): yi-1 = ai(xi-1)3 + bi(xi-1)2 + ci(xi-1) + di yi = ai(xi)3 + bi(xi)2 + ci(xi) + di 2 equations/interval = 2n 1st and 2nd derivatives must be equal at interior knots (xi,yi): 3ai(xi)2 + 2bixi + ci = 3ai+1(xi)2 + 2bi+1xi + ci+1 6aixi+ 2bi = 6ai+1xi+ 2bi+1 2 equations/interior knot = 2n-2 TOTAL: 2n + (2n - 2) = 4n - 2,Cubic Splines,We need an extra 2 conditions for Cubic Splines Natural Splines Setting the 2nd derivatives at exterior knots equal to zero allows the function to “relax“: 0 = 6a1x0 + 2b1 & 0 = 6anxn+1 + 2bn 1 equations/exterior knot = 2 TOTAL: 2n + (2n - 2) +2 = 4n,Cubic Splines,Alternatives for two extra conditions (instead of setting 2nd derivative = 0) 1) Specify 1st derivatives at exterior knots. f (x0) = 3a1x02 + b1x0 f (xn+1) = 3anxn+12 + bn xn+1 2) Add an extra point to the first and last intervals through which spline must pass: not-a-knot splines.,Cubic Splines Computation,If set up cleverly, the 4n x 4n system of equations can be reduced to solving an (n-1) x (n-1) tridiagonal system of equations. Define a new set of unknowns: Let si = f “(xi) be the second derivative of the cubic spline at interior point i, i = 1, ., n-1. First set up the n-1 equations to solve for curvatures, f “(x) at each of the interior knots (see C&C Box 18.3): (xi xi-1) s i-1 + 2 (x i+1 x i-1) s i + (x i+1 x i) s i+1,Natural Cubic Splines,Convenient tridiagonal equations for natural splines: These basic equations for the second derivatives can also be written in terms of the distances hi = (xi+1xi) between the points or knots and of yi = f(xi): For i = 0: s0 = 0 natural spline condition For i = 1: 2 (h0 +h1) s1 + h1s2 = RHS1 For i = 2 to n2: hi-1s i-1 + 2 (h i-1 + h i) s i + h is i+1 = RHSi For i = n1: h n-2 sn-2 + 2 (h n-2 + h n-1) s n-1 = RHSn-1 For i = n: sn = 0 natural spline condition,Natural Cubic Splines,Noting that this is a triangular banded system of equations of order n-1, we solve with a method which takes advantage of this, i.e., the Thomas algorithm given in C&C Section 11.1. Solve the triangular banded system, e.g., for n = 7 looks like: s1 s2 s3 s4 s5 s6 x x - - - - x x x - - - - x x x - - - - x x x - - - - x x x - - - - x x,Cubic Spline Interpolation,If we want to find specific interpolants, we do not need to determine the cubic functions in all of the intervals, rather just the interval in which the x lies. For the ith interval spanning xi-1,xi: fi(x) = ai(xxi-1)3 + bi(xxi-1)2 + ci(xxi-1) + di in which:,di = yi-1,Obtain a different cubic polynomial for each interval xi-1, xi. However, first need to solve for the values of all the si .,Natural Cubic Splines: Example,Set up equations to solve for the unknown curvatures at each interior knot for the following data:,(xi