Interpolation Sanzheng Qiao Department of Computing and Software - PowerPoint PPT Presentation

Intro Polynomial Piecewise Cubic Spline Software Summary Interpolation Sanzheng Qiao Department of Computing and Software McMaster University July, 2012

Intro Polynomial Piecewise Cubic Spline Software Summary Outline Introduction 1 Polynomial Interpolation 2 Piecewise Polynomial Interpolation 3 Natural Cubic Spline 4 Software Packages 5

Intro Polynomial Piecewise Cubic Spline Software Summary Introduction Problem setting: Given ( x 0 , y 0 ) , ( x 1 , y 1 ) , ..., ( x n , y n ) , x 0 < x 1 < · · · x n , for example, a set of measurements, construct a function f : f ( x i ) = y i , i = 0 , 1 , ..., n Desirable properties of f : smooth: analytic and | f ′′ ( x ) | not too large (the first and second derivatives are continuous). simple: polynomial of minimum degree, easy to evaluate.

Intro Polynomial Piecewise Cubic Spline Software Summary Example Measurements of the speed of sound in ocean water 5060 5040 5020 5000 4980 speed (ft/sec) 4960 4940 4920 4900 4880 4860 0 2000 4000 6000 8000 10000 12000 depth (ft)

Intro Polynomial Piecewise Cubic Spline Software Summary Example Interpolation 5060 5040 5020 5000 4980 speed (ft/sec) 4960 4940 4920 4900 4880 4860 0 2000 4000 6000 8000 10000 12000 depth (ft)

Intro Polynomial Piecewise Cubic Spline Software Summary Polynomial Interpolation Advantages: easy to evaluate and differentiate Weierstrass Approximation Theorem: If f is any continuous function on the finite closed interval [a,b], then for every ǫ > 0 there exists a polynomial p n ( x ) of degree n = n ( ǫ ) such that x ∈ [ a , b ] | f ( x ) − p n ( x ) | < ǫ. max Impractical (degree is often too high)

Intro Polynomial Piecewise Cubic Spline Software Summary A straightforward approach A polynomial of degree n is determined by its n + 1 coefficients. Given ( x 0 , y 0 ) , ..., ( x n , y n ) to be interpolated, we construct the linear system (Vandermonde matrix): x n       1 x 0 · · · a 0 y 0 0 x n 1 x 1 · · · a 1 y 1       1  =  . . .   .   .  . . . . .       · · · . . . . .      x n 1 x n · · · a n y n n solve for the coefficients of the polynomial p n ( y )( x ) = a 0 + a 1 x + · · · + a n x n

Intro Polynomial Piecewise Cubic Spline Software Summary Vandermonde matrix When x 0 , ..., x n are distinct, the Vandermonde matrix is nonsingular. Thus the system has a unique solution (coefficients of the interpolating polynomial). Example. Given two points ( 28 , 0 . 4695 ) and ( 30 , 0 . 5000 ) , we have the system � 1 � � a 0 � 0 . 4695 � � 28 = 1 30 a 1 0 . 5000 and the solution � a 0 � 30 � � 0 . 4695 � 0 . 04250 � = 1 � � − 28 = . a 1 − 1 1 0 . 5000 0 . 01525 2

Intro Polynomial Piecewise Cubic Spline Software Summary Vandermonde matrix problem: The coefficient (Vandermonde) matrix is often ill-conditioned question What is the condition number of the Vandermonde matrix constructed by x i = 2000 + i , i = 0 , 1 , ..., 7? Answer: 1 . 87 × 10 37

Intro Polynomial Piecewise Cubic Spline Software Summary Lagrange form (conceptually simple) Basis of polynomials: { l j ( x ) } ( j = 0 , 1 , ..., n ) of degree n such that � 1 , if i = j l j ( x i ) = otherwise 0 , construct x − x i � l j ( x ) = x j − x i i � = j Thus n � p n ( y )( x ) = l j ( x ) y j j = 0

Intro Polynomial Piecewise Cubic Spline Software Summary Example Given three points: ( 28 , 0 . 4695 ) , ( 30 , 0 . 5000 ) , ( 32 , 0 . 5299 ) , construct a second degree interpolating polynomial in the Lagrange form: ( x − 30 )( x − 32 ) p 2 ( x ) = ( 28 − 30 )( 28 − 32 ) 0 . 4695 ( x − 28 )( x − 32 ) + ( 30 − 28 )( 30 − 32 ) 0 . 5000 ( x − 28 )( x − 30 ) + ( 32 − 28 )( 32 − 30 ) 0 . 5299 p 2 ( 31 ) = 0 . 5150 ≈ sin ( 31 ◦ ) Hard to evaluate.

Intro Polynomial Piecewise Cubic Spline Software Summary Horner’s rule Evaluating a 0 x 3 + a 1 x 2 + a 2 x + a 3 Horner’s form: (( a 0 x + a 1 ) x + a 2 ) x + a 3 v = a(0); for (i = 1:n) v = v*x + a(i); end The optimal (most efficient and accurate) way of evaluating a 0 x n + ... + a n (not in the Lagrange form).

Intro Polynomial Piecewise Cubic Spline Software Summary Dangers of polynomial interpolation An example. Runge’s function (continuous derivatives of all order) 1 y ( x ) = on [ − 1 , 1 ] 1 + 25 x 2 equally spaces x 0 = − 1, x 1 , · · · , x n = 1 Equal Spacing (n = 13) 1 0.5 0 −0.5 −1 −1.5 −2 −2.5 −3 −3.5 −4 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 It is often best not to use global polynomial interpolation.

