## 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