Scientific Computing

1. Scientific Computing Interpolation – Polynomial Interpolation

2. Interpolation Definition: Given a set of (x,y) data points, Interpolation is the process of finding a function (usually a polynomial) that passes through these data points.

3. Polynomial Interpolation Definition: In polynomial interpolation, we find a polynomial (or set of polynomials) that passes through each data point. Simplest case: Linear Interpolation: Given two (different) data points (x0, y0) and (x1, y1) there is a unique line through these points.

4. Polynomial Interpolation Linear Interpolation: The equation of the line through (x0, y0) and (x1, y1) is Note: Consider the function Clearly, p(x0) = y0 and p(x1)=y1. Also, p(x) is linear (degree = 1). Thus, p(x) must also be the equation of the line through (x0, y0) and (x1, y1). The terms in p(x) are called Lagrange Polynomials

5. Polynomial Interpolation Definition: For a given set of n + 1 nodes xi the Lagrange polynomials are the n + 1 polynomials Li defined by Definition: We define a Lagrange Interpolating polynomial pn(x) by

6. Polynomial Interpolation Example: Consider the data (1,3), (1/2, -10), (3,2). There should be a quadratic polynomial passing through these three points. To find this polynomial, construct the three Lagrange terms:

7. Polynomial Interpolation Example: The Lagrange Interpolating polynomial is then Or,

8. Matlab Implementation function v = polyinterp2(x,y,u) % v = polyinterp(x,y,u) computes v = p(u) where p is the % polynomial of degree d = length(x)-1 with p(x(i)) = y(i). % Slightly modified from Moler textbook version n = length(x); v=0; % Stores p(u) for k = 1:n % Generate all terms of the Lagrange polynomial w = 1; % Stores one term for j = [1:k-1 k+1:n] % Computes the product of the fractions w = w*(u - x(j))/(x(k)-x(j)); end v = v + w*y(k); end

9. Matlab Example Example: We input the data (1,3), (1/2, -10), (3,2) by: x = [1 0.5 3]‘ y = [3 -10 2]‘ Then, we can test whether polyinterp2 works by checking its values on u=1, 0.5, 3 polyinterp2(x,y,1) ans = 3 polyinterp2(x,y,0.5) ans = -10 polyinterp2(x,y,3) ans = 2

10. Matlab Example Example: Let’s graph the x and y data and the values interpolated by polyinterp2. x = [1 0.5 3]'; y = [3 -10 2]'; u = -.25:.01:3.25; % A series of values along the x-axis [m n] = size(u); % polyvals stores the interpolating polynomial values polyvals = zeros(n); for i = 1:n % Compute each polynomial value polyvals(i) = polyinterp2(x,y,u(i)); end % Plot the x-y data as circles ('o') and the polynomial data as '-' plot(x,y,'o',u,polyvals,'-')

11. Matlab Example Example:

12. Matlab Example 2 Example: Here is an example with a larger data set x = 1:6; y = [16 18 21 17 15 12]; disp([x; y]) u = .75:.05:6.25; [m n] = size(u); polyvals = zeros(n); for i = 1:n polyvals(i) = polyinterp2(x,y,u(i)); end plot(x,y,'o',u,polyvals,'-')

13. Matlab Example 2 Example:

14. Polynomial Interpolation Theorem Theorem 5.2 (Interpolant Existence and Uniqueness). Let {xi} be n+1 distinct nodes. Then given y values at the nodes, {yi}, there is exactly one polynomial, p(x) of degree no greater than n such that p(xi) = yi Proof. We can construct the Lagrange Polynomial, so there is at least one polynomial p(x) with the desired properties. Suppose there were two polynomials that interpolated the points, call them p(x) and q(x), each of degree no greater than n. Let r(x) = p(x)-q(x). Then r(x) is a polynomial of degree n with (n+1) zeroes at the points (xi , yi ). The only polynomial with this property is the zero polynomial. Thus, 0=r(x) and so p(x) = q(x). QED

15. Newton Interpolating Method In Newton’s Interpolating method, we create a polynomial p(x) for data (xi, yi) such that p(xi) = yi (just as we did in the Lagrange Method). The difference is that we construct p(x) Iteratively. That is we construct a sequence of polynomials p0(x), p1(x), …, pn(x) such that pn(x) is the desired interpolating polynomial

16. Newton Interpolating Method Algorithm: Let p0(x) = y0 (p0(x) matches data at x0 ) Suppose we have calculated pk(x) Then, set pk+1 (x) = pk (x) + c(x – x0 )(x – x1 ) … (x – xk ) for some constant c. Note: pk+1 (xi ) = pk (xi ) for i = 0, 1, …, k (Why?). So, pk+1 (x) correctly interpolates all data up to xk. To make it match at xk+1 we set yk+1 = pk (xk+1) + c(xk+1 – x0 )(xk+1 – x1 ) … (xk+1 – xk ) and solve for c.

17. Newton Interpolating Method Example:Data = {(1,3), (1/2, -10), (3,2)} p0(x) = 3 p1(x) = p0(x) + c (x-1) = 3 + c (x-1) Set y1 = -10 = 3 + c (x1 -1) = 3 + c (1/2 – 1) Get c = 26. So, p1(x) = 3 + 26 (x-1) p2(x) = p1(x) + c (x-1)(x-1/2) = = 3+ 26 (x-1)+c (x-1)(x-1/2) Set y2 = 2 = 3+ 26 (3-1)+c (3-1)(3-1/2) Get c = -53/5 So, p(x) = 3+ 26 (x-1) – (53/5) (x-1)(x-1/2)

18. Newton Interpolating Method Example:p(x) = 3+ 26 (x-1) – (53/5) (x-1)(x-1/2)

19. Newton’s Method vs Lagrange Note: Both methods produce the same polynomial Newton’s Method is more flexible. It is easy to add new data points and get a new interpolating polynomial. The coefficients in Newton’s Method (determined by solving for the constants c) can be calculated very efficiently (explained next class period!).