Chebfun in Python

This is a demo following the original paper of Battles & Trefethen.

Initialisation

In principle, it suffices to use the class Chebfun, but we'll import all the objects in chebfun for convenience.

In [2]:
from pychebfun import *
In [3]:
x = Chebfun.identity()

This is only for this IPython notebook

In [85]:
from IPython.core.display import HTML

Chebfuns and barycentric interpolation

Creation of a chebfun

A typical way of initialising a Chebfun is to use the function chebfun. It will work for any function.

In [86]:
cube = chebfun(lambda x:x**3)

A slightly more pythonic way to achieve the same result would be using Chebfun.from_function, it will also work for any function.

In [5]:
cube = Chebfun.from_function(lambda x:x**3)

Another possibility is to use the variable x defined above by x = Chebfun.identity(). Note however, that this will only work if the function of x has been registered in chebfun. Here it works because chebfun knows how to compute arbitray powers.

In [6]:
cube = x**3

Other examples could include:

In [88]:
cos(x)
tan(x)
exp(sin(x))
Out[88]:
<Chebfun(array([ 2.31977682,  2.30256549,  2.25010157,  2.16063784,  2.03347888,
        1.8713815 ,  1.68216731,  1.47842802,  1.27510342,  1.08598295,
        0.92082477,  0.78425011,  0.67639411,  0.59447119,  0.53436459,
        0.49176808,  0.46282629,  0.44442438,  0.43429818,  0.43107595]))>

Display

The size of f is the number of interpolation point.

In [7]:
cube.size()
Out[7]:
4

Displaying a Chebfun gives the interpolation values at the Chebyshev interpolation points:

In [8]:
print repr(cube)
<Chebfun(array([ 1.   ,  0.125, -0.125, -1.   ]))>
In [10]:
Chebfun.from_function(lambda x: x**7 - 3*x + 5).size()
Out[10]:
8
In [11]:
f_sin5 = Chebfun.from_function(lambda x: sin(5*pi*x))
f_sin5.size()
Out[11]:
42

Convergence

For some functions, convergence is not achievable, and one must set a limit to the dichotomy algorithm.

In [12]:
f_abs = Chebfun.from_function(abs, N=1000)

Evaluations

Chebfun objects behave as ordinary function, which can be evaluated at any point.

In [15]:
cube(5)
Out[15]:
array(124.99999999999774)
In [16]:
cube(-.5)
Out[16]:
array(-0.1250000000000001)

It is also possible (and more efficient) to evaluate a chebfun at several point at once. For instance, we evaluate the Chebfun for $sin(5x)$ at several points at once:

In [17]:
f_sin5(linspace(0, .2, 5))
Out[17]:
array([  4.33257766e-16,   7.07106781e-01,   1.00000000e+00,
         7.07106781e-01,  -2.68873793e-15])

Operations

One can for instance add two chebfuns:

In [18]:
chebsum = cube+f_sin5
print chebsum.size()
chebsum.plot()
42

Or multiply them:

In [19]:
chebmult = cube*f_sin5
print chebmult.size()
chebmult.plot()
43
In [23]:
f_sin = Chebfun.from_function(lambda x:sin(x))
print f_sin.size()
print (f_sin*f_sin).size()
14
17

It is also possible to add and multiply chebfuns with scalars:

In [24]:
subplot(121)
(3*f_sin).plot()
subplot(122)
(3+f_sin).plot()
In [25]:
f_cubic = Chebfun.from_function(lambda x: x**3 - x)
f_cubic.plot()
In [26]:
f_sin5.plot()

2D Chebfun

In the paper, they suggest to create two Chebfuns as follows:

In [27]:
sin(16*x), sin(18*x)
Out[27]:
(<Chebfun(array([-0.28790332, -0.2426432 , -0.10433329,  0.12848036,  0.43822319,
        0.76298866,  0.97977108,  0.92715359,  0.49731325, -0.22071134,
       -0.86140282, -0.93924815, -0.27138199,  0.66515185,  0.97793397,
        0.25892868, -0.76748918, -0.8903491 ,  0.09519056,  0.96531778,
        0.57519601, -0.57519601, -0.96531778, -0.09519056,  0.8903491 ,
        0.76748918, -0.25892868, -0.97793397, -0.66515185,  0.27138199,
        0.93924815,  0.86140282,  0.22071134, -0.49731325, -0.92715359,
       -0.97977108, -0.76298866, -0.43822319, -0.12848036,  0.10433329,
        0.2426432 ,  0.28790332]))>,
 <Chebfun(array([-0.75098725, -0.78181709, -0.86309961, -0.95807035, -0.99912028,
       -0.89316158, -0.55644827,  0.0110737 ,  0.64115513,  0.99398134,
        0.74432845, -0.07297019, -0.86854365, -0.86676693,  0.04783357,
        0.93000121,  0.68652549, -0.45407882, -0.98729282, -0.12752248,
        0.92188239,  0.61105796, -0.61105796, -0.92188239,  0.12752248,
        0.98729282,  0.45407882, -0.68652549, -0.93000121, -0.04783357,
        0.86676693,  0.86854365,  0.07297019, -0.74432845, -0.99398134,
       -0.64115513, -0.0110737 ,  0.55644827,  0.89316158,  0.99912028,
        0.95807035,  0.86309961,  0.78181709,  0.75098725]))>)

It is certainly possible, but we can't create a 2D Chebfun from them, at the moment. Instead we have to do it like this:

In [29]:
def vector_func(x):
    return vstack([sin(16*x), sin(18*x)]).T
cheb_vec = Chebfun.from_function(vector_func)
cheb_vec.plot()

Elementary functions

Again, we explore the two ways to create a chebfun. First, the general way, which works for any smooth functions:

In [30]:
f_expsin = Chebfun.from_function(lambda x: exp(sin(x)))

Then, the “operator overloading” way. We take advantage of the fact that exp and sin are defined on chebfuns. Not all functions are, though. This is very similar to ufuncs in numpy.

In [32]:
g_expsin = exp(sin(x))

Of course, both chebfuns are equivalent, as we demonstrate here.

In [35]:
print f_expsin.size()
print g_expsin.size()
20
20
In [36]:
subplot(121)
f_expsin.plot()
subplot(122)
g_expsin.plot()
In [37]:
(f_expsin - g_expsin).norm()
Out[37]:
3.1401849173675503e-16

Approximation Theory

Gibbs phenomenon

By limiting the accuracy of the chebfun, one observes the celebrated Gibbs phenomenon.

In [38]:
f_sign = Chebfun.from_function(lambda x: sign(x), N=25)
f_sign.plot()

Pychebfun implements the method compare which plots a graph of the error.

In [40]:
def absfive(x):
    return abs(x)**5
#errors = array([(Chebfun(absfive, N) - exact).norm() for N in range(10, 20)])
#loglog(errors)
chebfun(absfive).compare(absfive)
Out[40]:
<matplotlib.axes.AxesSubplot at 0x106b88d50>

Interpolation of random data

A chebfun passing through random values at 50 Chebyshev points. The method to initialise a Chebfun from data at Chebyshev points is from_data:

In [50]:
rand_chebfun = Chebfun.from_data(rand(51))
rand_chebfun.plot()
rand_chebfun.size()
Out[50]:
51
In [51]:
f_walk = Chebfun.from_data(randn(100,2))
f_walk.plot()

Extrapolation

Clearly, as chebfuns are based on polynomial interpolation, there is no guarantee that extrapolation outside the interval $[-1,1]$ will give any sensible result at all. This is an illustration of that phenomenon.

In [53]:
xx = linspace(-1.4, 1.4, 100)
error = log10(abs(f_sin5(xx)-sin(5*pi*xx)))
subplot(211)
plot(xx, f_sin5(xx), marker='.')
ylabel('interpolant')
subplot(212)
plot(xx,error, marker='.')
ylabel('log_10(error)')
Out[53]:
<matplotlib.text.Text at 0x106c08cd0>

Chebyshev expansions and FFT

The Chebyshev coefficients of a chebfun are quickly computed using a fast Fourier transform. The corresponding method is chebyshev_coefficients

In [54]:
chebfun(exp(x)).chebyshev_coefficients()
Out[54]:
array([  1.26606588e+00,   1.13031821e+00,   2.71495340e-01,
         4.43368498e-02,   5.47424044e-03,   5.42926312e-04,
         4.49773230e-05,   3.19843646e-06,   1.99212481e-07,
         1.10367718e-08,   5.50589647e-10,   2.49795740e-11,
         1.03930539e-12,   4.00363503e-14])
In [55]:
Chebfun.from_chebcoeff(array([3,2,1.]))
Out[55]:
<Chebfun(array([ 6.,  2.,  2.]))>

