Homework 0: Interpolation examples (Python version)¶

Developed by Randy LeVeque for a course on Approximation Theory and Spectral Methods at the University of Washington.

See http://faculty.washington.edu/rjl/classes/am590a2013/codes.html for more IPython Notebook examples.

This is a demonstration of using an IPython notebook to combine code with descriptions.

It also demonstrates how to use pychebfun from within the notebook.

Sample problem¶

Determine the quadratic polynomial $p(x)$ that interpolates the data

$(-1,3),~~(0,-1),~~(1,2)$ first using the Vandermonde matrix and then using the Lagrange form of interpolating polynomial.

Note that if the notebook is started with --pylab option then numpy and matplotlib are already imported.

Interpolation data for all the examples:¶

In [1]:

xj = array([-1., 0., 1.])
yj = array([3., -1., 2.])

Using monomials¶

For a quadratic interpolation we use basis functions $1,~x,$ and $x^2$ .

Define the Vandermonde matrix: The columns are the basis functions evaluated at the interpolation points.

In [2]:

A = vstack((xj**0, xj, xj**2)).T
print A

[[ 1. -1.  1.]
 [ 1.  0.  0.]
 [ 1.  1.  1.]]

Solve the system for the monomial coefficients:

In [3]:

from numpy.linalg import solve
c = solve(A,yj)
print c

[-1.  -0.5  3.5]

Plot the resulting polynomial on a fine grid:

In [4]:

x = linspace(-1, 1, 1001)
p = c[0] + c[1]*x + c[2]*x**2
plot(x,p)
hold(True)
plot(xj,yj,'b.',markersize=10)

Out[4]:

[<matplotlib.lines.Line2D at 0x10f0b0850>]

Lagrange form¶

Evaluate the Lagrange basis functions on the fine grid for plotting:

$\begin{align} L_1(x) &= \frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)},\\\\ L_2(x) &= \frac{(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x_3)},\\\\ L_3(x) &= \frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)} \end{align}$

Evaluate the Lagrange basis functions on the fine grid for plotting:

In [5]:

L1 = (x-xj[1])*(x-xj[2]) / ((xj[0]-xj[1])*(xj[0]-xj[2]))
L2 = (x-xj[0])*(x-xj[2]) / ((xj[1]-xj[0])*(xj[1]-xj[2]))
L3 = (x-xj[0])*(x-xj[1]) / ((xj[2]-xj[0])*(xj[2]-xj[1]))

In this form, the data values are the coefficients

In [6]:

p = yj[0]*L1 + yj[1]*L2 + yj[2]*L3

Plot this version (should look the same as before!)

In [7]:

clf()
plot(x,p)
hold(True)
plot(xj,yj,'b.',markersize=10)

Out[7]:

[<matplotlib.lines.Line2D at 0x10f14ed90>]

Using pychebfun¶

To use pychebfun, you must download or clone from https://github.com/olivierverdier/pychebfun and add the location to the python path. See the README file at the bottom of the github page.

In [8]:

import pychebfun as P

Use pychebfun to plot the function $f(x) = \left(\frac{1}{1 + 25x^2}\right)^{10}$ and the function $g(x) = 5 \int_0^x f(s)\, ds.$

What degree polynomial approximation is used for each?

In [9]:

x = P.Chebfun.identity()
f = (1/(1+25*x**2))**10;

In [10]:

f.plot()

g = 5 * f.integrate()
g = g - g(-1.)  # because f.integrate gives integral from 0 rather than from -1.
g.plot()

legend(['f','g'])
print "The degree of f is %i " % (len(f.chebyshev_coefficients()) - 1)
print "The degree of g is %i " % (len(g.chebyshev_coefficients()) - 1)

The degree of f is 234 
The degree of g is 201

See http://nbviewer.ipython.org/6724986 for more pychebfun examples.

In [10]: