Description
1 Introduction
This lab is concerned with interpolating data with polynomials and with trigonometric functions. There is
a difference between interpolating and approximating. Interpolating means that the function passes through
all the given points with no regard to the points between them. Approximating means that the function is
close to all the given points and and has some additional mathematical properties often making it close to
nearby points as well. It might seem that interpolating is always better than approximating, but there are
many situations in which approximation is the better choice.
Polynomials are commonly used to interpolate data. There are several methods for defining, constructing,
and evaluating the polynomial of a given degree that passes through data values. One common way is to
solve for the coefficients of the polynomial by setting up a linear system called Vandermonde’s equation.
A second way is to represent the polynomial as the sum of simple Lagrange polynomials. Other ways are
available and are discussed in the text and the lectures. In addition, we will discuss interpolation using
trigonometric polynomials.
We will examine these two ways of finding the interpolating polynomial and a way of finding the interpolating trigonometric polynomial. We will see examples showing that interpolation does not necessarily mean
good approximation and that one way that a polynomial interpolant can fail to approximate is because of a
bad case of “the wiggles.” Trigonometric polynomial interpolation does better, but can also break down.
1.1 Notation
This lab is focussed on finding functions that pass through certain specified points. In customary notation,
you are given a set of n points (xk, yk) and you want to find a function f(x) that interpolates this data. That
is so say, f(xk) = yk for all k. The pairs of (x, y) values to be interpolated will be called the data points.
When possible, we will use notation that distinguishes the data values with a special name. Thus the
names xdata and ydata will be preferred to denote the vectors containing xk and yk. We plan to construct
an interpolating function f(x), and plot it. To do so, we will need to create a second set of vectors, which
we might think of as sample points, for which we will prefer names like xval and yval. Thus, in a typical
computation, we will be given the information xdata, ydata, we will construct a function y = f(x), we
will choose a set of sample x values xval at which to sample the function, and do evaluations like yval =
f(xval).
Once we have the function f(x), we can also verify the interpolation property by evaluating it at each x
data value and comparing to the specified y data value: f(xdata)
?= ydata
2 Computing Tips
At some point in this lab, you will need to determine if an integer m is not equal to an integer n. Such a
condition might be expressed as:
i f ( m != n ) :
do something
1
Often, we will write a function file f.py which takes as input an argument x and evaluates f(x). If x is a
vector, we may need to know how big it is, so we can create a result vector just as big. One way to do this
is with the command
n = len ( x )
This works fine if x is actually a vector, that is, an np.array(). However, this will fail if x is a scalar, that
is, if we simply called the function with a command like
y = f ( 2 . 0 )
which will cause the error message
TypeError: object of type ’float’ has no len()
We can avoid this problem by modifying the text of the function to read:
x = np . a t l e a s t 1 d ( x )
n = len ( x )
The np.atleast 1d() command “promotes” any scalar value to a vector of length 1. The function now will
work the same way, whether we give it a scalar or a vector input. Don’t worry about this peculiar problem.
In the cases where it might affect you, we have indicated how to use np.atleast 1d() to avoid it!
In some cases, you will also need to solve a linear system, which we can represent as: A ∗ x = b, where A
is an n by n matrix, and b and x are column vectors. This can be done via the solve() command:
x = np . l i n a l g . s o l v e ( A, b )
A polynomial is represented by the vector of its coefficients. There are several functions for dealing with
polynomials, available through the import statement. If you need several of these functions, you can just
import them all, by specifying * rather than a list of specific names.
from numpy . p ol yn omi al . p ol yn omi al import ∗
c=polyfromroots(r) Finds the coefficients c of the polynomial whose roots are given by r.
r=polyroots(c) Finds the roots r of the polynomial whose coefficients are given by c.
c=polyfit( xdata, ydata, d) Finds the coefficients c of the polynomial of degree d passing through, or
approaching as closely as possible, the points xdata,ydata.
yval=polyval(xval,c) Evaluates a polynomial with coefficients c at the values xval.
Note that polyfit() can interpolate or approximate data. What happens depends on the relationship
between ndata, the number of data pairs, and the value of d. For our interpolation work, we will always
choose the degree to be d = ndata – 1. In other words, a constant (degree 0) polynomial can interpolate 1
value, a linear (degree 1) polynomial can interpolate 2 values, and so on.
In Python, the coefficients ck of a polynomial of degree d are defined as a vector of length d + 1, indexed
starting at 0, with ck the coefficient of x
k
. Therefore, a polynomial of degree d can be evaluated by:
p(x) = X
d
k=0
ck x
k
(1)
The ordering of coefficients is the reverse of what Matlab uses.
2
3 Exercise 1
Let’s just practice using Python’s functions that support working with polynomials. Create a Python script
exercise1.py and gradually carry out each of the following steps:
1. Begin your script with these useful import statements:
import m a t pl o tli b . p y pl o t a s p l t
import numpy . p ol yn omi al . p ol yn omi al a s p ol yn omi al
import numpy a s np
2. Set roots r, and use polyfromroots() to get coefficients c, use polyval() to evaluate the polynomial,
and then plot the polynomial and the roots! Once this works, you can add the next step to your script.
r = np . a r r a y ( [ 0 . 0 , 2 . 0 , 2 . 0 , 3 . 5 , 5. 0 ] )
c = p ol yn omi al . p ol y f r om r o o t s ( r )
x v al = np . l i n s p a c e ( −0.1 , 5 . 1 , 100 )
y v al = p ol yn omi al . p ol y v al ( xval , c )
p l t . pl o t ( xval , y v al )
p l t . pl o t ( xval , 0. 0 ∗ xval , ’ k−−’ ) # Se t 0 y v a l u e s f o r t h e x a x i s .
p l t . pl o t ( r , 0. 0 ∗ r , ’ r . ’ , m a r k e r si z e = 15 ) # Se t 0 y v a l u e s f o r t h e r o o t s .
p l t . g ri d ( True )
p l t . s a v e f i g ( ’ e x e r c i s e 1 p l o t 1 . jp g ’ )
p l t . show ( )
3. Set coeficients c, then use polyroots() to get roots. Again, evaluate polynomial, and display plot of
polynomial, x axis, and roots. Once this works, you can add the next step to your script.
c = np . a r r a y ( [ 1 0. 0 , −4.0 , −7.0 , 1. 0 ] )
r = ?
x v al = np . l i n s p a c e ( ? )
y v al = p ol yn omi al . p ol y v al ( ? )
p l t . pl o t ( ? )
p l t . pl o t ( ? )
p l t . pl o t ( ? )
p l t . g ri d ( ? )
p l t . s a v e f i g ( ? )
p l t . show ( )
4. Instead of using roots, supply data values, and ask polyfit() to find the coefficients of the polynomial
that goes through these values. Plot the polynomial, the x axis, and the pairs of (xdata,ydata)
points.
xdata = np . a r r a y ( [ −1.0 , 0 . 0 , 1 . 0 , 2. 0 ] )
ydata = np . a r r a y ( [ 8 . 0 , 4 . 0 , 4 . 0 , 1 4. 0 ] )
d e g r e e = 3
. . .NOW WHAT? . . .
Now you should be comfortable with the idea that a polynomial can be described by
• a set of coefficients c, or
• a set of roots r, or
• a set of data pairs (xdata,ydata);
and you should realize that Python gives you tools to convert from one description to another, and to use
plots to verify that your results are reasonable.
3
In this and some later labs, you will be writing codes analogous to polyfit() and polyval(), using
several different methods. Functions with the prefix coef will generate a vector of coefficients, as does
polyfit(), and functions with the prefix eval will evaluate a polynomial (or other function) at values of
xval, as does polyval().
To test our codes, we will start with a known polynomial and use it to generate “data” values. Now, using
only the data, our codes are asked to generate function values at intermediate points. These interpolated
values can then be compared to the exact values associated with the original polynomial.
4 Vandermonde’s Equation
Here’s one way to see how to organize the computation of the polynomial that passes through a set of data.
Suppose we wanted to determine the linear polynomial p(x) = c0 + c1x that passes through the n = 2
data points (x0, y0) and (x1, y1). We simply have to solve a set of linear equations for c0 and c1 constructed
by plugging in the two data points into the general linear polynomial. These equations are
c0 + c1x0 = y0
c0 + c1x1 = y1
or, written as a linear system,
1 x0
1 x1
c0
c1
=
y0
y1
which has the solution
c0 = (y1 − y0)/(x1 − x0)
c1 = (x1y0 − x0y1)/(x1 − x0)
Now suppose we want to determine the quadratic polynomial p(x) = c0 + c1x + c2x
2
that passes through
n = 3 data values. Then we have to solve the following set of linear equations for the polynomial coefficients
c:
1 x0 x
2
0
1 x1 x
2
1
1 x2 x
2
2
c0
c1
c2
=
y0
y1
y2
These are examples of second and third order Vandermonde Equations. Row i of the Vandermonde
matrix is formed by a sequence of powers of the i-th x data value.
You should be able to see that, for any collection of abscissas and ordinates, it is possible to define a
linear system whose solution yields the desired polynomial coefficients.
Assuming we have n data values xdata and ydata, the matrix A can be defined as follows:
A = np . z e r o s ( [ n , n ] )
fo r i in range ( 0 , n ) :
fo r j in range ( 0 , n ) :
A[ i , j ] = xdata [ i ] ∗ ∗ j
The right hand side of the linear system is simply the y data values. Thus, to solve for the polynomial
coefficients, we write:
c = np . l i n a l g . s o l v e ( A, ydata )
Once we have the coefficients c, we can evaluate our polynomial using polyval().
4
5 Exercise 2
We will use the Vandermonde matrix to find the coefficients of the polynomial that passes through a given
set of data points.
1. Write a code, coef vander.py that accepts a pair of vectors xdata and ydata and computes the coefficients c of the interpolating polynomial using the Vandermonde matrix. Your function should have
the signature:
def c o e f v a n d e r ( xdata , ydata ) :
””” c = c o e f v a n d e r ( xda ta , y d a t a )
. . . comments
”””
? ? ?
return c
Think carefully about how to determine ndata, the length of xdata, how this relates to the degree
degree of the polynomial, and which value to use when dimensioning the Vandermonde matrix.
2. Computing the coefficients of the polynomial through the following data. This polynomial is y = x
2
,
so you can check your coefficient vector by inspection.
xdata = np . a r r a y ( [ 0 , 1 , 2 ] )
ydata = np . a r r a y ( [ 0 , 1 , 4 ] )
3. Compute the coefficients of the polynomial through this data:
xdata = np . a r r a y ( [ −3, −2, −1, 0 , 1 , 2 , 3 ] )
ydata = np . a r r a y ( [ 1 6 3 6 , 2 4 7 , 2 8 , 7 , 4 , 3 1 , 412 ] )
4. Import polyval() using the statement from numpy.polynomial.polynomial import polyval and
confirm that your polynomial passes through these data points.
5. Using the np.polyfit() function to recompute the values of c, and comparing with your results.
Please include both the full command you used and the coefficient vector it returned in your summary.
6 Exercise 3
Now you will construct a polynomial using your code coef vander() to interpolate data points and then
you will see what happens between the interpolation points.
1. Import the polynomial library using the statement
import numpy . p ol yn omi al . p ol yn omi al a s p ol yn omi al
Use polyfromroots() to find the coefficients of the polynomial whose roots are [-2, -1, 1, 2, 3].
Call these coefficients cTrue.
2. This polynomial obviously passes through zero at each of these five “data” points. We want to see
if our coef vander() function can reproduce it. To use our coef vander function, we need a sixth
point. You can “read off” the value of the polynomial at x=0 from its coefficients cTrue. What is this
value?
3. Set up arrays xdata and ydata which list the polynomial roots, and the value at 0 (indicated here by
???).
5
xdata = np . a r r a y ( [ −2, −1, 0 , 1 , 2 , 3 ] )
ydata = np . a r r a y ( [ 0 , 0 , ? ? ? , 0 , 0 , 0 ] )
Use your coef vander() function to find the coefficients of this polynomial, and call these coefficients
cVander.
4. Use the following code to compute and plot the values of the true and interpolant polynomials on the
interval [-2,3]. Print err,the maximum difference between the two curves.
def t e s t p l o t ( cTrue , cVander ) :
import numpy . p ol yn omi al . p ol yn omi al a s p ol yn omi al
import m a t pl o tli b . p y pl o t a s p l t
import numpy a s np
x v al = np . l i n s p a c e ( −2.0 , 3 . 0 , 101 )
yvalTrue = p ol yn omi al . p ol y v al ( xval , cTrue )
yvalVander = p ol yn omi al . p ol y v al ( xval , cVander )
e r r = max ( abs ( ( yvalTrue − yvalVander ) ) ) / max ( abs ( yvalTrue ) )
p l t . pl o t ( xval , yvalTrue , ’ b−’ )
p l t . pl o t ( xval , yvalVander , ’ r o ’ )
p l t . s a v e f i g ( ’ t e s t p l o t . jp g ’ )
p l t . show ( )
return e r r
Please include a copy of this plot with your summary.
7 Lagrange Polynomials
Suppose we fix the set of N distinct abscissas xk, and think about the problem of constructing a polynomial
that has (not yet specified) values yk at these points. Now suppose there is a polynomial `7(x) whose value
is zero at each xk, k 6= 7, and is 1 at x7. Then the polynomial y7 · `7(x) would have the value y7 at x7,
and be 0 at all the other xi
. Doing the same for each abscissa and adding these intermediate polynomials
together results in the polynomial that interpolates the data, and we didn’t have to solve any equations to
do this.
In fact, the Lagrange polynomials `k are easily constructed for any set of N abscissas. There will be N
Lagrange polynomials, one per abscissa, each of degree N − 1, and the k
th polynomial `k(x) will have a
special relationship with the abscissa xk, namely, it will be 1 there, and 0 at the others.
This means that we can represent the interpolating polynomial to a set of data as:
p(x) = y1`1(x) + y2`2(x) + . . . + yN `N (x) (2)
=
X
N
k=1
yk `k(x)
This is a second way to define the interpolating polynomial to a set of data.
To construct the k-th polynomial, `k(x), which is 1 at the k-th data abscissa, and zero at the others, we
write:
`k(x) = f1(x) f2(x)· · · fk−1(x)· · · fk+1(x)· · · fN (x) (3)
skipping the k
th factor, where each factor has the form
fj (x) = (x − xj )
(xk − xj )
. (4)
6
Now we will show how to combine these Lagrange polynomials into a polynomial that interpolates a
given set of data.
8 Exercise 4
You will construct Lagrange polynomials based on the data set for y = x
2
:
xdata = np . a r r a y ( [ 0 1 2 ] )
ydata = np . a r r a y ( [ 0 1 4 ] )
1. Explain in a sentence or two why the formula (3) using the factors (4) yields a function
`k(xj ) =
1 j = k
0 j 6= k
(5)
2. Write a code lagrangep.py which defines the k−th Lagrange polynomial (3) associated with the abscissas
xdata and evaluates it at the point xval. Assume that xval may actually be a vector of length nval,
so that the result pval is also a vector of length nval:
def l a g r a n g e p ( k , xdata , x v al ) :
””” p v a l = l a g r a n g e p ( k , xda ta , x v a l )
. . . comments . . .
”””
x v al = np . a t l e a s t 1 d ( x v al ) # ( So we can t a k e l e n g t h even i f x v a l i s s c a l a r )
n v al = ? # l e n g t h o f x v a l
p v al = np . one s ( n v al ) # s t a r t each p r o d uc t as 1
ndata = ? # l e n g t h o f x d a t a
fo r j in range ( 0 , ndata ) :
i f ( j != k )
p v al = p v al ∗ ? ? ? # m u l t i p l y by t h e j−t h f a c t o r , b u t s k i p k
return p v al
3. From the definition in (5), what do you expect for the values of lagrangep( 0, xdata, xdata[0]),
lagrangep( 0, xdata, xdata[1]), lagrangep( 0, xdata, xdata[2])? Does your code agree with
your expectation?
4. Your code lagrangep() can evaluate a Lagrange polynomial at multiple points. So let xval be set to
xdata. What results do you compute for lagrangep( 0, xdata, xdata)? lagrangep( 1, xdata,
xdata)? lagrangep( 2, xdata, xdata)?
9 Exercise 5
Now we want to use our lagrangep() function as a helper for a second replacement for the polyfit() and
polyval() pair, called eval lag(), that implements Equation (2). Unlike coef vander(), the coefficient
vector of the polynomial does not need to be generated separately because it is so easy, and that is why
eval lag() both fits and evaluates the Lagrange interpolating polynomial.
1. Write a code eval lag.py with the signature
def e v a l l a g ( xdata , ydata , x v al ) :
””” y v a l = e v a l l a g ( xda ta , yda ta , x v a l )
7
. . . comments
”””
x v al = np . a t l e a s t 1 d ( x v al ) # ( So we can t a k e l e n g t h even i f x v a l i s s c a l a r )
n v al = ?
y v al = np . z e r o s ( n v al )
ndata = ?
fo r k in range ( 0 , ndata ) :
y v al = y v al + . . .
return y v al
This function should take the data values xdata and ydata, and compute the value of the interpolating
polynomial at xval according to (2), using your lagrangep() function for the Lagrange polynomials.
2. Test eval lag() on the simplest data set we have been using.
xdata = np . a r r a y ( [ 0 , 1 , 2 ] )
ydata = np . a r r a y ( [ 0 , 1 , 4 ] )
by evaluating it at xval=xdata. Of course, you should get ydata back.
3. Test your function by interpolating the polynomial that passes through the following points, again by
evaluating it at xval=xdata.
xdata = np . a r r a y ( [ −3, −2, −1, 0 , 1 , 2 , 3 ] )
ydata = np . a r r a y ( [ 1 6 3 6 , 2 4 7 , 2 8 , 7 , 4 , 3 1 , 412 ] )
4. Consider the function f(x) = tan−1
(40x − 15), over the interval [0,1]. Write a Python file hill.py
which evaluates this function. Let xdata be 7 equally spaced points in the interval, define ydata =
f(xdata), and let l(x) be the resulting Lagrange interpolant to this data.
Using 101 sample points in the interval [0,1], evaluate f(x) and l(x). Compute and print the maximum
absolute value of their difference. Instead of being tiny, this difference will probably be substantial.
Can you explain why?
Plot f(x) and l(x) together. Also display the 7 data values, where the two functions should agree,
using a command like
p l t . pl o t ( xdata , ydata , ’ r o ’ , m a r k e r si z e = 15 )
Does your Lagrange function interpolate f(x) well? Does it approximate f(x) well?
10 Exercise 6
When the function to be interpolated is not a polynomial, then away from the data points, the function and
interpolant will disagree. We usually have the choice of increasing the number of interpolation points, but
is this likely or guaranteed to make the interpolant closer to the function?
You will construct interpolants for the hyperbolic sine function sinh(x) and see that it and its polynomial
interpolant are quite close.
1. We would like to interpolate the function y = sinh(x) on the interval [−π, π], so use the following
outline to examine the behavior of the polynomial interpolant for five evenly-spaced points. Put these
commands into a file exercise6.py.
8
import m a t pl o tli b . p y pl o t a s p l t
import numpy a s np
# c o n s t r u c t da ta p o i n t s
ndata = 5
xdata = np . l i n s p a c e ( −np . pi , np . pi , ndata )
ydata = np . si n h ( xdata )
# c o n s t r u c t many t e s t p o i n t s
n v al = 4001
x v al = np . l i n s p a c e ( ? ? ? , ? ? ? , n v al )
# c o n s t r u c t t h e t r u e t e s t p o i n t v al u e s , f o r r e f e r e n c e
yvalTrue = np . si n h ( ? ? ? )
# Ev al u a te t h e i n t e r p o l a n t a t t h e t e s t p o i n t s
y v al = e v a l l a g ( ? ? ? , ? ? ? , x v al )
# e s t im a t e t h e r e l a t i v e a p pr ox im a t i on e r r o r o f t h e i n t e r p o l a n t
e r r = max ( abs ( yvalTrue − y v al ) ) / max ( abs ( yvalTrue ) )
# p l o t r e f e r e n c e v a l u e s and i n t e r p o l a t i o n da ta .
p l t . pl o t ( xval , yvalTrue , ’ g ’ , li n e wi d t h = 4 )
p l t . pl o t ( xdata , ydata , ’ r o ’ , m a r k e r si z e = 15 )
p l t . pl o t ( xval , yval , ’ k ’ )
p l t . s a v e f i g ( ’ e x e r c i s e 6 . jp g ’ )
p l t . show ( )
Please send me this plot.
2. By zooming, confirm visually that the function and its interpolant agree at the interpolation points.
3. Increasing ndata gives higher degree interpolation polynomials. Does this reduce the interpolation
error? Fill in the following table.
ndata Err
5 ________
11 ________
21 ________
You should observe that the approximation error becomes quite small as ndata is increased. Technically, all polynomials, as well as sinh(), are entire functions, which turns out to mean they can be well
approximated with our approach. In the following exercise, you will interpolate a function that is not entire.
11 Exercise 7
1. Write a file runge.py which evaluates the Runge function y =
1
1+x2 .
2. Copy exercise6.py to exercise7.py and modify it to use the Runge function, over the interval [-4,+4],
starting as before with ndata=5. Examine the plot, and verify that the interpolant and the Runge
function agree at the data points. Please send me the plot it generates.
3. Using more data points gives higher degree interpolation polynomials. Let’s see if higher degree means
lower error. Fill in the following table.
ndata Err
5 ________
9
11 ________
21 ________
Many people expect that an interpolating polynomial p(x) gives a good approximation to the function
f(x) everywhere, no matter what function we choose. If the approximation is not good, we expect it to get
better if we increase the number of data points. These expectations will be fulfilled only when the function
does not exhibit some “essentially non-polynomial” behavior. You will see why the Runge function cannot
be approximated well by polynomials.
12 Exercise 8
The Runge function has Taylor series
1
1 + x
2
= 1 − x
2 + x
4 − x
6 + . . . (6)
which has the (infinite) coefficient vector [1,0,-1,0,1,0,-1,0,…].
This series has a radius of convergence of 1 in the complex plane. Polynomials, on the other hand, are
entire functions, meaning that their Taylor series converge everywhere in the complex plane. Any finite sum
of polynomials must be entire, and no entire function can interpolate the Runge function on a disk with
radius larger than one about the origin in the complex plane. If there were one, it would have to agree with
the series (6) inside the unit disk, but the series diverges at x = i and an entire function cannot have an
infinite value (a pole).
1. Make a copy of exercise7.py called exercise8.py that uses coef vander() and polyval() to evaluate
the interpolating polynomial rather than eval lag().
2. Look at the tendency of the coefficients ck of the interpolating polynomials by filling in the following
table.
ndata=5 c[ 0]= _____ c[ 2]= _____ c[ 4]= _____
ndata=11 c[ 0]= _____ c[ 2]= _____ c[ 4]= _____ c[ 6]= _____ c[ 8]= _____ c[10]= _____
ndata=21 c[ 0]= _____ c[ 2]= _____ c[ 4]= _____ c[ 6]= _____ c[ 8]= _____ c[10]= _____
Taylor +1 -1 +1 -1 +1 -1
You should see that the interpolating polynomials are “trying” to reproduce the Taylor series (6). These
polynomials cannot agree with the Taylor series at all points, though, because the Taylor series does not
converge at all points.
13 Trigonometric interpolation
Instead of using powers of x to build our approximation, trigonometic interpolation uses combinations of
sines and cosines for the interpolation. Mathematically, these combinations are efficiently expressed using
exponential format, using the relationship
e
iθ = cos θ + isin θ
For the moment, we will use the mathematical convention that the complex unit is denoted by i. When we
get into Python coding, we will have to change our conventions!
10
We plan to approximate a function f(x) using a trigonometric interpolant t(x). The basic expression for
a trigonometric interpolant t(x), using (2n + 1) functions over the interval [−π, π], is
t(x) = X
2n
k=0
ak e
i(k−n)x
. (7)
We consider (2n + 1) points because we want to use real functions and real interpolants. Thus one trigonometric basis function e
i(k−n)x
is the constant term when k = n, and all the rest come in complex conjugate
pairs.
For data values, we will use 2n + 1 evenly-spaced points in the interval (−π, π). These data points are
offset from the endpoints, and so are given by:
xl o = − 2 ∗ n / ( 2 ∗ n + 1 ) ∗ np . pi # a l i t t l e more than −p i
xhi = + 2 ∗ n / ( 2 ∗ n + 1 ) ∗ np . pi # a l i t t l e l e s s than +p i
x = np . l i n s p a c e ( xlo , xhi , 2∗n+1 )
For k = 0, 1, . . . , 2n, the coefficient ak can be determined by multiplying the trigonometric basis functions,
with a negated exponent, by the original function f(x) at the x data values:
ak =
1
2n + 1
X
2n
j=0
f(xj ) e
−i(k−n)xj
(8)
The trigonometric coefficients ak in (8) play the same role as the polynomial coefficients ck in (1). Now
we are going to write functions coef trig() to be analogous to polyfit() and coef vander(), and also
eval trig() to be analogous to polyval() and eval lag().
14 Exercise 9
1. Write a code coef trig.py to evaluate the trigonometric coefficients, with the signature
def c o e f t r i g ( func , ndata ) :
””” a = c o e f t r i g ( func , n )
. . . comments . . .
”””
xl o = − 2 ∗ n / ( 2 ∗ n + 1 ) ∗ np . pi # a l i t t l e more than −p i
xhi = + 2 ∗ n / ( 2 ∗ n + 1 ) ∗ np . pi # a l i t t l e l e s s than +p i
xdata = np . l i n s p a c e ( xlo , xhi , 2∗n+1 )
a = ? ? ?
return a
Replace the ??? text by defining a for loop that carries out the summation according to Equation (8).
Note that a mathematical quantity like value = e
−iθ is expressed in Python by
1 v al u e = np . exp ( − 1 j ∗ t h e t a )
Here 1j represents the imaginary unit in Python.
2. Test coef trig() by applying it to the function f(x) = e
ix, for n=10. You can write a separate
function file for f(x), or else use a lambda expression directly in the call to coef trig(), like this:
a = c o e f t r i g ( lambda x : np . exp ( − 1 j ∗ x ) , n )
By examining Equation (7), you should be able to see that ak = 1.0 for k = n + 1 and zero otherwise.
11
3. Test coef trig() again by applying it to f(x) = sin 4x with n=10. Again, you can create a function
file for this f(x), or use a lambda expression. You should see that ak is zero for all but two subscripts.
What subscripts (k) have non-zero ak and what are the values?
15 Exercise 10
We can compute the coefficients a of a trigonometric interpolant t(x) for a given function f(x). Now we
want to be able to evaluate t(x), plot it, or compare it to the original f(x).
1. Write a file eval trig.py to evaluate trigonometric interpolants, with the signature
def e v a l t r i g ( a , x v al ) :
””” t v a l = e v a l t r i g ( a , x v a l )
. . . comments . . .
”””
n v al = len ( x v al )
n = int ( ( len ( a ) − 1 ) / 2 )
t v a l = np . z e r o s ( n v al )
fo r k in range ( 0 , 2 ∗ n + 1 ) :
? ? ?
return t v a l
This function should evaluate Equation (7) at an arbitrary collection of points, xval.
2. Write a file exercise10.py which uses coef trig() to find the coefficients for sin 4x with n=10 and then
applies eval trig() to evaluate t(x) at 101 equally-spaced points in the interval [−π, π]. Compare
t(x) to f(x) = sin 4x by plotting them together: the lines should overlap. Your interpolated values
may appear to to be complex; if so, they might not be plotted the way you expect. In that case, plot
only the real parts; that is, instead of plotting tval plot tval.real. Please send me the plot.
16 Exercise 11
Trigonometric polynomial interpolation can do well even when the functions are not themselves trigonometric
polynomials.
1. Make a new file exercise11.py by copying exercise10.py. Now interpolate the function f(x) = x(π
2−x
2
),
over the interval [−π, +π], using trigonometric interpolation with n=5, 10, 15, 20. For each case,
compute err, the maximum absolute value of the error. You do not have to send me the plots, but
submit the following error table.
y = x*(pi**2-x**2)
n Err
5: ________
10: ________
15: ________
20: ________
2. Do the same thing for the Runge function, completing and turning in the following table. (You do not
need to submit the plots).
12
y = 1/(1+x^2)
n Err
2: ________
4: ________
5: ________
17 Exercise 12
The previous exercise shows that trigonometric polynomial interpolation does well for some functions. It
does much less well when the function is not continuous or not periodic in [−π, π]. Consider, for example,
the function
f(x) = (
x + π for x < 0
x − π for x ≥ 0
1. Create a code sawshape6.py
def sawshape6 ( x ) :
””” y = sawshape6 ( x )
. . . comments . . .
”””
y = ( x + np . pi ) ∗ ( x < 0 ) + ( x − np . pi ) ∗ ( 0 <= x )
return y
2. Write a file exercise12.py by copying and adjusting exercise11.py, which carries out trigonometric
interpolation on the sawshape6 function over the interval [−π, π], and compute and print the maximum
relative interpolation error. Turn in the table:
y = sawshape6()
n Err
5 ________
10 ________
100 ________
Please send me the plots for n=5 and n=100.
For polynomial interpolation, increasing the order can lead to divergence of approximation. However,
trigonometric interpolation does not diverge as n → ∞. However, if the function to be interpolated has
discontinuities, convergence of the interpolant is greatly slowed down.
13
