Page 168 - Numerical Methods for Chemical Engineering

Polynomial interpolation 157

The coefficients of the polynomial therefore must satisfy the linear system
 2 N 
1 x ... x    
x 0
0 0 a 0 f 0
 1 x 2 ... N   f 1

x 1
1
1 x   a 1  
     
Xa =  1 x 2 x 2 2 ... x   a 2  =  f 2  = f (4.22)



N  

2

 . . .   .   . 
.
.
 . . . . . . .   .   . 
.
. 
1 x N x 2 N ... x N a N f N
N
As long as no two support points are colocated, det(X) = 0, and there is a unique polynomial
of degree N that solves the interpolation problem. However, we would like to obtain the
interpolating polynomial without having to solve a system of N linear equations by elim-
3
ination, requiring ∼ N FLOPs. Also, X easily may have quite a high condition number
and thus round-off errors may be significant, as we multiply, add, and subtract numbers of
vastly different magnitudes during elimination.
Lagrange interpolation
From the drawbacks of solving (4.22) by elimination, it would be nice to have a direct way
to construct p(x) from the function values. In fact, it is possible to express p(x) as the linear
combination
N

p(x) = f j L j (x) (4.23)
j=0
The {L j (x)} that satisfy p j (x) = f j (x) are the Lagrange polynomials
(x − x 0 ) ··· (x − x j−1 )(x − x j+1 ) ··· (x − x N )
L j (x) =
(x j − x 0 ) ··· (x j − x j−1 )(x j − x j+1 ) ··· (x j − x N )
(4.24)
N
x − x k
0
=
x j − x k
k=0
k = j
L j (x k ) = δ jk . In Figure 4.1, we plot the Lagrange polynomials for support points at
√
{0, 0.5, 1.0}. We use these polynomials to approximate f (x) = x on [0, 1] by (4.23). Note
that outside of [0, 1], p(x) and f (x) diverge rapidly. This demonstrates the drastic loss in
accuracy that may occur when extrapolating outside of the range of data using polynomials.
Newton interpolation
One disadvantage of Lagrange interpolation is that if we add another support point x N+1 ,we
have to recompute all of the Lagrange polynomials again from scratch. Newton interpolation
avoids this difficulty, and expresses the interpolating polynomial as
P N (x) = a 0 + a 1 (x − x 0 ) + a 2 (x − x 0 )(x − x 1 )
+··· + a N (x − x 0 )(x − x 1 )(x − x 2 ) ... (x − x N−1 ) (4.25)

We can evaluate P N (x) rapidly through the factorization
P N (x) = a 0 + (x − x 0 )[a 1 + (x − x 1 )[a 2 + (x − x 2 )[··· + (x − x N−1 )a N ]]] (4.26)