Intro Polynomial Piecewise Cubic Spline Software Summary Piecewise Polynomial Interpolation Given the partition α = x 1 < x 2 < · · · < x n = β, interpolate on each [ x i , x i + 1 ] with a low degree polynomial. Linear L i ( z ) = a i + b i ( z − x i ) , z ∈ [ x i , x i + 1 ] b i = y i + 1 − y i a i = y i , , 1 ≤ i ≤ n − 1 x i + 1 − x i

Intro Polynomial Piecewise Cubic Spline Software Summary Algorithm. Piecewise linear interpolation Given vectors x and y with interpolating points, this function returns the piecewise linear interpolation coefficients in the vectors a and b . function [a,b] = pwL(x,y) n = length(x); a = y(1:n-1); b = diff(y)./diff(x);

Intro Polynomial Piecewise Cubic Spline Software Summary Evaluation Given the piecewise linear interpolation L ( z ) represented by the coefficient vectors a , b , how do we evaluate this function at z ∈ [ α, β ] ? First, we locate [ x i , x i + 1 ] such that z ∈ [ x i , x i + 1 ] . Then, we evaluate L ( z ) using L i ( z ) . Search method: binary search, since x i are sorted. Observation: If [ x i , x i + 1 ] is associated with the current z , then it is likely that this subinterval will be the one for the next value.

Intro Polynomial Piecewise Cubic Spline Software Summary Algorithm. Locate Idea: Use the previous subinterval as a guess. If not, do binary search. Given the vector x of breakpoints and a scalar z between x 1 and x n , this function locates i so that x i ≤ z ≤ x i + 1 . The optional g is a guess. function i = Locate(x,z,g) if nargin==3 % try the guess if (x(g)<=z)&(z<=x(g+1)) i = g; return % quick return end end

Intro Polynomial Piecewise Cubic Spline Software Summary Algorithm. Locate (cont.) n = length(x); if z==x(n) i = n-1; % quick return else % binary search left = 1; right = n; while right > left+1 mid = floor((left + right)/2); if z < x(mid) right = mid; else left = mid; end end i = left; end

Intro Polynomial Piecewise Cubic Spline Software Summary Algorithm. pwLEval Given a piecewise linear interpolation coefficient vectors a and b from pwL and its breakpoints in x , this function returns the values of the interpolation evaluated at the points in z . function v = pwLEval(a,b,x,z) m = length(z); v = zeros(m,1); g = 1; for j=1:m i = Locate(x,z(j),g); v(j) = a(i) + b(i)*(z(j) - x(i)); g = i; end

Intro Polynomial Piecewise Cubic Spline Software Summary Example 1 1 y = ( x − 0 . 3 ) 2 + 0 . 01 + ( x − 0 . 9 ) 2 + 0 . 04 − 6 Interpolation of humps(x) with PWL, n = 10 100 80 60 40 20 0 0 0.5 1 1.5 2 2.5 3

Intro Polynomial Piecewise Cubic Spline Software Summary Problem setting Given ( x 1 , y 1 ) , ( x 2 , y 2 ) , ..., ( x n , y n ) , find s ( x ) : in each subinterval [ x i , x i + 1 ] , s ( x ) is cubic s ( x i ) = y i , i = 1 , ..., n s ′ ( x ) and s ′′ ( x ) are continuous at x 2 , x 3 , ..., x n − 1 s ′′ ( x 1 ) = s ′′ ( x n ) = 0 The second derivative of s ( x ) is zero at the end points means that s ( x ) is linear at the end points.

Intro Polynomial Piecewise Cubic Spline Software Summary A straightforward approach Suppose a i + b i x + c i x 2 + d i x 3 on [ x i , x i + 1 ] , i = 1 , ..., n − 1. 4 ( n − 1 ) unknowns to be determined. Interpolation: a i + b i x i + c i x 2 i + d i x 3 i = y i , i = 1 , ..., n − 1 a i + b i x i + 1 + c i x 2 i + 1 + d i x 3 i + 1 = y i + 1 , i = 1 , ..., n − 1 Continuous first derivative (consider [ x i − 1 , x i ] and [ x i , x i + 1 ] ): b i − 1 + 2 c i − 1 x i + 3 d i − 1 x 2 i = b i + 2 c i x i + 3 d i x 2 i , i = 2 , ..., n − 1 Continuous second derivative: 2 c i − 1 + 6 d i − 1 x i = 2 c i + 6 d i x i , i = 2 , ..., n − 1 Two end conditions: 2 c 1 + 6 d 1 x 1 = 0 and 2 c n − 1 + 6 d n − 1 x n = 0 Total of 4 ( n − 1 ) equations, a dense system.

Intro Polynomial Piecewise Cubic Spline Software Summary A clever approach: Constructing s ( x ) In the subinterval [ x i , x i + 1 ] , let h i = x i + 1 − x i and introduce new variables: ¯ w = ( x − x i ) / h i , w = 1 − w . Note: w ( x i ) = 0, w ( x i + 1 ) = 1 and ¯ w ( x i ) = 1, ¯ w ( x i + 1 ) = 0, (linear Lagrange polynomials). Thus wy i + 1 + ¯ wy i is the (linear) Lagrange interpolation on [ x i , x i + 1 ] . Construct i [( w 3 − w ) σ i + 1 + ( ¯ w 3 − ¯ wy i + h 2 s ( x ) = wy i + 1 + ¯ w ) σ i ] where σ i to be determined, so that the properties (the first and second derivatives are continuous) are satisfied.