One way to create the basis Chebyshev polynomial of order 10.

In [57]:
T_10 = Chebfun.from_chebcoeff(array([0.]*10+[1]))
T_10.plot()

Incidentally, the same result is achieved using Chebfun.basis:

In [58]:
T_10_ = Chebfun.basis(10).plot()

As an illustration of how fast the fast Fourier transform is, we create a chebfun with 100000 points.

In [59]:
r = randn(100000)
f_randn = Chebfun.from_chebcoeff(r)
In [60]:
norm(r)
Out[60]:
315.72498126138117
In [61]:
f_randn.sum()
#f_randn.norm() # doesn't work yet
Out[61]:
2.7094293040817283
In [62]:
T_20 = Chebfun.basis(20)
T_20.norm()
Out[62]:
1.4142135623730951

Integration and Differentiation

Clenshaw Curtis Quadrature

The method sum computes the integral of the chebfun using the Chebyshev expansion and integral of the Chebyshev basis polynomial (it is known as the Clenshaw–Curtis quadrature)

In [66]:
k_odd = 5
k_even = 2
HTML(r'For odd $k$ we check that $\int T_k = %s$ ($k=%s$), otherwise, $\int T_k = \frac{2}{1-k^2} = %s \approx %s$ for $k = %s$.' % 
(Chebfun.basis(k_odd).sum(), 
 k_odd,
2./(1-k_even**2), 
Chebfun.basis(k_even).sum(), 
k_even))
Out[66]:
For odd $k$ we check that $\int T_k = 0.0$ ($k=5$), otherwise, $\int T_k = \frac{2}{1-k^2} = -0.666666666667 \approx -0.666666666667$ for $k = 2$.

Examples

Computing some integrals.

In [49]:
f1 = sin(pi*x)**2
print f1.size(), f1.sum()
HTML('$\int \sin(\pi x)^2 dx = %s$' % f1.sum())
25 1.0
Out[49]:
$\int \sin(\pi x)^2 dx = 1.0$
In [51]:
f2 = chebfun(lambda x: 1/(5+3*cos(pi*x)))
print f2.size(), f2.sum()
51 0.5
In [53]:
f3 = chebfun(lambda x: abs(x)**9 * log(abs(x)+1e-100))
print f3.size(), f3.sum()
91 -0.02
In [54]:
HTML(r'Computing the norm of $x^2$ gives: %s $\approx \sqrt{\int_{-1}^1 x^4 dx} = \sqrt{\frac{2}{5}} = %s$' % ((x**2).norm(), sqrt(2/5)))
Out[54]:
Computing the norm of $x^2$ gives: 0.632455532034 $\approx \sqrt{\int_{-1}^1 x^4 dx} = \sqrt{\frac{2}{5}} = 0.632455532034$

The dot method computes the Hilbert scalar product, so f.dot(g) corresponds to $$\int_{-1}^1 f(x) g(x) \mathrm{d} x$$

In [67]:
x.dot(x)
Out[67]:
0.66666666666666674

The integrate method computes a primitive $F$ of a chebfun $f$, which is zero at zero i.e., $$F(x) = \int_{0}^x f(y) \mathrm{d}y$$

In [57]:
f = chebfun(lambda t: 2/sqrt(pi) * exp(-t**2))
erf2_ = f.integrate()
erf2 = erf2_ - erf2_(0)
from scipy.special import erf
randx = rand()
#print erf2(randx) - erf(randx)

This allows to define continuous versions of the prod and cumprod functions for chebfuns:

In [69]:
def prod(f):
    return exp(log(f).sum())
def cumprod(f):
    return exp(log(f).integrate())
In [70]:
prod(exp(x))
Out[70]:
0.99999999999999911
In [71]:
prod(exp(exp(x)))
Out[71]:
10.489789833690235
In [72]:
cumprod(exp(x)).plot()

Differentiation

One can also differentiate chebfuns to arbitrary orders with the method differentiate

In [62]:
f_sin5 = sin(5*x)
fd_sin5 = f_sin5.differentiate(4)
print fd_sin5.norm()/f_sin5.norm()
625.000000029
In [63]:
f = sin(exp(x**2)).differentiate()
print f(1)
-4.95669946592
In [73]:
g = Chebfun.from_function(lambda x: 1/(2 + x**2))
h = g.differentiate()
print h(1)
-0.222222222229
In [65]:
#f.differentiate().integrate() == f - f(-1)
#f.integrate().differentiate() == f

