Chebfun in Python¶

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

Initialisation¶

In [1]:

%matplotlib inline
%config InlineBackend.figure_formats = {'svg', 'png'}
import numpy as np
import matplotlib.pyplot as plt

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 Chebfun

In [3]:

import pychebfun

In [4]:

x = Chebfun.identity()

This is only for this IPython notebook

In [5]:

from IPython.core.display import HTML

We define a convenience plotting function

In [6]:

def cplot(fun):
    pychebfun.plot(fun, with_interpolation_points=True)

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 [7]:

cube = pychebfun.chebfun(lambda x:x**3)

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

In [8]:

cube = Chebfun.from_function(lambda x:x**3)

Another possibility is to use the variable x defined by x = Chebfun.identity():

In [9]:

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 [10]:

cube = x**3

Other examples could include:

In [11]:

print(np.cos(x))
print(np.tan(x))
print(np.exp(np.sin(x)))

<Chebfun(13)>
<Chebfun(32)>
<Chebfun(20)>

Display¶

The size of f is the number of interpolation point.

In [12]:

cube.size()

Out[12]:

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

In [13]:

print(repr(cube))

Chebfun 
     domain        length     endpoint values
  [ -1.0,   1.0]         4       -1.00    1.00
 vscale = 1.00e+00

In [14]:

Chebfun.from_function(lambda x: x**7 - 3*x + 5).size()

Out[14]:

In [15]:

f_sin5 = Chebfun.from_function(lambda x: np.sin(5*np.pi*x))
f_sin5.size()

Out[15]:

Convergence¶

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

In [16]:

f_abs = Chebfun.from_function(abs, N=1000)

Evaluations¶

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

In [17]:

cube(5)

Out[17]:

array(125.)

In [18]:

cube(-.5)

Out[18]:

array(-0.125)

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 [19]:

f_sin5(np.linspace(0, .2, 5))

Out[19]:

array([ 3.17383112e-16,  7.07106781e-01,  1.00000000e+00,  7.07106781e-01,
       -1.26628627e-15])

Operations¶

One can for instance add two chebfuns:

In [20]:

chebsum = cube+f_sin5
print(chebsum.size())
cplot(chebsum)

Or multiply them:

In [21]:

chebmult = cube*f_sin5
print(chebmult.size())
cplot(chebmult)

In [22]:

f_sin = Chebfun.from_function(lambda x:np.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 [23]:

plt.subplot(121)
cplot(3*f_sin)
plt.subplot(122)
cplot(3+f_sin)

In [24]:

f_cubic = Chebfun.from_function(lambda x: x**3 - x)
cplot(f_cubic)

In [25]:

cplot(f_sin5)

2D Chebfun¶

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

In [26]:

np.sin(16*x), np.sin(18*x)

Out[26]:

(Chebfun 
      domain        length     endpoint values
   [ -1.0,   1.0]        42        0.29   -0.29
  vscale = 1.00e+00,
 Chebfun 
      domain        length     endpoint values
   [ -1.0,   1.0]        44        0.75   -0.75
  vscale = 1.00e+00)

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 [27]:

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

Elementary functions¶

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

In [28]:

f_expsin = Chebfun.from_function(lambda x: np.exp(np.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 [29]:

g_expsin = np.exp(np.sin(x))

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

In [30]:

print(f_expsin.size())
print(g_expsin.size())

20
20

In [31]:

plt.subplot(121)
cplot(f_expsin)
plt.subplot(122)
cplot(g_expsin)

In [32]:

(f_expsin - g_expsin).norm()

Out[32]:

6.280369834735101e-16

Approximation Theory¶

Gibbs phenomenon¶

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

In [33]:

f_sign = Chebfun.from_function(lambda x: np.sign(x), N=25)
cplot(f_sign)

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

In [34]:

def absfive(x):
    return np.abs(x)**5
#errors = array([(Chebfun(absfive, N) - exact).norm() for N in range(10, 20)])
#loglog(errors)
pychebfun.compare(Chebfun.from_function(absfive), absfive)

Out[34]:

<Axes: >

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 [35]:

rng = np.random.default_rng(0)

In [36]:

rand_chebfun = Chebfun.from_data(rng.random(51))
cplot(rand_chebfun)
rand_chebfun.size()

Out[36]:

In [37]:

f_walk = Chebfun.from_data(rng.normal(size=(100,2)))
cplot(f_walk)

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 [38]:

xx = np.linspace(-1.4, 1.4, 100)
error = np.log10(np.abs(f_sin5(xx)-np.sin(5*np.pi*xx))+1e-16)
plt.subplot(211)
plt.plot(xx, f_sin5(xx), marker='.')
plt.ylabel('interpolant')
plt.subplot(212)
plt.plot(xx,error, marker='.')
plt.ylabel('log_10(error)')

Out[38]:

Text(0, 0.5, 'log_10(error)')

Chebyshev expansions and FFT¶

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

In [39]:

Chebfun.from_function(np.exp(x)).coefficients()

Out[39]:

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.10367720e-08, 5.50589728e-10, 2.49795910e-11,
       1.03889547e-12, 3.97630646e-14])

In [40]:

Chebfun.from_coeff(np.array([3,2,1.]))

Out[40]:

Chebfun 
     domain        length     endpoint values
  [ -1.0,   1.0]         3        2.00    6.00
 vscale = 1.00e+00

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

In [41]:

T_10 = Chebfun.from_coeff(np.array([0.]*10+[1]))
cplot(T_10)

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

In [42]:

cplot(Chebfun.basis(10))

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

In [43]:

r = rng.normal(size=100000)
f_randn = Chebfun.from_coeff(r)

In [44]:

np.sqrt(np.sum(np.square(r)))

Out[44]:

316.21969834936846

In [45]:

f_randn.sum()
#f_randn.norm() # doesn't work yet

Out[45]:

0.35631543750789174

In [46]:

T_20 = Chebfun.basis(20)
T_20.norm()

Out[46]:

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 [47]:

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[47]:

For odd

$k$ we check that

$\int T_k = 0.0$ (

$k=5$ ), otherwise,

$\int T_k = \frac{2}{1-k^2} = -0.6666666666666666 \approx -0.6666666666666666$ for

$k = 2$ .

Examples¶

Computing some integrals.

In [48]:

f1 = np.sin(np.pi*x)**2
print(f1.size(), f1.sum())
HTML(r'$\int \sin(\pi x)^2 dx = %s$' % f1.sum())

25 0.9999999999999993

Out[48]:

$\int \sin(\pi x)^2 dx = 0.9999999999999993$

In [49]:

f2 = Chebfun.from_function(lambda x: 1/(5+3*np.cos(np.pi*x)))
print(f2.size(), f2.sum())

51 0.5

In [50]:

f3 = Chebfun.from_function(lambda x: np.abs(x)**9 * np.log(np.abs(x)+1e-100))
print(f3.size(), f3.sum())

91 -0.019999999999999997

In [51]:

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(), np.sqrt(2./5)))

