This is a demo following the original paper of Battles & Trefethen.
In principle, it suffices to use the class Chebfun
, but we'll import all the objects in chebfun
for convenience.
from pychebfun import *
x = Chebfun.identity()
This is only for this IPython notebook
from IPython.core.display import HTML
A typical way of initialising a Chebfun is to use the function chebfun
. It will work for any function.
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.
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.
cube = x**3
Other examples could include:
cos(x)
tan(x)
exp(sin(x))
<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]))>
The size of f
is the number of interpolation point.
cube.size()
4
Displaying a Chebfun gives the interpolation values at the Chebyshev interpolation points:
print repr(cube)
<Chebfun(array([ 1. , 0.125, -0.125, -1. ]))>
Chebfun.from_function(lambda x: x**7 - 3*x + 5).size()
8
f_sin5 = Chebfun.from_function(lambda x: sin(5*pi*x))
f_sin5.size()
42
For some functions, convergence is not achievable, and one must set a limit to the dichotomy algorithm.
f_abs = Chebfun.from_function(abs, N=1000)
Chebfun
objects behave as ordinary function, which can be evaluated at any point.
cube(5)
array(124.99999999999774)
cube(-.5)
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:
f_sin5(linspace(0, .2, 5))
array([ 4.33257766e-16, 7.07106781e-01, 1.00000000e+00, 7.07106781e-01, -2.68873793e-15])
One can for instance add two chebfuns:
chebsum = cube+f_sin5
print chebsum.size()
chebsum.plot()
42
Or multiply them:
chebmult = cube*f_sin5
print chebmult.size()
chebmult.plot()
43
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:
subplot(121)
(3*f_sin).plot()
subplot(122)
(3+f_sin).plot()
f_cubic = Chebfun.from_function(lambda x: x**3 - x)
f_cubic.plot()
f_sin5.plot()
In the paper, they suggest to create two Chebfuns as follows:
sin(16*x), sin(18*x)
(<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:
def vector_func(x):
return vstack([sin(16*x), sin(18*x)]).T
cheb_vec = Chebfun.from_function(vector_func)
cheb_vec.plot()
Again, we explore the two ways to create a chebfun. First, the general way, which works for any smooth functions:
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.
g_expsin = exp(sin(x))
Of course, both chebfuns are equivalent, as we demonstrate here.
print f_expsin.size()
print g_expsin.size()
20 20
subplot(121)
f_expsin.plot()
subplot(122)
g_expsin.plot()
(f_expsin - g_expsin).norm()
3.1401849173675503e-16
By limiting the accuracy of the chebfun, one observes the celebrated Gibbs phenomenon.
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.
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)
<matplotlib.axes.AxesSubplot at 0x106b88d50>
A chebfun passing through random values at 50 Chebyshev points.
The method to initialise a Chebfun from data at Chebyshev points is from_data
:
rand_chebfun = Chebfun.from_data(rand(51))
rand_chebfun.plot()
rand_chebfun.size()
51
f_walk = Chebfun.from_data(randn(100,2))
f_walk.plot()
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.
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)')
<matplotlib.text.Text at 0x106c08cd0>
The Chebyshev coefficients of a chebfun are quickly computed using a fast Fourier transform.
The corresponding method is chebyshev_coefficients
chebfun(exp(x)).chebyshev_coefficients()
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])
Chebfun.from_chebcoeff(array([3,2,1.]))
<Chebfun(array([ 6., 2., 2.]))>
One way to create the basis Chebyshev polynomial of order 10.
T_10 = Chebfun.from_chebcoeff(array([0.]*10+[1]))
T_10.plot()
Incidentally, the same result is achieved using Chebfun.basis
:
T_10_ = Chebfun.basis(10).plot()
As an illustration of how fast the fast Fourier transform is, we create a chebfun with 100000 points.
r = randn(100000)
f_randn = Chebfun.from_chebcoeff(r)
norm(r)
315.72498126138117
f_randn.sum()
#f_randn.norm() # doesn't work yet
2.7094293040817283
T_20 = Chebfun.basis(20)
T_20.norm()
1.4142135623730951
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)
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))
Computing some integrals.
f1 = sin(pi*x)**2
print f1.size(), f1.sum()
HTML('$\int \sin(\pi x)^2 dx = %s$' % f1.sum())
25 1.0
f2 = chebfun(lambda x: 1/(5+3*cos(pi*x)))
print f2.size(), f2.sum()
51 0.5
f3 = chebfun(lambda x: abs(x)**9 * log(abs(x)+1e-100))
print f3.size(), f3.sum()
91 -0.02
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)))
The dot
method computes the Hilbert scalar product, so f.dot(g)
corresponds to
$$\int_{-1}^1 f(x) g(x) \mathrm{d} x$$
x.dot(x)
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$$
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:
def prod(f):
return exp(log(f).sum())
def cumprod(f):
return exp(log(f).integrate())
prod(exp(x))
0.99999999999999911
prod(exp(exp(x)))
10.489789833690235
cumprod(exp(x)).plot()
One can also differentiate chebfuns to arbitrary orders with the method differentiate
f_sin5 = sin(5*x)
fd_sin5 = f_sin5.differentiate(4)
print fd_sin5.norm()/f_sin5.norm()
625.000000029
f = sin(exp(x**2)).differentiate()
print f(1)
-4.95669946592
g = Chebfun.from_function(lambda x: 1/(2 + x**2))
h = g.differentiate()
print h(1)
-0.222222222229
#f.differentiate().integrate() == f - f(-1)
#f.integrate().differentiate() == f
As chebfun are polynomial it is possible to find all the roots in the interval $[-1,1]$ using the method roots
:
print (x - cos(x)).roots()
[ 0.73908513]
print (x - cos(4*x)).roots()
[-0.89882622 -0.53333306 0.31308831]
Zeros of the Bessel function $J_0$ on the interval $[0,20]$:
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]
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
.
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
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.
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])
[-0.1, 0.1, -0.7, -0.6]
We apply quadrature to the function $$ f(x) = \tan(x+\frac{1}{4}) + \cos(10x^2 + \exp(\exp(x)))$$
def s(x):
return tan(x+1/4) + cos(10*x**2 + exp(exp(x)))
f = chebfun(s)
f.plot()
title(str(f))
<matplotlib.text.Text at 0x106f13b90>
Using the built-in quad
quadrature integrator of scipy.integrate
:
import scipy.integrate
scipy.integrate.quad(s, -1, 1, epsabs=1e-14)[0]
0.2954776762437714
Using Chebyshev integration, we obtain exactly the same value:
f.sum()
0.2954776762437713
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:
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
<matplotlib.text.Text at 0x106f0ac50>
One solves the boundary problem $$u''(x) = \mathrm{e}^{4x}$$ with boundary conditions $$u(-1) = 0 \qquad u(1) = 0 $$
f = exp(4*x)
u_ = f.integrate().integrate()
u = u_ - (u_(1)*(x+1) + u_(-1)*(1-x))/2
u.plot()
grid()