Operations based on rootfinding

As chebfun are polynomial it is possible to find all the roots in the interval $[-1,1]$ using the method roots:

In [74]:
print (x - cos(x)).roots()
[ 0.73908513]
In [75]:
print (x - cos(4*x)).roots()
[-0.89882622 -0.53333306  0.31308831]

Zeros of the Bessel function $J_0$ on the interval $[0,20]$:

In [76]:
from scipy.special import jn
def J_0(x):
    return jn(0,x)
f = chebfun(lambda x: J_0(10*(1+x)))
f.plot()
froots = f.roots()
roots = 10*(1+froots)
print 'roots:', roots
plot(froots, zeros_like(froots), color='red', marker='.', linestyle='', markersize=10)
grid()
roots: [  2.40482556   5.52007811   8.65372791  11.79153444  14.93091771
  18.07106397]

Extrema

The methods min and max are not implemented but it is not difficult to do so, by finding the roots of the derivative. This may come in a future version of pychebfun.

In [77]:
f = x - x**2
#print 'min', f.min()
#print 'max', f.max()
#print 'argmax', f.argmax()
#print 'norm inf', f.norm('inf')
#print 'norm 1', f.norm(1)

Total variation

In [78]:
def total_variation(f):
    return f.differentiate().norm(1)
#total_variation(x)
#total_variation(sin(5*pi*x))

Here is an example of computation of extrema of a function.

In [82]:
r = f.roots()
fd = f.differentiate()
fdd = fd.differentiate()
e = fd.roots()
ma = e[fdd(e) <= 0]
mi = e[fdd(e) > 0]
def plot_all():
    f.plot(with_interpolation_points=False)
    plot(r, zeros_like(r), linestyle='', marker='o', color='white', markersize=8)
    plot(ma, f(ma), linestyle='', marker='o', color='green', markersize=8)
    plot(mi, f(mi), linestyle='', marker='o', color='red', markersize=8)
    grid()
plot_all()
figure()
plot_all()
axis([-.1, .1,-.7,-.6])
Out[82]:
[-0.1, 0.1, -0.7, -0.6]

Applications in numerical analysis

Quadrature

We apply quadrature to the function $$ f(x) = \tan(x+\frac{1}{4}) + \cos(10x^2 + \exp(\exp(x)))$$

In [79]:
def s(x):
    return tan(x+1/4) + cos(10*x**2 + exp(exp(x)))
f = chebfun(s)
f.plot()
title(str(f))
Out[79]:
<matplotlib.text.Text at 0x106f13b90>

Using the built-in quad quadrature integrator of scipy.integrate:

In [80]:
import scipy.integrate
scipy.integrate.quad(s, -1, 1, epsabs=1e-14)[0]
Out[80]:
0.2954776762437714

Using Chebyshev integration, we obtain exactly the same value:

In [81]:
f.sum()
Out[81]:
0.2954776762437713

ODE Solving

One solves the ODE $$u'(x) = \mathrm{e}^{-2.75 x u}$$ $$u(-1)=0$$

The problem is rewritten first as an integral problem, namely $$u(x) = T(u)$$ with $$T(u) := \int_{-1}^{x} \mathrm{e}^{-2.75 yu(y)} dy$$ This suggest using the following Picard iteration:

In [83]:
uold = Chebfun(0.)
du = 1.
for i in xrange(100):
    integrand = exp(-2.75*x*uold)
    uint = integrand.integrate()
    u = uint - uint(-1)
    du = u-uold
    uold = u
    if du.norm() < 1e-13:
        break
else:
    print 'no convergence'
print 'Converged in {} steps'.format(i)
print 'u(1) = {:.14f}'.format(float(u(1)))
u.plot(with_interpolation_points=False)
grid()
title(str(u))
Converged in 32 steps
u(1) = 5.07818302387461
Out[83]:
<matplotlib.text.Text at 0x106f0ac50>

Boundary value problems

One solves the boundary problem $$u''(x) = \mathrm{e}^{4x}$$ with boundary conditions $$u(-1) = 0 \qquad u(1) = 0 $$

In [76]:
f = exp(4*x)
u_ = f.integrate().integrate()
u = u_ - (u_(1)*(x+1) + u_(-1)*(1-x))/2
u.plot()
grid()