Expansiones Bessel-Fourier.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bases de funciones\n",
    "\n",
    "\n",
    "Las funciones ${\\mathbb R}\\rightarrow {\\mathbb R}$ (${\\mathbb R}\\rightarrow {\\mathbb C}$) pueden pensarse como vectores en un espacio vectorial, definiendo\n",
    "$$\n",
    "(a f_1+ b f_2)(x) =a f_1(x)+b f_2(x)\n",
    "$$\n",
    "de manera que obviamente, $0(x) := 0 f(x)$ es también un elemento de ese espacio. Lo interesante de definir a las funciones como vectores es que ahora podemos construir bases de funciones $\\{\\varphi_i(x)\\, i=0, 1, \\ldots \\}$, de forma que, cualquier $f$ se exprese formalmente como\n",
    "$$\n",
    "f(x)=\\sum_i a_i \\varphi_i(x)\n",
    "$$\n",
    "con $a$ ciertos coeficientes.\n",
    "\n",
    "Propiedades como la continuidad, la diferenciabilidad, analiticidad o la integrabilidad de una función, o la condición de ser solución de una ecuación diferencial lineal, definen subespacios en ese espacio de funciones, que también pueden expresarse en términos de ciertas bases. Por ejemplo, las funciones analíticas en el intervalo $(-1,1)$ se expresan como series de potencias, y por lo tanto, los *monomios* $x^k$ representan una base para este espacio de funciones. Veremos entonces que esta construcción es muy útil entre otras cosas, para resolver ecuaciones diferenciales lineales.\n",
    "\n",
    "Como las bases de funciones en general tienen infinitos elementos, pero la vida de los físicos y matemáticos es finita, para que este tipo de expansiones sea útil para evaluarlas, es necesario en algún punto \"truncar\" la base. En un espacio finito, esto es equivalente a *proyectar* los vectores sobre un hiperplano.\n",
    "Con esta idea, dada una función expresada como\n",
    "$$\n",
    "f(x)=\\sum_{k=0}^\\infty a_k \\varphi_k(x)\n",
    "$$\n",
    "definimos la función truncada a los primeros $n+1$ elementos de la base como\n",
    "$$\n",
    "f_n(x)=\\sum_{k=0}^n\\varphi_k(x)\n",
    "$$\n",
    "La medida de cuan bien aproximamos función la función $f_n$ con la función proyectada, se puede definir en términos de la distancia entre el vector original y el proyectado. \n",
    "$$\n",
    "|f-f_n|\n",
    "$$\n",
    "\n",
    "Existen muchas formas de definir estas distancias, y su elección dependerá de la aplicación concreta. Por ejemplo, si queremos asegurarnos de que para cualquier valor de $x$ el error cometido esté acotado, podemos definir la distancia según \n",
    "$$\n",
    "|f-f_n|_{\\infty}=\\max_x |f(x)-f_n(x)|\n",
    "$$\n",
    "Si la expansión que define a $f(x)$ tiene sentido, $\\lim_{n\\rightarrow \\infty} |f-f_n|\\rightarrow 0$ y decimos que $f_n$ converge **uniformemente** a $f$.\n",
    "\n",
    "\n",
    "\n",
    "Otras alternativas, menos restrictivas y más prácticas de calcular vienen dadas por\n",
    "$$\n",
    "|f-f_n|_1 = \\int |f(x)-f_n(x)| dx\n",
    "$$\n",
    "y\n",
    "$$\n",
    "|f-f_n|_2 = \\sqrt{\\int |f(x)-f_n(x)|^2 dx}\n",
    "$$\n",
    "Estas definiciones acotan el \"error promedio\" al estimar $f$ por $f_n$. De estas dos definiciones, la segunda tiene una forma particular: corresponde a la generalización \"continua\" de la distancia euclídeana\n",
    "$$\n",
    "|\\vec{v}|=\\sqrt{\\sum_i |v_i|^2}\n",
    "$$\n",
    "Esta distancia puede definirse en términos de un *producto interno* o métrica:\n",
    "$$\n",
    "|\\vec{v}|=\\sqrt{\\vec{v}^*\\cdot \\vec{v}}\n",
    "$$\n",
    "con $\\vec{v}^*\\cdot \\vec{v}=\\sum_i |\\vec{v}_i|^2$.\n",
    "\n",
    "Definiendo $f^* \\cdot g \\equiv (f,g)=\\int f^*(x) g(x) dx$ dotamos al espacio vectorial de las funciones (más rigurosamente, al subespacio de las funciones *cuadrado integrables*) de un producto interno. El (sub)espacio de las funciones cuadrado integrable equipado con este producto es lo que se conoce como un *Espacio de Hilbert*.\n",
    "\n",
    "Disponer de un producto interno en un espacio vectorial es muy práctico, ya que nos permite definir **bases ortogonales**, y recuperar los coeficientes de la expansión en esa base para cualquier vector.\n",
    "Consideremos por ejemplo las funciones suaves en el intervalo $(-1,1)$, la base de monomios $x^k$, y la función\n",
    "$f(x)=exp(-1/x^2)$. Si se tratase de una función analítica, podríamos identificar los coeficientes de la expansión en esta base con los coeficientes de su serie de potencias:\n",
    "$$\n",
    "a_n =\\frac{1}{n!} \\frac{\\partial^n}{ \\partial x^n}f(x)\n",
    "$$\n",
    "Sin embargo, en este caso los coeficientes se anulan indenticamente: como la función no es analítica en el dominio, la serie de potencias falla en representar correctamente a las funciones. \n",
    "\n",
    "Para construir la expansión, conviene realizar un cambio de base, a una base \"ortogonalizada\" (vía Gramm-Schmidt):\n",
    "\n",
    "$$\n",
    "p_l(x)= x^l - \\sum_{m=0}^{l-1} \\frac{\\int_{-1}^{1} p_{m}(x') (x')^l dx'}{\\int_{-1}^{1} p_{m}(x')^2 dx'} p_{m}(x)\n",
    "$$\n",
    "Estos polinomios (proporcionales a los que más adelante llamaremos **polinomios de Legendre**) permiten dar una expresión simple a la expansión, en términos de integrales:\n",
    "\n",
    "$$\n",
    "f(x)=\\sum_{l=0}^{\\infty} a_l p_l(x)\n",
    "$$\n",
    "con \n",
    "$$\n",
    "a_l = \\frac{\\int_{-1}^{1} p_{m}(x') f(x)  dx'}{\\int_{-1}^{1} p_{m}(x')^2 dx'}\n",
    "$$\n",
    "\n",
    "\n",
    "De manera semejante, podemos definir otras bases ortogonales respecto a otros productos. Por ejemplo, si definimos\n",
    "$$\n",
    "(f,g)=\\int f^*(x)g(x) w(x)dx\n",
    "$$\n",
    "con $w(x)>0$, podemos definir bases que sean más \"sensibles\" en las regiones donde $w(x)$ tome mayores valores. \n",
    "Por ejemplo, si elegimos $w(x)=\\frac{1}{\\sqrt{1-x^2}}$ obtenemos los llamados **polinomios de Tchevicheff**, que satisfacen \n",
    "$$\n",
    "T_{l}(\\cos(u))=\\cos(l u)\n",
    "$$\n",
    "\n",
    "Series de Fourier\n",
    "-------------------------\n",
    "\n",
    "Otra base de funciones ortogonales muy útil es la de las funciones armónicas:\n",
    "$$\n",
    "\\varphi_{+, n}(x)=\\cos(n \\pi \\theta)\\;\\;\\varphi_{-, n}(x)=\\sin(n \\pi \\theta)\n",
    "$$\n",
    "que cumplen\n",
    "$$\n",
    "\\int_{-1}^{1} \\varphi_{\\pm,n}(x)\\varphi_{\\pm,m}(x) dx = \\left\\{\n",
    "\\begin{array}{lr}\n",
    "1 & m=n\\\\\n",
    "0 & m\\neq n\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "y $\\int_{-1}^{1} \\varphi_{\\pm,n}(x)\\varphi_{\\mp,m}(x) dx = 0$\n",
    "\n",
    "\n",
    "Mostremos que este conjunto es un conjunto \"completo\". Para eso, debemos probar que\n",
    "\n",
    "$$\n",
    "S(x,x')=\\frac{1}{2}+\\sum_{m=1}^\\infty \\sin(\\pi m x)\\sin(\\pi m x')+\n",
    "\\sum_{m=1}^\\infty \\cos(\\pi m x)\\cos(\\pi m x') =\\delta(x-x')\n",
    "$$\n",
    "en el intervalo. Para eso, primero vamos a reducir un poco las sumas, notando que\n",
    "\n",
    "$$\n",
    "\\sin(\\pi m x)\\sin(\\pi m x')+\\cos(\\pi m x)\\cos(\\pi m x')  =\\cos(\\pi m (x-x'))\n",
    "$$\n",
    "de manera que podemos juntar ambas sumatorias y escribir\n",
    "$$\n",
    "\\sum_{m=1}^\\infty \\sin(\\pi m x)\\sin(\\pi m x')+\n",
    "\\sum_{m=1}^\\infty \\cos(\\pi m x)\\cos(\\pi m x')=\\sum_{m=1}^\\infty \\cos(\\pi m (x-x'))= \\Re\\sum_{m=0}^\\infty e^{{\\bf i}\\pi m (x-x')}-1\n",
    "$$\n",
    "Para construir la suma de exponenciales, vamos a evaluar la suma de los primeros $N$ términos, y luego  tomar el límite $N\\rightarrow \\infty$:\n",
    "$$\n",
    "\\sum_{m=0}^N e^{{\\bf i}\\pi m (x-x')}=\\frac{e^{{\\bf i}\\pi (N+1)(x-x')}-1}{e^{{\\bf i}\\pi (x-x')}-1}=\n",
    "e^{{\\bf i}\\pi (N)(x-x')/2}\\frac{\\sin(\\pi (N+1)(x-x')/2)}{\\sin(\\pi (x-x')/2)}\n",
    "$$\n",
    "Tomamos ahora la parte real, y obtenemos\n",
    "$$\n",
    "\\frac{\\cos(\\pi N(x-x')/2)\\sin(\\pi (N+1)(x-x')/2)}{\\sin(\\pi (x-x')/2)}\n",
    "$$\n",
    "ó, reduciendo el numerador,\n",
    "$$\n",
    "\\frac{\\sin(\\pi (2 N+1)(x-x')/2 )}{2\\sin(\\pi (x-x')/2)}+\\frac{1}{2}\n",
    "$$\n",
    "Remplazando encontramos que\n",
    "$$\n",
    "S(x,x')=\\lim_{N\\rightarrow \\infty}S_N(x-x')\n",
    "$$\n",
    "con \n",
    "$$\n",
    "S_N(u)=\\frac{\\sin(\\pi (2 N+1)u/2 )}{2\\sin(\\pi u/2)}\n",
    "$$\n",
    "\n",
    "$S_N$ es una función que conocemos de estudiar patrones de interferencia en óptica. Tiene un máximo en $u=0$ de valor $2N+1$, y fuera del origen oscila con frecuencia $\\approx N/2$. Esta oscilación implica que la integral de esta función sobre un intervalo que no contenga al origen se anula como $1/N$.  Sin embargo, (por construcción) su integral sobre todo el intervalo vale 1, para cualquier valor de $N$. Por lo tanto,\n",
    "$$\n",
    "S(x-x')=\\delta(x-x')\n",
    "$$\n",
    "con lo que probamos la completitud.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Gráfica de $S_n(u)$ para diferentes valores de $N$](sn.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hide_input": false,
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "\n",
    "\n",
    "\n",
    "func0 = lambda x: 3*np.exp(-1./(x**2+1.e-20)) if abs(x)<1 else 0.\n",
    "func1 = lambda x: np.sqrt(1.-x**2) if abs(x)<1 else 0.\n",
    "func2 = lambda x: 1./np.sqrt(1.-x**2) if abs(x)<1 else 0.\n",
    "func3 = lambda x: 1. if x>0 else 0.\n",
    "func4 = lambda x: x\n",
    "func5 = lambda x: np.exp(-x) if abs(x)>-1 else 0.\n",
    "func6 = lambda x: np.exp(-8.*(x-.5)**2)\n",
    "\n",
    "\n",
    "funcs = {\"3 e^{-1/x^2}\": func0,\n",
    "         \"\\sqrt{1-x^2}\": func1,\n",
    "         \"1/\\sqrt{1-x^2}\": func2,\n",
    "         \"x\": func4,\n",
    "         \"H(x)\": func3,\n",
    "         \"H(x)exp(-x)\": func5,\n",
    "         \"exp(-8.*(x-.5)^2)\": func6\n",
    "        }\n",
    "\n",
    "\n",
    "############## Bases ##################################\n",
    "\n",
    "def generate_interpolation_legendre(n, f, paridad=None):\n",
    "    from scipy.special import legendre\n",
    "    from scipy.integrate import quad\n",
    "    legendres_p = [legendre(l) for l in range(n+1)]\n",
    "    coeffs = []\n",
    "    for p in legendres_p:\n",
    "        a = quad(lambda x: f(x)*p(x),-1.,1.)[0]/quad(lambda x: p(x)**2, -1., 1.)[0]\n",
    "        coeffs.append(a)\n",
    "\n",
    "    if paridad:\n",
    "        for i in range(n+1):\n",
    "            if paridad == \"par\" and i%2==1:\n",
    "                coeffs[i] = 0\n",
    "            elif paridad == \"impar\" and i%2==0:\n",
    "                coeffs[i] = 0\n",
    "\n",
    "    return sum(legendres_p[l]*coeffs[l] for l in range(n+1))\n",
    "\n",
    "\n",
    "\n",
    "def generate_interpolation_tchevicheff(n, f, paridad=None):\n",
    "    from numpy.polynomial import chebyshev\n",
    "    from scipy.integrate import quad\n",
    "    tchevicheff = lambda l, u: np.cos(l * np.arccos(u))\n",
    "    coeffs = []\n",
    "    for l in range(n+1):\n",
    "        a = quad(lambda x: f(x)*tchevicheff(l, x)/np.sqrt(1-x**2),-1.,1.)[0]/quad(lambda x: tchevicheff(l, x)**2/np.sqrt(1-x**2), -1., 1.)[0]\n",
    "        coeffs.append(a)\n",
    "\n",
    "    if paridad:\n",
    "        for i in range(n+1):\n",
    "            if paridad == \"par\" and i%2==1:\n",
    "                coeffs[i] = 0\n",
    "            elif paridad == \"impar\" and i%2==0:\n",
    "                coeffs[i] = 0\n",
    "\n",
    "    return lambda x: sum(tchevicheff(l,x)*coeffs[l] for l in range(n+1))\n",
    "\n",
    "##########################################################################################################\n",
    "\n",
    "def generate_interpolation_fourier(n, f, paridad=None):\n",
    "    from scipy.integrate import quad\n",
    "    coeffs_s = []\n",
    "    coeffs_c = []\n",
    "    kmax = int((n-1)/2)   \n",
    "    a = quad(lambda x: f(x),-1.,1.)[0]/quad(lambda x: 1., -1., 1.)[0]\n",
    "    coeffs_c.append(a)\n",
    "    for k in range(1, kmax):\n",
    "        a = quad(lambda x: f(x)*np.cos(np.pi*k*x),-1.,1.)[0]/quad(lambda x: np.cos(np.pi*k*x)**2, -1., 1.)[0]\n",
    "        coeffs_c.append(a)\n",
    "        b = quad(lambda x: f(x)*np.sin(np.pi*k*x),-1.,1.)[0]/quad(lambda x: np.sin(np.pi*k*x)**2, -1., 1.)[0]\n",
    "        coeffs_s.append(b)\n",
    "        \n",
    "    if paridad:\n",
    "        if paridad == \"par\":\n",
    "            return lambda x: sum([coeffs_c[k]*np.cos(np.pi*k*x) for k in range(kmax)])\n",
    "        if paridad == \"impar\":\n",
    "            return lambda x: sum([coeffs_s[k]*np.sin(np.pi*(k+1)*x) for k in range(kmax-1)])\n",
    "    return lambda x: (sum([coeffs_c[k]*np.cos(np.pi*k*x) for k in range(kmax)])+ \\\n",
    "                         sum([coeffs_s[k]*np.sin(np.pi*(k+1)*x) for k in range(kmax-1)]))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def generate_interpolation_bernstein(n, f, paridad=None):\n",
    "    from scipy.integrate import quad\n",
    "    from scipy.special import binom\n",
    "    fx = [f(x) for x in np.linspace(-1,1,n+1)]\n",
    "    if paridad:\n",
    "        if paridad==\"par\":\n",
    "            return lambda x: .5 * sum([(fx[k]+fx[n-k])*binom(n,k)*(x+1)**k * (x-1)**(n-k) for k in range(n+1)])\n",
    "        elif paridad==\"impar\":\n",
    "            return lambda x: .5 * sum([(fx[k]-fx[n-k])*binom(n,k)*(x+1)**k * (x-1)**(n-k) for k in range(n+1)])\n",
    "\n",
    "    return lambda x: sum([ fx[k]*binom(n,k)*(x+1)**k * (x-1)**(n-k) for k in range(n+1)])\n",
    "    \n",
    "\n",
    "bases = {\"Legendre\": generate_interpolation_legendre,\n",
    "         \"Tchevicheff\": generate_interpolation_tchevicheff,\n",
    "         \"Fourier\": generate_interpolation_fourier,\n",
    "        \"Bernstein\": generate_interpolation_fourier}\n",
    "\n",
    "###########################################################\n",
    "\n",
    "\n",
    "def plot_decomposition(fun, base, n, paridad=None):\n",
    "    from scipy.integrate import quad\n",
    "    basegen = bases[base]\n",
    "    f = funcs[fun]\n",
    "    xs = np.linspace(-1.1,1.1,100)\n",
    "    fapprox = basegen(n, f, paridad)\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.set_ylim(-1,3)\n",
    "    ax.plot(xs,[f(x) for x in xs],label=fun)\n",
    "    ax.plot(xs,[fapprox(x) for x in xs],label=\"$\"+fun +\"_{\"+str(n)+\"}$\" )\n",
    "    ax.plot(xs,[abs(f(x)-fapprox(x)) for x in xs],label=\"error\" )\n",
    "    l2err = quad(lambda x: (f(x)-fapprox(x))**2,-1.,1.)[0]\n",
    "    ax.text(-1., -.5, \"$|f-f_{n}|^2=\"+str(l2err)+\"$\")\n",
    "    ax.legend()\n",
    "    plt.show()\n",
    "\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hide_input": false,
    "init_cell": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "65029aa47997422b88038ea914c5d862",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='función:', options=('3 e^{-1/x^2}', '\\\\sqrt{1-x^2}', '1/\\\\sqrt{1-x…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ejemplo = widgets.Dropdown(options=[k for k in funcs.keys()],\n",
    "    description='función:',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "base = widgets.Dropdown(options=[k for k in bases.keys()],\n",
    "    description='Base de funciones:',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "paridad = widgets.Dropdown(options=['ambas', 'par', 'impar'],\n",
    "    description='Paridad:',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "\n",
    "nslider = widgets.IntSlider(min=0,max=20,step=1.,\n",
    "    description='n',\n",
    "    ensure_option=True,\n",
    "    disabled=False)\n",
    "\n",
    "\n",
    "interact = widgets.interact(plot_decomposition,fun=ejemplo, base=base,paridad=paridad, n=nslider)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Initialization Cell",
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}