Out[51]:

Computing the norm of

$x^2$ gives: 0.6324555320336759

$\approx \sqrt{\int_{-1}^1 x^4 dx} = \sqrt{\frac{2}{5}} = 0.6324555320336759$

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 [52]:

x.dot(x)

Out[52]:

0.6666666666666667

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 [53]:

f = Chebfun.from_function(lambda t: 2/np.sqrt(np.pi) * np.exp(-t**2))
erf2_ = f.integrate()
erf2 = erf2_ - erf2_(0)
from scipy.special import erf
randx = rng.random()
print(erf2(randx) - erf(randx))

-4.218847493575595e-15

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

In [54]:

def prod(f):
    return np.exp(np.log(f).sum())
def cumprod(f):
    return np.exp(np.log(f).integrate())

In [55]:

prod(np.exp(x))

Out[55]:

0.9999999999999991

In [56]:

prod(np.exp(np.exp(x)))

Out[56]:

10.489789833690239

In [57]:

cplot(cumprod(np.exp(x)))

Differentiation¶

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

In [58]:

f_sin5 = np.sin(5*x)
fd_sin5 = f_sin5.differentiate(4)
print(fd_sin5.norm()/f_sin5.norm())

625.0000000268575

In [59]:

f = np.sin(np.exp(x**2)).differentiate()
print(f(1))

-4.956699465919201

In [60]:

g = Chebfun.from_function(lambda x: 1/(2 + x**2))
h = g.differentiate()
print(h(1))

-0.22222222222884166

In [61]:

#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 [62]:

print((x - np.cos(x)).roots())

[0.73908513]

In [63]:

print((x - np.cos(4*x)).roots())

[ 0.31308831 -0.89882622 -0.53333306]

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

In [64]:

import scipy.special
def J_0(x):
    return scipy.special.jn(0,x)
#f = chebfun(lambda x: J_0(10*(1+x)))
f = Chebfun.from_function(J_0, domain=[0,20])
cplot(f)
roots = f.roots()
print('roots:', roots)
plt.plot(roots, np.zeros_like(roots), color='red', marker='.', linestyle='', markersize=10)
plt.grid()

roots: [ 2.40482556  5.52007811  8.65372791 18.07106397 11.79153444 14.93091771]

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 [65]:

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 [66]:

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 [67]:

f = np.sin(6*x) + np.sin(30*np.exp(x))
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():
    pychebfun.plot(f, with_interpolation_points=False)
    plt.plot(r, np.zeros_like(r), linestyle='', marker='o', color='gray', markersize=8)
    plt.plot(ma, f(ma), linestyle='', marker='o', color='green', markersize=8)
    plt.plot(mi, f(mi), linestyle='', marker='o', color='red', markersize=8)
    plt.grid()
plot_all()

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 [68]:

def s(x):
    return np.tan(x+1./4) + np.cos(10*x**2 + np.exp(np.exp(x)))
tancos = Chebfun.from_function(s)
cplot(tancos)
plt.title(str(tancos))

Out[68]:

Text(0.5, 1.0, '<Chebfun(73)>')

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

In [69]:

import scipy.integrate
scipy.integrate.quad(s, -1, 1, epsabs=1e-14)[0]

Out[69]:

0.29547767624377147

Using Chebyshev integration, we obtain exactly the same value:

In [70]:

tancos.sum()

Out[70]:

0.2954776762437714

ODE Solving¶

One solves the ODE

$\begin{align} u'(x) &= \mathrm{e}^{-2.75 x u(x)} \\ u(-1)&= 0 \end{align}$

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 fixed-point iteration:

In [71]:

uold = Chebfun(0.)
du = 1.
for i in range(100):
    integrand = np.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))))
pychebfun.plot(u)
plt.grid()
plt.title(str(u))

Converged in 32 steps
u(1) = 5.07818302385441

Out[71]:

Text(0.5, 1.0, '<Chebfun(149)>')

Boundary value problems¶

One solves the boundary problem

$u''(x) = \mathrm{e}^{4x}$ with boundary conditions

$u(-1) = 0 \qquad u(1) = 0$

In [72]:

f = np.exp(4*x)
u_ = f.integrate().integrate()
u = u_ - (u_(1)*(x+1) + u_(-1)*(1-x))/2.
cplot(u)
plt.grid()