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.
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.
xj = array([-1., 0., 1.])
yj = array([3., -1., 2.])
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.
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:
from numpy.linalg import solve
c = solve(A,yj)
print c
[-1. -0.5 3.5]
Plot the resulting polynomial on a fine grid:
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)
[<matplotlib.lines.Line2D at 0x10f0b0850>]
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:
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
p = yj[0]*L1 + yj[1]*L2 + yj[2]*L3
Plot this version (should look the same as before!)
clf()
plot(x,p)
hold(True)
plot(xj,yj,'b.',markersize=10)
[<matplotlib.lines.Line2D at 0x10f14ed90>]
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.
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?
x = P.Chebfun.identity()
f = (1/(1+25*x**2))**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.