{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chebfun in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a demo following [the original paper of Battles & Trefethen][1].\n",
    "\n",
    "[1]: http://people.maths.ox.ac.uk/trefethen/publication/PDF/2004_107.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_formats = {'svg', 'png'}\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It mostly suffices to use the class `Chebfun`, but we also import the whole `pychebfun` module for convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from pychebfun import Chebfun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pychebfun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "x = Chebfun.identity()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is only for this IPython notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a convenience plotting function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def cplot(fun):\n",
    "    pychebfun.plot(fun, with_interpolation_points=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chebfuns and barycentric interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creation of a chebfun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A typical way of initialising a Chebfun is to use the function `chebfun`. It will work for any function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cube = pychebfun.chebfun(lambda x:x**3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A more pythonic way to achieve the same result would be using `Chebfun.from_function`, it will also work for any function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cube = Chebfun.from_function(lambda x:x**3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another possibility is to use the variable `x` defined by `x = Chebfun.identity()`:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "x = Chebfun.identity()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cube = x**3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other examples could include:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print(np.cos(x))\n",
    "print(np.tan(x))\n",
    "print(np.exp(np.sin(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The _size_ of `f` is the number of interpolation point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cube.size()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Displaying a Chebfun gives the interpolation values at the Chebyshev interpolation points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print(repr(cube))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "Chebfun.from_function(lambda x: x**7 - 3*x + 5).size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_sin5 = Chebfun.from_function(lambda x: np.sin(5*np.pi*x))\n",
    "f_sin5.size()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convergence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For some functions, convergence is not achievable, and one must set a limit to the dichotomy algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_abs = Chebfun.from_function(abs, N=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Chebfun` objects behave as ordinary function, which can be evaluated at any point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cube(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cube(-.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible (and more efficient) to evaluate a chebfun at several point at once.\n",
    "For instance, we evaluate the Chebfun for $\\sin(5x)$ at several points at once:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_sin5(np.linspace(0, .2, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can for instance add two chebfuns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "chebsum = cube+f_sin5\n",
    "print(chebsum.size())\n",
    "cplot(chebsum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or multiply them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "chebmult = cube*f_sin5\n",
    "print(chebmult.size())\n",
    "cplot(chebmult)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_sin = Chebfun.from_function(lambda x:np.sin(x))\n",
    "print(f_sin.size())\n",
    "print((f_sin*f_sin).size())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to add and multiply chebfuns with scalars:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plt.subplot(121)\n",
    "cplot(3*f_sin)\n",
    "plt.subplot(122)\n",
    "cplot(3+f_sin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_cubic = Chebfun.from_function(lambda x: x**3 - x)\n",
    "cplot(f_cubic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cplot(f_sin5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2D Chebfun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the paper, they suggest to create two Chebfuns as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "np.sin(16*x), np.sin(18*x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is certainly possible, but we can't create a 2D Chebfun from them, at the moment.\n",
    "Instead we have to do it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def vector_func(x):\n",
    "    return np.vstack([np.sin(16*x), np.sin(18*x)]).T\n",
    "cheb_vec = Chebfun.from_function(vector_func)\n",
    "cplot(cheb_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elementary functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we explore the two ways to create a chebfun.\n",
    "First, the general way, which works for any smooth functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_expsin = Chebfun.from_function(lambda x: np.exp(np.sin(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, the “operator overloading” way.\n",
    "We take advantage of the fact that `exp` and `sin` are defined on chebfuns.\n",
    "Not all functions are, though.\n",
    "This is very similar to ufuncs in numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "g_expsin = np.exp(np.sin(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, both chebfuns are equivalent, as we demonstrate here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print(f_expsin.size())\n",
    "print(g_expsin.size())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plt.subplot(121)\n",
    "cplot(f_expsin)\n",
    "plt.subplot(122)\n",
    "cplot(g_expsin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "(f_expsin - g_expsin).norm()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approximation Theory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gibbs phenomenon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By limiting the accuracy of the chebfun, one observes the celebrated [Gibbs phenomenon](http://en.wikipedia.org/wiki/Gibbs_phenomenon)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_sign = Chebfun.from_function(lambda x: np.sign(x), N=25)\n",
    "cplot(f_sign)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pychebfun implements the method `compare` which plots a graph of the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def absfive(x):\n",
    "    return np.abs(x)**5\n",
    "#errors = array([(Chebfun(absfive, N) - exact).norm() for N in range(10, 20)])\n",
    "#loglog(errors)\n",
    "pychebfun.compare(Chebfun.from_function(absfive), absfive)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolation of random data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A chebfun passing through random values at 50 Chebyshev points.\n",
    "The method to initialise a Chebfun from data _at Chebyshev points_ is `from_data`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "rand_chebfun = Chebfun.from_data(rng.random(51))\n",
    "cplot(rand_chebfun)\n",
    "rand_chebfun.size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_walk = Chebfun.from_data(rng.normal(size=(100,2)))\n",
    "cplot(f_walk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extrapolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "This is an illustration of that phenomenon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(-1.4, 1.4, 100)\n",
    "error = np.log10(np.abs(f_sin5(xx)-np.sin(5*np.pi*xx))+1e-16)\n",
    "plt.subplot(211)\n",
    "plt.plot(xx, f_sin5(xx), marker='.')\n",
    "plt.ylabel('interpolant')\n",
    "plt.subplot(212)\n",
    "plt.plot(xx,error, marker='.')\n",
    "plt.ylabel('log_10(error)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chebyshev expansions and FFT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Chebyshev coefficients of a chebfun are quickly computed using a [fast Fourier transform](https://en.wikipedia.org/wiki/Fast_fourier_transform).\n",
    "The corresponding method is `coefficients`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "Chebfun.from_function(np.exp(x)).coefficients()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "Chebfun.from_coeff(np.array([3,2,1.]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to create the basis Chebyshev polynomial of order 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "T_10 = Chebfun.from_coeff(np.array([0.]*10+[1]))\n",
    "cplot(T_10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Incidentally, the same result is achieved using `Chebfun.basis`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cplot(Chebfun.basis(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an illustration of how fast the fast Fourier transform is, we create a chebfun with 100000 points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "r = rng.normal(size=100000)\n",
    "f_randn = Chebfun.from_coeff(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "np.sqrt(np.sum(np.square(r)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_randn.sum()\n",
    "#f_randn.norm() # doesn't work yet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "T_20 = Chebfun.basis(20)\n",
    "T_20.norm()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration and Differentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clenshaw Curtis Quadrature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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](http://en.wikipedia.org/wiki/Clenshaw%E2%80%93Curtis_quadrature))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "k_odd = 5\n",
    "k_even = 2\n",
    "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$.' % \n",
    "(Chebfun.basis(k_odd).sum(), \n",
    " k_odd,\n",
    "2./(1-k_even**2), \n",
    "Chebfun.basis(k_even).sum(), \n",
    "k_even))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing some integrals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f1 = np.sin(np.pi*x)**2\n",
    "print(f1.size(), f1.sum())\n",
    "HTML(r'$\\int \\sin(\\pi x)^2 dx = %s$' % f1.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f2 = Chebfun.from_function(lambda x: 1/(5+3*np.cos(np.pi*x)))\n",
    "print(f2.size(), f2.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f3 = Chebfun.from_function(lambda x: np.abs(x)**9 * np.log(np.abs(x)+1e-100))\n",
    "print(f3.size(), f3.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "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)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `dot` method computes the Hilbert scalar product, so `f.dot(g)` corresponds to\n",
    "\\\\[\\int_{-1}^1 f(x) g(x) \\mathrm{d} x\\\\]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "x.dot(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `integrate` method computes a primitive $F$ of a chebfun $f$, which is zero at zero i.e.,\n",
    "\\\\[F(x) = \\int_{0}^x f(y) \\mathrm{d}y\\\\]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = Chebfun.from_function(lambda t: 2/np.sqrt(np.pi) * np.exp(-t**2))\n",
    "erf2_ = f.integrate()\n",
    "erf2 = erf2_ - erf2_(0)\n",
    "from scipy.special import erf\n",
    "randx = rng.random()\n",
    "print(erf2(randx) - erf(randx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This allows to define continuous versions of the `prod` and `cumprod` functions for chebfuns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def prod(f):\n",
    "    return np.exp(np.log(f).sum())\n",
    "def cumprod(f):\n",
    "    return np.exp(np.log(f).integrate())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "prod(np.exp(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "prod(np.exp(np.exp(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cplot(cumprod(np.exp(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Differentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can also differentiate chebfuns to arbitrary orders with the method `differentiate`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_sin5 = np.sin(5*x)\n",
    "fd_sin5 = f_sin5.differentiate(4)\n",
    "print(fd_sin5.norm()/f_sin5.norm())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = np.sin(np.exp(x**2)).differentiate()\n",
    "print(f(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "g = Chebfun.from_function(lambda x: 1/(2 + x**2))\n",
    "h = g.differentiate()\n",
    "print(h(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "#f.differentiate().integrate() == f - f(-1)\n",
    "#f.integrate().differentiate() == f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operations based on rootfinding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As chebfun are polynomial it is possible to find all the roots in the interval $[-1,1]$ using the method `roots`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print((x - np.cos(x)).roots())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print((x - np.cos(4*x)).roots())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zeros of the Bessel function $J_0$ on the interval $[0,20]$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import scipy.special\n",
    "def J_0(x):\n",
    "    return scipy.special.jn(0,x)\n",
    "#f = chebfun(lambda x: J_0(10*(1+x)))\n",
    "f = Chebfun.from_function(J_0, domain=[0,20])\n",
    "cplot(f)\n",
    "roots = f.roots()\n",
    "print('roots:', roots)\n",
    "plt.plot(roots, np.zeros_like(roots), color='red', marker='.', linestyle='', markersize=10)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extrema"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The methods `min` and `max` are not implemented but it is not difficult to do so, by finding the roots of the derivative.\n",
    "This may come in a future version of `pychebfun`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = x - x**2\n",
    "#print('min', f.min())\n",
    "#print('max', f.max())\n",
    "#print('argmax', f.argmax())\n",
    "#print('norm inf', f.norm('inf'))\n",
    "#print('norm 1', f.norm(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Total variation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def total_variation(f):\n",
    "    return f.differentiate().norm(1)\n",
    "#total_variation(x)\n",
    "#total_variation(sin(5*pi*x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example of computation of extrema of a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = np.sin(6*x) + np.sin(30*np.exp(x))\n",
    "r = f.roots()\n",
    "fd = f.differentiate()\n",
    "fdd = fd.differentiate()\n",
    "e = fd.roots()\n",
    "ma = e[fdd(e) <= 0]\n",
    "mi = e[fdd(e) > 0]\n",
    "def plot_all():\n",
    "    pychebfun.plot(f, with_interpolation_points=False)\n",
    "    plt.plot(r, np.zeros_like(r), linestyle='', marker='o', color='gray', markersize=8)\n",
    "    plt.plot(ma, f(ma), linestyle='', marker='o', color='green', markersize=8)\n",
    "    plt.plot(mi, f(mi), linestyle='', marker='o', color='red', markersize=8)\n",
    "    plt.grid()\n",
    "plot_all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applications in numerical analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quadrature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We apply quadrature to the function\n",
    "$$ f(x) = \\tan(x+\\frac{1}{4}) + \\cos(10x^2 + \\exp(\\exp(x)))$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def s(x):\n",
    "    return np.tan(x+1./4) + np.cos(10*x**2 + np.exp(np.exp(x)))\n",
    "tancos = Chebfun.from_function(s)\n",
    "cplot(tancos)\n",
    "plt.title(str(tancos))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the built-in `quad` quadrature integrator of `scipy.integrate`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import scipy.integrate\n",
    "scipy.integrate.quad(s, -1, 1, epsabs=1e-14)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Chebyshev integration, we obtain exactly the same value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "tancos.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ODE Solving"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One solves the ODE\n",
    "\\begin{align}\n",
    "u'(x) &= \\mathrm{e}^{-2.75 x u(x)} \\\\\n",
    "u(-1)&= 0\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem is rewritten first as an integral problem, namely\n",
    "$$u(x) = T(u)$$\n",
    "with\n",
    "$$T(u) := \\int_{-1}^{x} \\mathrm{e}^{-2.75 yu(y)} dy$$\n",
    "This suggest using the following [fixed-point iteration](https://en.wikipedia.org/wiki/Fixed-point_iteration):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "uold = Chebfun(0.)\n",
    "du = 1.\n",
    "for i in range(100):\n",
    "    integrand = np.exp(-2.75*x*uold)\n",
    "    uint = integrand.integrate()\n",
    "    u = uint - uint(-1)\n",
    "    du = u-uold\n",
    "    uold = u\n",
    "    if du.norm() < 1e-13:\n",
    "        break\n",
    "else:\n",
    "    print('no convergence')\n",
    "print('Converged in {} steps'.format(i))\n",
    "print('u(1) = {:.14f}'.format(float(u(1))))\n",
    "pychebfun.plot(u)\n",
    "plt.grid()\n",
    "plt.title(str(u))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boundary value problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One solves the boundary problem \n",
    "$$u''(x) = \\mathrm{e}^{4x}$$\n",
    "with boundary conditions\n",
    "\\\\[u(-1) = 0 \\qquad u(1) = 0 \\\\]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = np.exp(4*x)\n",
    "u_ = f.integrate().integrate()\n",
    "u = u_ - (u_(1)*(x+1) + u_(-1)*(1-x))/2.\n",
    "cplot(u)\n",
    "plt.grid()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}