Download Polynomial Interpolation: Conditions for Unique Solution and Error Analysis and more Assignments Mathematics in PDF only on Docsity!
Numerical Computing S’10 MATH/CSCI 4800- Homework-4 Solution Outline
PROBLEMS
- (Pencil-and-paper)
(a) When, if ever, is it possible for two distinct polynomials to interpolate the same set of data points? (b) i. In interpolating a continuous function by a polynomial, what key features determine the approximation error? Support your answer by means of an appropriate mathematical expres- sion. ii. Under what circumstances does the error increase as the number of interpolation points is increased? Give an example. iii. Under what circumstances does the error decrease as the number of interpolation points is increased? Give an example. iv. Given a fixed number of data points, what effect does their placement have on the error?
(a) If there are n distinct data points, then there is a unique interpolating polynomial of degree at most n − 1. However, there is an infinity of higher-degree polynomials that interpolate the same data. For example, the straight line y = 1 interpolates the points (− 1 , 1) and (1, 1), but so does the parabola y = x^2 , the quartic y = x^4 , etc. (b) i. Interpolation of f (x) at n data points (xi, f (xi)), i = 1, 2 ,... n by the polynomial pn− 1 (x) in the interval [a, b] leads to the following expression for the error:
f (x) − pn− 1 (x) =
∏^ n
j=
(x − xj )
f (n)(cx) n! , cx ∈ [a, b].
Thus the error bound is determined by n, the number of points, the maximum size of the∏ n th derivative f (n)(x) in the interval [a, b], and by the maximum of the polynomial n j=1(x^ −^ xj^ ). ii. The error will increase with n if f (n)(x) increases rapidly with n. An example is the function f (x) = 1/(1 + 10x^2 ) discussed in class. iii. The error will decrease with n if f (n)(x) remains moderate. An example is the function f (x) = exp(x). iv. Uniformly-spaced data points lead to a larger error as compared to points that are more densely distributed near the ends of the interval than in the middle.
- (Pencil-and-paper) Consider the discrete data
x y 1 2 2 12 3 38 4 86
(a) Determine the polynomial interpolant to the data using the monomial basis. (b) Determine the Lagrange polynomial interpolant to the same data and show that it is equivalent to that obtained in part (a). (c) Determine the Newton polynomial interpolant to the same data, and show that it is equivalent to that obtained in part (a).
(d) For the Newton interpolant only, add the extra point (5, 100) to the data set and determine the revised interpolant.
(a) Let the interpolant be p(x) = c 1 + c 2 x + c 3 x^2 + c 4 x^3. Then the interpolation requirement leads to the linear system (^)
c 1 c 2 c 3 c 4
Matlab’s backslash command yields the solution
c 1 = 2, c 2 = − 3 , c 3 = 2, c 4 = 1,
yielding the interpolant p = 2 − 3 x + 2x^2 + x^3.
(b) The Lagrange form of the interpolant is
q(x) =
∑^4
i=
yiLi(x).
Here,
L 1 (x) =
(x − 2)(x − 3)(x − 4) (1 − 2)(1 − 3)(1 − 4)
(x^3 − 9 x^2 + 26x − 24),
L 2 (x) =
(x − 1)(x − 3)(x − 4) (2 − 1)(2 − 3)(2 − 4)
(x^3 − 8 x^2 + 19x − 12),
L 3 (x) =
(x − 1)(x − 2)(x − 4) (3 − 1)(3 − 2)(3 − 4)
(x^3 − 7 x^2 + 14x − 8),
L 4 (x) =
(x − 1)(x − 2)(x − 3) (4 − 1)(4 − 2)(4 − 3)
(x^3 − 6 x^2 + 11x − 6).
Then,
q(x) = 2 L 1 (x) + 12L 2 (x) + 38L 3 (x) + 86L 4 (x) = 2 − 3 x + 2x^2 + x^3.
Thus we obtain the same polynomial as p(x). (c) The Newton form of the interpolant is
r(x) = a 1 + a 2 (x − 1) + a 3 (x − 1)(x − 2) + a 4 (x − 1)(x − 2)(x − 3).
The interpolation conditions lead to the linear system
a 1 a 2 a 3 a 4
The solution to this lower triangular system is
a 1 = 2, a 2 = 10, a 3 = 8, a 4 = 1,
leading to
r(x) = 2 + 10(x − 1) + 8(x − 1)(x − 2) + (x − 1)(x − 2)(x − 3) = 2 − 3 x + 2x^2 + x^3.
S′′(x) =
4 − x on [0, 1] 2 b − (x − 1) on [1, 2] 2 − (x − 2) on [2, 3]. Now we check the various conditions. First, continuity of S and S′^ at the knots x = 1 and x = 2 requires the following four conditions.
4 + a + 2 −
a + 4 −
= c, 3 = 2 b, 2 b − 1 = 2.
All of these conditions are satisfied by a = − 29 / 6 , b = 3/ 2 , c = 7/ 6. Normally one would also check the interpolation conditions, but since we are not supplied with the data to be interpolated, no such conditions are given. Let us now check the remaining conditions. First, the not-a-knot condition. This requires continuity of S′′′(x) at the knots x = 1 and x = 2. The third derivative is given by
S′′′(x) =
− 1 on [0, 1] − 1 on [1, 2] − 1 on [2, 3].
Since S′′′^ = −1 everywhere, the not-a-knot condition is satisfied. Since there are only two knots, this means that the spline is in fact the same cubic over all the subintervals, given by
x + 2x^2 −
x^3.
- (Pencil-and-paper and Matlab) Consider the three equations
x^21 + x^22 + x^23 = 1 , 2 x^21 + x^22 − 4 x 3 = 0 , 3 x^21 − 4 x 2 + x^23 = 0.
We seek to find a solution to this set of equations by using Newton’s method for nonlinear systems. As we discussed in class, (also see section 2.7 of the text), Newton’s method applies to nonlinear systems of form fi(x 1 , x 2 , · · · xn) = 0, i = 1, 2 , · · · n. Given the n th iterate (^)
x( 1 n) x( 2 n) · · · x( nn)
the method computes the increment (^)
∆x( 1 n) ∆x( 2 n) · · · ∆x( nn)
by solving a related linear system J(n)∆x(n)^ = −f (n),
where J is the Jacobian matrix. The n + 1 th iterate is then given by x(n+1)^ = x(n)^ + ∆x(n).
(a) Compute by hand the Jacobian matrix of the given 3 × 3 system. (b) Implement Newton’s method numerically to solve the given system. Your code must print at each iteration the following information: the iteration number, the solution at that iteration, and the norm of the increment to the solution at that iteration. Terminate the solution when the norm is smaller than 10−^5. Use the initial guess x 1 = x 2 = x 3 = 0. 5. You may use as a guide the Newton script, available on line under the title M-files for root finding. This script is for finding the root of a single equation, and will need to be extended to a system.
(a) The system can be written as
f 1 (x 1 , x 2 , x 3 ) = x^21 + x^22 + x^23 − 1 = 0 , f 2 (x 1 , x 2 , x 3 ) = 2x^21 + x^22 − 4 x 3 = 0 , f 3 (x 1 , x 2 , x 3 ) = 3x^21 − 4 x 2 + x^23 = 0.
The Jacobian matrix is ∂f 1 ∂x 1
∂f 1 ∂x 2
∂f 1 ∂x 3 ∂f 2 ∂x 1
∂f 2 ∂x 2
∂f 2 ∂x 3 ∂f 3 ∂x 1
∂f 3 ∂x 2
∂f 3 ∂x 3
2 x 1 2 x 2 2 x 3 4 x 1 2 x 2 − 4 6 x 1 − 4 2 x 3
(b) The Newton iteration for the increment is found by solving the linear system
2 x 1 2 x 2 2 x 3 4 x 1 2 x 2 − 4 6 x 1 − 4 2 x 3
(n)
∆x 1 ∆x 2 ∆x 3
(n) = −
x^21 + x^22 + x^23 − 1 2 x^21 + x^22 − 4 x 3 3 x^21 − 4 x 2 + x^23
(n) .
The fresh iterate is given by
x 1 x 2 x 3
(n+1)
x 1 x 2 x 3
(n)
∆x 1 ∆x 2 ∆x 3
(n) .
The Matlab implementation and the results are given below. Note that three different Matlab scripts have been used, and convergence occurs in 7 iterations.
% Calling program for Newton3d tol=1e-5; m=20; x0=[0.5;0.5;0.5]; x=Newton3d(’f3d’,x0,tol,m) %%%%%% %function x= Newton3d(fun,x0,tol,m) % function x = Newton3d(fun,x0,tol,m) % % Solves a nonlinear system f(x)=0 by Newton’s method %
k = 3 x =
k = 4 x =
k = 5 x =
(Computer) One reason why polynomial interpolants are often used to approximate more complicated functions is because polynomials are easy to differentiate, integrate, and otherwise manipulate. This problem explores these ideas. Let f (x) be our example of a ‘complicated’ function on the interval I : [a, b].
(a) With the aid of the pseudocodes discussed in class, or otherwise, write one or more Matlab scripts to do the following. i. Find the coefficients of the polynomial
pn− 1 (x) = a 1 + a 2 x + a 3 x^2 + · · · anxn−^1
that interpolates f at n data points on I, supplied as the vector x = [x 1 x 2 x 3 · · · xn]. ii. Evaluate both pn− 1 and its derivative p′ n− 1 at a set of m values, supplied as the vector t = [t 1 t 2 t 3 · · · tm]. This evaluation should use the method of nested multiplication discussed in class and also in section 0.1 of the Sauer Text. iii. Evaluate exactly the integral (^) ∫ b
a
pn− 1 (t) dt.
These scripts should be of your own design, and should not employ Matlab’s built-in commands such as polyfit and polyval used in interpolation. (b) Let f (x) = tan x and let the interval I be [0, π/4]. Consider n = 5, 10 , 15 , and 20. i. Using Matlab’s subplot command, plot a set of four graphs. Each graph should correspond to a value of n specified above, and should show the plots of the derivatives p′ n− 1 (t) and f ′(t) = sec^2 (t) on the interval I using a grid of 100 equidistant points. Comment on the influence of n on the accuracy with which p′ n− 1 approximates f ′. ii. For each n compute the integral (^) ∫ π/ 4
0
pn− 1 (t) dt (1)
and compare it with the exact integral ∫ (^) π/ 4
0
tan x dx =
ln 2. (2)
Comment on the influence of n on the accuracy with which (1) approximates (2).
(a) See the Matlab codes main.m, interp.m and horner.m reproduced below.
% script main.m % close all clc format long e for k=1: n=5k; % supply interpolation data x=linspace(0,pi/4,n); y=tan(x); % supply plotting points t=linspace(0,pi/4,101); f=tan(t); fp=(sec(t)).^2; % compute coefficients a of polynomial interpolant p a=interpm(x,y); % compute coefficients ap of derivative p’ ap=ones(1,n-1); for i=1:n- ap(i)=ia(i+1); end % compute and plot f’ and p’ Fp=horner(ap,t); % plot figure subplot(2,2,k), plot(t,fp,’r’,t,Fp,’b’); legend(’Exact’,’Approximate’); % compute coefficients c of integral of p c=zeros(1,n+1); for i=1:n c(i+1)=a(i)/i; end % evaluate integral of p between 0 and pi/ Fint=horner(c,pi/4)-horner(c,0); % evaluate exact integral exact = 0.5*log(2); % print results fprintf(’n = %5i \n’,n); fprintf(’approx integral = %20.15e \n’,Fint); fprintf(’exact integral = %20.15e \n’,exact); fprintf(’error = %10.5e \n’,abs(Fint-exact)); end %%%%%%%%% %%%%%%%%% function a = interpm(x,y) % % Finds the coefficients a(n) of the polynomial that interpolates n given % points (x_i, y_i). Coefficients are found by using the monomial basis.
0.8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.
1
2 Exact Approximate
(^10) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.
2 Exact Approximate
0.8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.
1
2 Exact Approximate
(^10) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.
2 Exact Approximate
Figure 1: Plots of p′ n− 1 (approximate) and f ′^ (exact) for n = 5 (top left), n = 10 (top right), n = 15 (bottom left) and n = 20 (bottom right).
n = 5 approx integral = 3.466003539490644e- exact integral = 3.465735902799726e- error = 2.67637e-
n = 10 approx integral = 3.465736124394301e- exact integral = 3.465735902799726e- error = 2.21595e-
n = 15 approx integral = 3.465735902831348e- exact integral = 3.465735902799726e- error = 3.16219e-
n = 20 approx integral = 3.465735902799779e- exact integral = 3.465735902799726e- error = 5.21805e-
- (Computer) The following table lists the average high temperature for each month, January through
December, measured at the airport. Month Temp 1 54. 6 2 54. 4 3 67. 1 4 78. 3 5 85. 3 6 88. 7 7 96. 9 8 97. 6 9 84. 1 10 80. 1 11 68. 8 12 61. 1
(a) Find and plot the polynomial interpolant to the data. You may use the Matlab commands polyfit and polyval to do so. (b) Find and plot the not-a-knot cubic spline interpolant to the data. You may use the Matlab commands spline and ppval to do so. Compare your graph with that of part (a) and comment on which of the two is a more reasonable representation of the data. (c) Find and plot the clamped cubic spline interpolant to the data. This will require a slight modi- fication of the way in which the command spline is invoked. Compare your graph with that of part (b) and comment on the difference between the two plots. (d) Rearrange the data, so that the list goes from April through March rather than from January through December. Then repeat parts (b) and (c) and comment on your results, especially as regards the differences between the graphs of not-a-knot and clamped splines.
Be sure that each graph has the axes labelled, has a legend and a title, and displays the discrete data as circles.
See the Matlab script below. It requires the x and y values for the data as input. For parts (a - c) the data are supplied via the following Matlab commands.
x=1:12; y=[54.6 54.4 67.1 78.3 85.3 88.7 96.9 97.6 84.1 80.1 68.8 61.1];
For part (d) the rearranged data are
x=1:12; y=[78.3 85.3 88.7 96.9 97.6 84.1 80.1 68.8 61.1 54.6 54.4 67.1];
The script computes and plots the polynomial and spline interpolants.
%%%%%%%%%%%%% % Script File HW4P7.m % % INPUT: % data in the form of vectors x and y, each of length n % % OUTPUT: % polynomial and spline interpolants % % Check data dimension and consistency n = length(x); m = length(y); if m ~= n
(^300 2 4 6 8 10 )
40
50
60
70
80
90
100
Polynomial Interpolation
(^500 2 4 6 8 10 )
55
60
65
70
75
80
85
90
95
100 Not−a−Knot Spline Interpolation
(^500 2 4 6 8 10 )
55
60
65
70
75
80
85
90
95
100
Clamped Spline Interpolation
Figure 2: Interpolants for original data.
(^00 2 4 6 8 10 )
20
40
60
80
100
120
140
Polynomial Interpolation
(^500 2 4 6 8 10 )
55
60
65
70
75
80
85
90
95
100 Not−a−Knot Spline Interpolation
(^500 2 4 6 8 10 )
55
60
65
70
75
80
85
90
95
100
Clamped Spline Interpolation
Figure 3: Interpolants for rearranged data.
