En este capítulo vamos a estudiar técnicas para aproximar una función cualquiera por una combinación lineal de funciones simples o conocidas como pueden ser las funciones polinómicas o las funciones trigonométricas.

La Teoría de la aproximación intenta resolver dos problemas:

• cuando la función a aproximar \(f\) viene dada de forma explícita, queremos hallar como hemos señalado antes otra función \(\hat{f}\), escrita en términos de funciones simples o más conocidas, de tal forma que \(\hat{f}\) esté lo más próxima a \(f\).
• cuando nos dan unos valores discretos \((x_k,y_k)\), \(k=0,\ldots,m\), queremos hallar la “mejor” función \(\hat{f}\) escrita también en términos de funciones simples o más conocidas, que aproxime la nube de puntos anterior.

En los problemas anteriores, no hemos concretado un punto importante: ¿qué queremos decir cuando decimos que \(\hat{f}\) esté lo más próxima a \(f\)?

Por ejemplo, en el tema dedicado a interpolación del curso Métodos numéricos con Python. Cálculo Numérico estudiamos cómo hallar el polinomio de interpolación de grado \(m-1\), \(p_{m-1}(x)\), dada una nube de puntos \((x_k,y_k),\ k=1,\ldots,m\), es decir, hallar un polinomio de grado menor o igual que \(m-1\) tal que \(p_{m-1}(x_k)=y_k\), \(k=1,\ldots,m\).

Cuando los puntos anteriores \((x_k,y_k)\) forman parte de la gráfica de una función \(f\), \(f(x_k)=y_k\), hallar el polinomio de interpolación sería una manera muy particular de aproximar dicha función \(f\).

Por tanto, en este caso, decir que \(\hat{f}=p_{m-1}\) esté cerca de \(f\) significa que la función interpoladora \(\hat{f}=p_{m-1}\) pase por los puntos de la nube: \(\hat{f}(x_k)=f(x_k)\).

Veremos otras formas de aproximar aparte de imponer que la función de aproximación \(\hat{f}\) pase por \(m\) puntos de la gráfica de la función a aproximar \(f\).

La aproximación por mínimos cuadrados consiste en, dado un conjunto de \(m\) puntos \((x_k,y_k)\), \(k=1,\ldots,m\), donde los \(x_k\) son diferentes \(2\) a \(2\), hallar una función \(\hat{f}(x)\) tal que la función error definida por: \[ E(\hat{f})=\sum_{k=1}^m (y_k-\hat{f}(x_k))^2, \] sea mínima.

La razón de que se llame aproximación por mínimos cuadrados es la expresión de la función error que consiste en la suma de los cuadrados de los errores cometidos en cada ordenada de cada punto \((x_k,y_k): (y_k-\hat{f}(x_k))^2\).

El siguiente paso es escoger el espacio de funciones donde elegiremos la función aproximación \(\hat{f}\) de la que hallaremos el mínimo.

Cuando el espacio de funciones considerado es el conjunto de funciones lineales, es decir, funciones del tipo \[ \hat{f}(x)=a_0+a_1 x, \] donde \(a_0\) y \(a_1\) son coeficientes libres que se tienen que calcular a partir de la nube de puntos \((x_k,y_k)\), \(k=1,\ldots,m\), la aproximación por mínimos cuadrados se denomina aproximación lineal por mínimos cuadrados.

Si la nube de puntos tiene una tendencia lineal como la del gráfico siguiente, hallar la aproximación lineal por mínimos cuadrados es hallar la mejor recta (en azul) tal que la función error antes definida es mínima, es decir, hallamos los coeficientes \(a_0\) y \(a_1\) tal que la suma de los errores (en rojo) al cuadrado entre las ordenadas \(y_k\) y los valores dados por la recta \(a_0+a_1 x_k\), \((y_k-a_0-a_1 x_k)^2\) es mínima.

Para hallar los coeficientes \(a_0\) y \(a_1\), tenemos que minimizar la función error que depende de dichos coeficientes: \[ E(a_0,a_1)=\sum_{k=1}^m (y_k-\hat{f}(x_k))^2 = \sum_{k=1}^m (y_k-a_0-a_1 x_k)^2. \] Para hallar los extremos de la función error \(E(a_0,a_1)\) de dos variables, tenemos que resolver el sistema siguiente: \[ \frac{\partial E(a_0,a_1)}{\partial a_0}=0,\quad \frac{\partial E(a_0,a_1)}{\partial a_1}=0. \]

Es decir: \[ \begin{align*} \frac{\partial E(a_0,a_1)}{\partial a_0} = & \frac{\partial}{\partial a_0}\sum_{k=1}^m (y_k-a_0-a_1 x_k)^2=-2\sum_{k=1}^m (y_k-a_0-a_1 x_k)\\ = & -2\left(\sum_{k=1}^m y_k-ma_0-a_1\sum_{k=1}^m x_k\right)=0,\\ \frac{\partial E(a_0,a_1)}{\partial a_1} = & \frac{\partial}{\partial a_1}\sum_{k=1}^m (y_k-a_0-a_1 x_k)^2=-2\sum_{k=1}^m (y_k-a_0-a_1 x_k)x_k\\ = & -2\left(\sum_{k=1}^m x_k y_k-a_0\sum_{k=1}^m x_k -a_1\sum_{k=1}^m x_k^2\right)=0. \end{align*} \]

Simplificando el sistema anterior, tenemos que los coeficientes \(a_0\) y \(a_1\) deben verificar el sistema de ecuaciones siguientes llamado sistema de ecuaciones normales: \[ \begin{align*} m a_0 + \left(\sum_{k=1}^m x_k\right) a_1 = & \sum_{k=1}^m y_k,\\ \left(\sum_{k=1}^m x_k\right) a_0+\left(\sum_{k=1}^m x_k^2\right) a_1 = & \sum_{k=1}^m x_k y_k. \end{align*} \]

Si llamamos \[ \Sigma_x =\sum_{k=1}^m x_k,\ \Sigma_y =\sum_{k=1}^m y_k,\ \Sigma_{x2} =\sum_{k=1}^m x_k^2,\ \Sigma_{xy} =\sum_{k=1}^m x_k y_k, \] el sistema anterior puede escribirse como: \[ \begin{align*} m a_0 + \Sigma_x a_1 = & \Sigma_y,\\ \Sigma_x a_0+\Sigma_{x2} a_1 = & \Sigma_{xy}, \end{align*} \] cuya solución puede obtenerse por el método de Cramer: \[ a_0 = \frac{\left|\begin{matrix}\Sigma_y & \Sigma_x\\ \Sigma_{xy}& \Sigma_{x2}\end{matrix}\right|}{\left|\begin{matrix}m & \Sigma_x\\ \Sigma_{x}& \Sigma_{x2}\end{matrix}\right|}=\frac{\Sigma_y\Sigma_{x2}-\Sigma_x\Sigma_{xy}}{m\Sigma_{x2}-\Sigma_x^2},\ a_1 = \frac{\left|\begin{matrix}m & \Sigma_y\\ \Sigma_{x}& \Sigma_{xy}\end{matrix}\right|}{\left|\begin{matrix}m & \Sigma_x\\ \Sigma_{x}& \Sigma_{x2}\end{matrix}\right|}=\frac{m\Sigma_{xy}-\Sigma_x\Sigma_{y}}{m\Sigma_{x2}-\Sigma_x^2}. \]

Usando las siguientes deficiones de estadística descriptiva: \[ \overline{x}=\frac{\Sigma_x}{m},\ \overline{y}=\frac{\Sigma_y}{m},\ s_{x}^2 =\frac{\Sigma_{x2}}{m}-\overline{x}^2,\ s_{xy}=\frac{\Sigma_{xy}}{m}-\overline{x}\overline{y}, \] las expresiones para \(a_0\) y \(a_1\) pueden escribirse como: \[ \begin{align*} a_1 =& \frac{m\Sigma_{xy}-\Sigma_x\Sigma_{y}}{m\Sigma_{x2}-\Sigma_x^2}=\frac{\frac{1}{m^2}\left(m\Sigma_{xy}-\Sigma_x\Sigma_{y}\right)}{\frac{1}{m^2}\left(m\Sigma_{x2}-\Sigma_x^2\right)}=\frac{\frac{\Sigma_{xy}}{m}-\overline{x}\overline{y}}{\frac{\Sigma_{x2}}{m}-\overline{x}^2}=\frac{s_{xy}}{s_x^2},\\ a_0 = & \frac{\Sigma_y\Sigma_{x2}-\Sigma_x\Sigma_{xy}}{m\Sigma_{x2}-\Sigma_x^2}=\frac{\frac{1}{m^2}\left(\Sigma_y\Sigma_{x2}-\Sigma_x\Sigma_{xy}\right)}{\frac{1}{m^2}\left(m\Sigma_{x2}-\Sigma_x^2\right)}=\frac{\overline{y}\frac{\Sigma_{x2}}{m}-\overline{x}\frac{\Sigma_{xy}}{m}}{s_x^2}\\ =& \frac{\overline{y}\left(\frac{\Sigma_{x2}}{m}-\overline{x}^2+\overline{x}^2\right)-\overline{x}\frac{\Sigma_{xy}}{m}}{s_x^2}=\frac{\overline{y}\left(s_x^2+\overline{x}^2\right)-\overline{x}\frac{\Sigma_{xy}}{m}}{s_x^2}\\ = & \overline{y}+\overline{x}\left(\frac{\overline{x}\overline{y}-\frac{\Sigma_{xy}}{m}}{s_x^2}\right)=\overline{y}-a_1\overline{x}. \end{align*} \]

Para calcular los coeficientes \(a_0\) y \(a_1\), hemos de realizar los pasos siguientes:

• calculamos \(\Sigma_x\), \(\Sigma_y\), \(\Sigma_{x2}\) y \(\Sigma_{xy}\),
• calculamos \(\overline{x},\overline{y},s_x^2\) y \(s_{xy}\) a partir de los valores anteriores,
• calculamos \(a_1=\frac{s_{xy}}{s_x^2}\) y \(a_0=\overline{y}-a_1\overline{x}\).

Observación.

Vemos que la recta por aproximación lineal por mínimos cuadrados pasa por el punto medio de la nube de puntos \((\overline{x},\overline{y})\) ya que como \(a_0=\overline{y}-a_1\overline{x}\) tenemos que \(a_0+a_1\overline{x}=\overline{y}\).

Observación.

Resolviendo las ecuaciones normales, sólo hemos comprobado que los valores \(a_0\) y \(a_1\) anteriores son un extremo de la función error. Deberíamos comprobar que los valores hallados son un mínimo.

Para ver que un valor candidato a extremo \((a_0,a_1)\) de una función de dos variables \(E(a_0,a_1)\) corresponde a un mínimo, hemos de ver que la matriz Hesiana evaluada en \((a_0,a_1)\) es definida positiva o que sus valores propios son positivos.

La matriz Hesiana es la siguiente: \[ H(a_0,a_1)=\begin{bmatrix}\frac{\partial^2 E}{\partial a_0^2} &\frac{\partial^2 E}{\partial a_0\partial a_1}\\ \frac{\partial^2 E}{\partial a_0\partial a_1}& \frac{\partial^2 E}{\partial a_1^2}\end{bmatrix} =\begin{bmatrix}2m&2\sum\limits_{k=1}^m x_k \\ 2\sum\limits_{k=1}^m x_k & 2\sum\limits_{k=1}^m x_k^2\end{bmatrix}. \]

Observación (continuación).

Antes de ver que la matriz anterior corresponde a un mínimo, necesitamos el lema siguiente:

Lema.

Dado un conjunto de datos \(x_1,\ldots,x_m\), donde hay como mínimo dos diferentes, su varianza poblacional \(s_x^2 = \frac{\sum\limits_{k=1}^m x_k^2}{m}-\overline{x}^2\) es estrictamente positiva, \(s_x^2> 0\).

Observación (continuación).

Para ver que la matriz \(H(a_0,a_1)\) es definida positiva podemos usar la regla de Sylvester, que dice que si los determinantes de orden \(1\) y \(2\) son positivos, la matriz \(H(a_0,a_1)\) es definida positiva: \[ \begin{align*} H(a_0,a_1)_1 =& 2m>0,\\ H(a_0,a_1)_2 = & 4m\sum_{k=1}^m x_k^2-4 \left(\sum_{k=1}^m x_k\right)^2 =\frac{4}{m^2}\left(\frac{\sum\limits_{k=1}^m x_k^2}{m}-\overline{x}^2\right)=\frac{4}{m^2}\cdot s_x^2 >0. \end{align*} \] En la última igualdad hemos aplicado el lema anterior.

Por tanto, concluimos que el valor \((a_0,a_1)\) hallado corresponde a un mínimo de la función error \(E(a_0,a_1)\).

• INPUT número de puntos de la nube \(m\), puntos de la nube \((x_k,y_k)\), \(k=1,\ldots,m\)
• Set \(Sigmax=0\).
• Set \(Sigmay=0\).
• Set \(Sigmax2=0\).
• Set \(Sigmaxy=0\).
• For k=1,...,m
  • Set \(Sigmax=Sigmax+x_k\).
  • Set \(Sigmay=Sigmay+y_k\).
  • Set \(Sigmax2=Sigmax2+x_k^2\).
  • Set \(Sigmaxy=Sigmaxy+x_k\cdot y_k\).

• Set \(\overline{x}=\frac{Sigmax}{m}\).
• Set \(\overline{y}=\frac{Sigmay}{m}\).
• Set \(s_x^2 = \frac{Sigmax2}{m}-\overline{x}^2\).
• Set \(s_{xy}= \frac{Sigmaxy}{m}-\overline{x}\cdot \overline{y}\).
• Set \(a_1=\frac{s_{xy}}{s_x^2}\).
• Set \(a_0=\overline{y}-a_1\cdot\overline{x}\).
• Print \(a_0,a_1\).

Vamos a generalizar el estudio anterior.

En lugar de considerar como espacio de funciones el conjunto de funciones lineales, vamos a considerar el conjunto de polinomios de grado máximo \(n\): \[ P_n(x)=a_n x^n+a_{n-1}x^{n-1}+\cdots +a_1 x+a_0, \] es decir, el problema que nos planteamos es el siguiente:

Problema de aproximación polinomial por mínimos cuadrados

Consideramos una nube de puntos \((x_k,y_k)\), \(k=1,\ldots,m\), donde los \(x_k\) son diferentes \(2\) a \(2\). Queremos hallar el polinomio \(P_n(x)\) de grado menor o igual que \(n\), con \(n< m-1\), tal que la función error siguiente sea mínima: \[ \begin{align*} E(a_0,a_1,\ldots,a_n)= & \sum_{k=1}^m (y_k-\hat{f}(x_k))^2 \\ =& \sum_{k=1}^m (y_k-(a_n x_k^n+a_{n-1}x_k^{n-1}+\cdots +a_1x_k+a_0))^2. \end{align*} \]

Observación.

La condición \(n<m-1\) se impone ya que sabemos que si \(n=m-1\) existe un polinomio denominado polinomio interpolador que cumple \(P_n(x_k)=y_k\).

Para este polinomio, tendremos que \(E(a_0,a_1,\ldots,a_n)=0\) ya que los sumandos que hay en el sumatorio de la función error son todos nulos.

En este caso, el problema quedaría resuelto usando dicho polinomio interpolador.

En resumen, intentamos solucionar el problema de aproximar de la mejor manera posible la nube de puntos por un polinomio cuyo grado no supere \(m-1\), siendo \(m\) el número de puntos de la nube.

La función error es una función de \(n+1\) variables: \(a_0,a_1,\ldots,a_n\).

Para hallar el mínimo de la misma hemos de resolver el sistema de ecuaciones siguiente: \[ \begin{align*} \frac{\partial E}{\partial a_0}(a_0,a_1,\ldots,a_n)= & 0,\\ \frac{\partial E}{\partial a_1}(a_0,a_1,\ldots,a_n)= & 0,\\ &\vdots\\ \frac{\partial E}{\partial a_n}(a_0,a_1,\ldots,a_n)= & 0. \end{align*} \]

El valor de \(\displaystyle\frac{\partial E}{\partial a_j}(a_0,a_1,\ldots,a_n)\), \(j=0,\ldots,n\) vale: \[ \begin{align*} & \frac{\partial E}{\partial a_j}(a_0,a_1,\ldots,a_n)=-2\sum_{k=1}^m (y_k-(a_n x_k^n+a_{n-1}x_k^{n-1}+\cdots +a_1 x_k+a_0))\cdot x_k^j\\ & =-2\left(\sum_{k=1}^m y_k\cdot x_k^j-a_n\sum_{k=1}^m x_k^{n+j}-a_{n-1}\sum_{k=1}^m x_k^{n-1+j}-\cdots - a_1\sum_{k=1}^m x_k^{1+j}-a_0\sum_{k=1}^m x_k^{j}\right), \end{align*} \] \(j=0,\ldots,n\).

La condición \(\displaystyle\frac{\partial E}{\partial a_j}(a_0,a_1,\ldots,a_n)=0\), \(j=0,\ldots,n\) puede escribirse de la forma siguiente: \[ \begin{align*} & \left(\sum_{k=1}^m x_k^{j}\right)a_0+\left(\sum_{k=1}^m x_k^{1+j}\right)a_1+\cdots +\left(\sum_{k=1}^m x_k^{n-1+j}\right)a_{n-1}+\left(\sum_{k=1}^m x_k^{n+j}\right)a_n\\ & \quad =\sum_{k=1}^m y_k\cdot x_k^j. \end{align*} \]

Fijémonos que el sistema anterior es lineal con \(n+1\) ecuaciones y \(n+1\) incógnitas \(a_0,a_1,\ldots,a_n\), que puede escribirse de la forma \(\mathbf{A}\mathbf{a}=\mathbf{b}\), donde \[ \begin{align*} \mathbf{A}= & (a_{ij})_{i,j=1,\ldots,n+1}=\left(\sum_{k=1}^m x_k^{i+j-2}\right)_{i,j=1,\ldots,n+1},\ \mathbf{a}=\begin{bmatrix}a_0\\ a_1\\\vdots\\ a_n\end{bmatrix},\\ \mathbf{b}= & \left(\sum_{k=1}^m y_k\cdot x_k^{i-1}\right)_{i=1,\ldots,n+1}. \end{align*} \] Para resolverlo, puede usarse un método directo como el método de Gauss con pivotaje o la descomposición \(LU\) (ver el curso Métodos numéricos con Python. Álgebra lineal numérica.)

Proposición.

El sistema anterior \(\mathbf{A}\mathbf{a}=\mathbf{b}\) tiene solución única.

• INPUT número de puntos de la nube \(m\), puntos de la nube \((x_k,y_k)\), \(k=1,\ldots,m\), \(n\) grado del polinomio aproximador
• Set \(\mathbf{A}=\mathbf{0}\). (Inicializamos la matriz \(\mathbf{A}\) del sistema \((n+1)\times (n+1)\) y el vector independiente \(\mathbf{b}\))
• Set \(\mathbf{b}=\mathbf{0}\).
• For k=1,...,m (Calculamos \(a_{ij}\), las componentes de la matriz del sistema \(\mathbf{A}\) y las componentes del vector \(\mathbf{b}\))
  • For i=1,...,n+1
    • Set \(b[i]=b[i]+y_k\cdot x_k^{i-1}\).
    • For j=1,...,n+1
      • Set \(A[i,j]=A[i,j]+x_k^{i+j-2}\)

• Solve \(\mathbf{A}\mathbf{a}=\mathbf{b}\) (Resolvemos el sistema para hallar los coeficientes del polinomio de aproximación \(\mathbf{a}\))
• Print \(a_0,a_1,\ldots,a_n\).

En ocasiones, la mejor forma de aproximar la nube de puntos no tiene por qué ser lineal e incluso polinómica.

Imaginemos una situación en que la relación entre los puntos sea exponencial con una función de tipo \(y=b\cdot\mathrm{e}^{ax}\) o potencial con una función del tipo \(y=b\cdot x^a\), donde las constantes \(a\) y \(b\) se tienen que hallar de forma que minimicen la función error \[ E(a,b)=\sum_{k=1}^m \left(y_k-b\mathrm{e}^{a x_k}\right)^2, \] en el caso de relación exponencial o, \[ E(a,b)=\sum_{k=1}^m \left(y_k-b\cdot x_k^a\right)^2, \] en el caso de relación potencial.

En el caso de relación exponencial de la nube de puntos, para hallar las “mejores” \(a\) y \(b\) tenemos que resolver el siguiente sistema de ecuaciones siguiendo el mismo procedimiento del método de mínimos cuadrados: \[ \begin{align*} & \frac{\partial E}{\partial a}(a,b)= -2\sum_{k=1}^m \left(y_k-b\mathrm{e}^{ax_k}\right)\cdot b\cdot x_k\cdot\mathrm{e}^{a x_k}=0,\ \stackrel{b\neq 0}{\Rightarrow} \\ & b\sum_{k=1}^m x_k\cdot \mathrm{e}^{2ax_k}=\sum_{k=1}^m x_k\cdot y_k\cdot\mathrm{e}^{ax_k}, \\ & \frac{\partial E}{\partial b}(a,b)= -2\sum_{k=1}^m \left(y_k-b\mathrm{e}^{ax_k}\right)\cdot \mathrm{e}^{a x_k}=0,\ \Rightarrow b\sum_{k=1}^m \mathrm{e}^{2ax_k}=\sum_{k=1}^m y_k\cdot \mathrm{e}^{ax_k}. \end{align*} \] Para hallar los coeficientes \(a\) y \(b\) que minimicen la función error, tendremos que resolver el sistema anterior.

Fijémonos que el sistema anterior es una sistema no lineal. Su resolución es en general compleja y no siempre tenemos asegurada la solución del mismo.

En el caso de relación potencial de la nube de puntos, para hallar las “mejores” \(a\) y \(b\) tenemos que resolver el siguiente sistema de ecuaciones siguiendo el mismo procedimiento del método de mínimos cuadrados: \[ \begin{align*} & \frac{\partial E}{\partial a}(a,b)= -2\sum_{k=1}^m \left(y_k-b\cdot x_k^a\right)\cdot b\cdot x_k^a\cdot\ln (x_k)=0,\ \stackrel{b\neq 0}{\Rightarrow} \\ & b\sum_{k=1}^m x_k^{2a}\cdot \ln (x_k)= \sum_{k=1}^m x_k^a\cdot y_k\cdot\ln(x_k),\\ & \frac{\partial E}{\partial b}(a,b)= -2\sum_{k=1}^m \left(y_k-b\cdot x_k^a\right)\cdot x_k^a=0,\ \Rightarrow \\ & b\sum_{k=1}^m x_k^{2a}= \sum_{k=1}^m x_k^a\cdot y_k. \end{align*} \] En este caso, para hallar los coeficientes \(a\) y \(b\) que minimicen la función error, tendremos que resolver el sistema anterior.

Fijémonos que el sistema anterior es también un sistema no lineal. Su resolución es en general compleja y no siempre tenemos asegurada la solución del mismo.

Una manera de solventar la dificultad anterior es realizar los cambios siguientes:

• Relación exponencial: en lugar de considerar la nube de puntos \((x_k,y_k),\ k=1,\ldots,m\), considerar la nube de puntos \((x_k,\ln (y_k))\) ya que si la relación buscada es del tipo \(y=b\cdot \mathrm{e}^{ax}\), tomando logaritmos en la expresión anterior, tenemos que \(\ln y=\ln b+a x\). Por tanto, buscamos una relación lineal en la nube de puntos \((x_k,\ln (y_k))\), \(k=1,\ldots,m\). Hallamos la aproximación lineal por mínimos cuadrados en la nube de puntos \((x_k,\ln (y_k))\), \(k=1,\ldots,m\) y hallamos \((a_0=\ln b,a_1=a)\). Entonces, la \(b\) buscada será \(b=\mathrm{e}^{a_0}\).

• Relación potencial: en lugar de considerar la nube de puntos \((x_k,y_k),\ k=1,\ldots,m\), considerar la nube de puntos \((\ln(x_k),\ln (y_k))\) ya que si la relación buscada es del tipo \(y=b\cdot x^{a}\), tomando logaritmos en la expresión anterior, tenemos que \(\ln y=\ln b+a\ln (x)\). Por tanto, buscamos una relación lineal en la nube de puntos \((\ln(x_k),\ln (y_k))\), \(k=1,\ldots,m\). Hallamos la aproximación lineal por mínimos cuadrados en la nube de puntos \((\ln(x_k),\ln (y_k))\), \(k=1,\ldots,m\) y hallamos \((a_0=\ln b,a_1=a)\). Entonces, la \(b\) buscada será \(b=\mathrm{e}^{a_0}\).

Como ya comentamos al principio del tema, vamos a intentar resolver el problema siguiente:

Problema de aproximación de funciones continuas.

Dada una función \(f\in {\cal C}[a,b]\) una función continua en un intervalo \([a,b]\), vamos a hallar una función \(\hat{f}\) en un determinado espacio de funciones tal que minimice el error cuadrático siguiente: \[ E(\hat{f})=\int_a^b (f(x)-\hat{f}(x))^2\, dx. \]

Entonces, en lugar de considerar una nube de puntos discreta \((x_k,y_k)\), \(k=1,\ldots,m\), consideramos todo un “conjunto continuo” de puntos \((x,f(x))\), \(x\in [a,b]\), y queremos hallar una aproximación \(\hat{f}(x)\), \(x\in [a,b]\) que minimice la función error anterior.

Como podemos observar, hemos “cambiado” los sumatorios por integrales ya que una integral puede definirse intuitivamente como un sumatorio continuo.

En el caso en que el espacio de funciones sea el conjunto de polinomios de grado menor o igual que un cierto \(n\), tendremos el problema de aproximación polinomial:

Problema de aproximación polinomial de funciones continuas.

Dada una función \(f\in {\cal C}[a,b]\) una función continua en un intervalo \([a,b]\), vamos a hallar un polinomio \(P(x)=a_n x^n+a_{n-1}x^{n-1}+\cdots +a_1 x+a_0\) de grado menor o igual que \(n\) tal que minimice el error cuadrático siguiente: \[ E(a_0,a_1,\ldots,a_n)=\int_a^b (f(x)-P(x))^2\, dx =\int_a^b \left(f(x)-\sum_{k=0}^n a_k x^k\right)^2\, dx. \]

La función error puede escribirse como: \[ E(a_0,a_1,\ldots,a_n)=\int_a^b f(x)^2\, dx-2\sum_{k=0}^n a_k\int_a^b x^k f(x)\, dx+\int_a^b\left(\sum_{k=0}^n a_k x^k\right)^2\, dx. \]

Como la función error \(E\) depende de \(n+1\) variables, el mínimo debe verificar el sistema de ecuaciones siguiente: \[ \begin{align*} \frac{\partial E}{\partial a_j}(a_0,a_1,\ldots,a_n)= 0,\ \Rightarrow & -2\int_a^b x^j f(x)\, dx+2\sum_{k=0}^n a_k\int_a^b x^{j+k}\, dx =0,\ \Rightarrow \\ \sum_{k=0}^n a_k\int_a^b x^{j+k}\, dx = & \int_a^b x^j f(x)\, dx,\ \Rightarrow \\ \sum_{k=0}^n a_k\left[\frac{x^{j+k+1}}{j+k+1}\right]_a^b = & \int_a^b x^j f(x)\, dx,\ \Rightarrow \\ \sum_{k=0}^n \left(\frac{b^{j+k+1}-a^{j+k+1}}{j+k+1}\right)\cdot a_k = & \int_a^b x^j f(x)\, dx, \end{align*} \] para \(j=0,\ldots,n\).

Obtenemos que el vector de coeficientes \(\mathbf{a}=(a_0,a_1,\ldots,a_n)\) verifica el sistema lineal siguiente denominado sistema de ecuaciones normales: \[ \mathbf{A}\mathbf{a}=\mathbf{b}, \] donde \(\mathbf{A}=(a_{ij})_{i,j=1,\ldots,n+1}\) es una matriz \((n+1)\times (n+1)\) de componentes \(a_{ij}=\frac{b^{i+j-1}-a^{i+j-1}}{i+j-1}\) y \(\mathbf{b}=(b_i)_{i=1,\ldots,n+1}\) es el término independiente de componentes \(b_i =\int_a^b x^{i-1} f(x)\, dx\):

\[ \mathbf{A}=\begin{bmatrix}b-a&\frac{b^2-a^2}{2}&\ldots&\frac{b^{n+1}-a^{n+1}}{n+1} \\\frac{b^2-a^2}{2}&\frac{b^3-a^3}{3}&\ldots&\frac{b^{n+2}-a^{n+2}}{n+2} \\\vdots&\vdots&\ddots&\vdots \\\frac{b^{n+1}-a^{n+1}}{n+1}&\frac{b^{n+2}-a^{n+2}}{n+2}&\ldots&\frac{b^{2n+1}-a^{2n+1}}{2n+1} \\\end{bmatrix},\ \mathbf{b}=\begin{bmatrix}\int_a^b f(x)\,dx\\ \int_a^b xf(x)\,dx\\\vdots\\\int_a^b x^{n}f(x)\,dx\end{bmatrix} \]

La matriz del sistema lineal anterior se denomina matriz de Hilbert.

Dicha matriz tiene un número de condición muy elevado, hecho que hace que aumenten los errores de redondeo de la solución. Ver el curso Métodos numéricos con Python - Álgebra lineal numérica para más detalles.

La tabla siguiente muestra el número de condición de la matriz de Hilbert para \(a=0\) y \(b=1\) usando la norma euclídea:

Vamos a resolver la aproximación polinomial de una forma computacional más eficiente.

Antes de hacerlo, necesitamos unas definiciones previas:

Definición de conjunto de funciones linealmente independientes.

Sean \(f_0,f_1\ldots,f_n\), \(n+1\) funciones definidas en el intervalo \([a,b]\). Diremos que el conjunto de funciones \(\{f_0,f_1,\ldots,f_n\}\) es linealmente independiente si la igualdad siguiente: \[ \lambda_0 f_0(x)+\lambda_1 f_1(x)+\cdots +\lambda_n f_n(x)=0, \] para cualquier \(x\in [a,b]\) implica que \(\lambda_0 =\lambda_1 =\cdots =\lambda_n =0\). En caso contrario, se dice que el conjunto de funciones es linealmente dependiente.

El resultado siguiente dice que cualquier conjunto de \(n+1\) polinomios de grado diferente es linealmente independiente:

Proposición: \(n+1\) polinomios de grados diferentes son linealmente independientes.

Sean \(P_0(x),P_1(x),\ldots,P_n(x)\), \(n+1\) polinomios donde el polinomio \(P_j(x)\) tiene grado \(j\), \(j=0,1,\ldots,n\). Entonces el conjunto de funciones polinómicas \(\{P_0(x),P_1(x),\ldots,P_n(x)\}\) es linealmente independiente.

Una consecuencia de la proposición anterior es la siguiente:

Proposición. Base de polinomios para el espacio vectorial de los polinomios de grado menor o igual que \(n\).

Sean \(P_0(x),P_1(x),\ldots,P_n(x)\), \(n+1\) polinomios donde el polinomio \(P_j(x)\) tiene grado \(j\), \(j=0,1,\ldots,n\). Entonces el conjunto de funciones polinómicas \(\{P_0(x),P_1(x),\ldots,P_n(x)\}\) forman una base para el espacio vectorial formado por los polinomios de grado menor o igual que \(n\).

Vamos a generalizar la aproximación polinomial de una función en términos de funciones peso.

Una función peso es una función que asocia más o menos importancia a los puntos del intervalo \([a,b]\):

Definición de función peso \(w\).

Diremos que una función \(w\) definida en el intervalo \((a,b)\) (¡ojo!, basta abierto) es una función peso si

• la función \(w(x)\) es integrable en el intervalo \([a,b]\). Por tanto, admitimos integración impropia convergente. Ver el curso Curso Completo de Cálculo en una variable: de cero a experto para más detalles,
• la función \(w(x)\) debe ser positiva, \(w(x)\geq 0\),
• la función \(w(x)\not\equiv 0\) no puede ser idénticamente nula en cualquier subintervalo \(I\subseteq [a,b]\).

Si tenemos una función peso \(w(x)\) en el intervalo \([a,b]\) y una base \(P_0(x),P_1(x),\ldots,P_n(x)\) de polinomios del espacio vectorial de polinomios de grado menor o igual que \(n\), podemos generalizar el problema de aproximación polinomial de la forma siguiente:

Problema de aproximación polinomial de funciones continuas.

Dada una función \(f\in {\cal C}[a,b]\) una función continua en un intervalo \([a,b]\) y una función peso \(w(x)\) vamos a hallar un polinomio \(P(x)=a_n P_n(x)+a_{n-1}P_{n-1}(x)+\cdots +a_1 P_1(x)+a_0 P_0(x)\) de grado menor o igual que \(n\) tal que minimice el error cuadrático siguiente: \[ \begin{align*} E(a_0,a_1,\ldots,a_n)= & \int_a^b w(x) (f(x)-P(x))^2\, dx \\ = & \int_a^b w(x)\left(f(x)-\sum_{k=0}^n a_k P_k(x)\right)^2\, dx. \end{align*} \]

Observación.

El problema de aproximación original se obtiene considerando la función peso \(w(x)\equiv 1\) y como base de polinomios del espacio vectorial de polinomios de grado menor o igual que \(n\), \(P_j(x)=x^j\), \(j=0,\ldots,n\).

Siguiendo un procedimiento similar, podemos escribir la función error de la forma siguiente: \[ \begin{align*} E(a_0,a_1,\ldots,a_n)= & \int_a^b w(x)f(x)^2\, dx-2\sum_{k=0}^n a_k\int_a^b w(x)P_k(x) f(x)\, dx\\ & \quad +\int_a^b w(x)\left(\sum_{k=0}^n a_k P_k(x)\right)^2\, dx. \end{align*} \]

A continuación, recordemos que para hallar los coeficientes \(a_j\), \(j=0,1,\ldots,n\), hemos de resolver el siguiente sistema de ecuaciones: Como la función error \(E\) depende de \(n+1\) variables, el mínimo debe verificar el sistema de ecuaciones siguiente: \[ \begin{align*} & \frac{\partial E}{\partial a_j}(a_0,a_1,\ldots,a_n)= 0,\ \Rightarrow \\ & -2\int_a^b w(x)P_j(x) f(x)\, dx+2\sum_{k=0}^n a_k\int_a^b w(x) P_j(x)P_k(x)\, dx = 0,\ \Rightarrow \\ & \sum_{k=0}^n a_k\int_a^b w(x) P_j(x)P_k(x)\, dx = \int_a^b w(x)P_j(x) f(x)\, dx, \end{align*} \] donde \(j=0,1,\ldots,n\).

Obtenemos que el vector de coeficientes \(\mathbf{a}=(a_0,a_1,\ldots,a_n)\) verifica el sistema lineal siguiente denominado sistema de ecuaciones normales: \[ \mathbf{A}\mathbf{a}=\mathbf{b}, \] donde \(\mathbf{A}=(a_{ij})_{i,j=1,\ldots,n+1}\) es una matriz \((n+1)\times (n+1)\) de componentes \(a_{ij}=\int_a^b w(x)P_{i-1}(x)P_{j-1}(x)\, dx\) y \(\mathbf{b}=(b_i)_{i=1,\ldots,n+1}\) es el término independiente de componentes \(b_i =\int_a^b w(x) P_{i-1}(x) f(x)\, dx\):

\[ \mathbf{A}=\begin{bmatrix}\int_a^b w(x)P_{0}(x)^2\, dx&\int_a^b w(x)P_{0}(x)P_1(x)\, dx&\ldots&\int_a^b w(x)P_{0}(x)P_n(x)\, dx \\\int_a^b w(x)P_{0}(x)P_1(x)\, dx&\int_a^b w(x)P_{1}(x)^2\, dx&\ldots&\int_a^b w(x)P_{1}(x)P_n(x)\, dx \\\vdots&\vdots&\ddots&\vdots \\\int_a^b w(x)P_{0}(x)P_n(x)\, dx&\int_a^b w(x)P_{1}(x)P_n(x)\, dx&\ldots&\int_a^b w(x)P_n(x)^2\, dx \\\end{bmatrix},\ \mathbf{b}=\begin{bmatrix}\int_a^b w(x)P_0(x)f(x)\,dx\\ \int_a^b w(x)P_1(x)f(x)\,dx\\\vdots\\\int_a^b w(x)P_n(x)f(x)\,dx\end{bmatrix}. \]

Aquí viene la idea clave: si pudiésemos elegir la base de polinomios \(P_0(x),P_1(x),\ldots,P_n(x)\) tal que: \[ \int_a^b w(x)P_i(x)P_j(x)=0, \] para \(i\neq j\) y \(i,j=0,\ldots,n+1\), la matriz \(\mathbf{A}\) del sistema lineal anterior sería diagonal y la resolución del sistema sería inmediata: \[ a_j = \frac{b_{j+1}}{a_{j+1,j+1}}=\frac{\int_a^b w(x)P_{j}(x)f(x)\,dx}{\int_a^b w(x)P_{j}(x)^2\, dx}, \] \(j=0,\ldots,n\).

La base de funciones que cumplan la condición anterior \(\int_a^b w(x)P_i(x)P_j(x)=0,\) para \(i\neq j\), \(i,j=0,\ldots,n\) se denominan base de funciones ortogonales:

Definición de funciones ortogonales.

Diremos que un conjunto de funciones \(\{f_0,f_1,\ldots,f_n\}\) es un conjunto ortogonal de funciones en el intervalo \([a,b]\) con respecto la función peso \(w(x)\) si: \[ \int_a^b w(x)f_i(x)f_j(x)\, dx=\begin{cases}0, & \mbox{si }i\neq j,\\ c_j>0, & \mbox{si } i=j.´\end{cases} \] Si además, \(c_0=\cdots = c_n=1\), el conjunto de funciones se denomina conjunto ortonormal de funciones.

El siguiente resultado dice cómo calcular la base de polinomios ortogonales dado un intervalo \([a,b]\) y una función peso \(w(x)\).

Teorema: construcción de polinomios ortogonales.

Sea \(\{P_0,P_1,\ldots,P_n\}\) el conjunto de polinomios definidos de la forma siguiente: \[ P_0(x)\equiv 1,\ P_1(x)=x-\alpha_1, \] para \(x\in [a,b]\) y \(\alpha_1 = \frac{\int_a^b x w(x)\, dx}{\int_a^b w(x)\, dx},\) y para \(k\geq 2\), \[ P_k(x)=(x-\alpha_k)P_{k-1}(x)-\beta_k P_{k-2}(x), \] para \(x\in [a,b]\) y \(\alpha_k=\frac{\int_a^b x w(x)P_{k-1}(x)^2\, dx}{\int_a^b w(x)P_{k-1}(x)^2\, dx},\) \(\beta_k=\frac{\int_a^b x w(x)P_{k-1}(x)P_{k-2}(x)\, dx}{\int_a^b w(x)P_{k-2}(x)^2\, dx}.\)

Entonces el conjunto anterior de polinomios es un conjunto ortogonal con respecto a la función peso \(w(x)\).

Observación.

Los polinomios ortogonales son únicos excepto por una constante.

Es decir, si \(\{P_0,P_1\ldots,P_n\}\) son \(n\) polinomios ortogonales definidos por el Teorema, también los son \(\{\gamma_0 P_0,\gamma_1 P_1\ldots,\gamma_n P_n\}\) donde \(\gamma_0,\gamma_1,\ldots,\gamma_n\) son \(n\) constantes arbitrarias..

El Teorema ha realizado una “normalización” eligiendo los polinomios \(P_i(x)\) mónicos, es decir, el coeficiente principal o el coeficiente de \(x^i\) siempre vale \(1\), \(i=0,1,\ldots,n\).

Observación (continuación).

Ahora bien, las soluciones de las ecuaciones normales \((a_i)_{i=0,\ldots,n}\) dependen de los polinomios ortogonales elegidos, es decir, dependen de las \(\gamma_i\), \(i=0,1,\ldots,n\) de la forma siguiente:

si en lugar de considerar \(P_i\), consideramos \(\gamma_i P_i\), \(i=0,1,\ldots,n\): \[ \tilde{a}_i = \frac{\int_a^b w(x)\gamma_i P_{i}(x)f(x)\,dx}{\int_a^b w(x)\gamma_i^2 P_{i}(x)^2\, dx} =\frac{1}{\gamma_i} \frac{\int_a^b w(x)P_{i}(x)f(x)\,dx}{\int_a^b w(x)P_{i}(x)^2\, dx}=\frac{1}{\gamma_i}a_i, \] siendo \(a_i\) las soluciones de las ecuaciones normales considerando los polinomios ortogonales \(P_i\) y \(\tilde{a}_i\), las soluciones de las ecuaciones normales considerando los polinomios ortogonales \(\gamma_i P_i\), \(i=0,1,\ldots,n\).

Usando la demostración del teorema anterior, tenemos el corolario siguiente:

Corolario sobre un conjunto de polinomios ortogonales.

Sea \(\{P_0,P_1,\ldots,P_n\}\) el conjunto de polinomios ortogonales definidos según el Teorema anterior. Entonces:

• El conjunto de polinomios \(\{P_0,P_1,\ldots,P_n\}\) es linealmente independiente.
• Cualquier polinomio \(P(x)\) de grado menor que \(n\) (grado menor o igual que \(n-1\)) es ortogonal al polinomio \(P_n(x)\) del conjunto anterior: \[ \int_a^b w(x)P_n(x)P(x)\, dx=0. \]

Una de las aproximaciones más usadas es la denominada aproximación minimax que consiste en resolver el problema siguiente:

Problema: aproximación minimax.

Sea \(f\in {\cal C}[a,b]\) una función continua en el intervalo \([a,b]\). Queremos hallar el polinomio \(\hat{P}_n(x)\) de grado menor o igual que \(n\) tal que minimice la norma infinito entre los polinomios de grado menor o igual que \(n\), \(P_n\) es decir: \[ \|f-\hat{P}_n\|_\infty =\max_{x\in [a,b]}|f(x)-\hat{P}_n(x)|=\min_{P_n}\|f-P_n\|_\infty =\min_{P_n}\max_{x\in [a,b]}|f(x)-P_n(x)|. \] Dicho en otras palabras, queremos hallar el polinomio \(\hat{P}_n\) de grado menor o igual que \(n\) que minimice el valor máximo del error \(|f(x)-\hat{P}_n(x)|\) dentro del intervalo \([a,b]\) de entre todos los polinomios de grado menor o igual que \(n\).

La razón que se llame aproximación minimax está en la propia definición del problema: \(\min_{P_n}\max_{x\in [a,b]}|f(x)-P_n(x)|\).

El siguiente resultado nos caracteriza el polinomio buscado \(\hat{P}_n\):

Teorema: caracterización de la aproximación minimax.

El polinomio \(\hat{P}_n\) que minimiza la norma infinito \(\|f-P_n\|_\infty\) de entre todos los polinomios de grado menor o igual que \(n\) verifica lo siguiente:

la función error \(e_n(x)=f(x)-\hat{f}(x)\) tiene como mínimo \(n+2\) abscisas \(\xi_i\), \(i=0,1,\ldots,n+1\) ordenadas en el intervalo \([a,b]\) de la forma \(a\leq\xi_0<\xi_1<\cdots <\xi_{n+1}\leq b\) verificando la propiedad de equioscilación: \[ |e_n(\xi_j)|=\|e_n\|_\infty =\max_{x\in [a,b]}|e_n(x)|=\max_{x\in [a,b]}|f(x)-\hat{P}_n(x)|, \] para \(j=0,1,\ldots,n+1\), con alternancia de signos: \[ e_n(\xi_0)=-e_n(\xi_1)=\cdots = (-1)^{n+1}e_n(\xi_{n+1}). \]

Observación.

El Teorema anterior nos dice que si hallamos un polinomio de grado menor o igual que \(n\), \(\hat{P}_n\), tal que existen \(n+1\) valores reales en el intervalo \([a,b]\), \(a\leq\xi_0<\xi_1<\cdots <\xi_{n+1}\leq b\), cumpliendo \(|e_n(\xi_j)|=\|e_n\|_\infty\), \(j=0,\ldots,n+1\) y \(e_n(\xi_j)=-e_n(\xi_{j-1})\), \(j=1,\ldots,n+1\), entonces dado un polinomio cualquiera \(P_n\) de grado menor o igual que \(n\), se cumple: \[ \|f-\hat{P}_n\|_\infty\leq \|f-P_n\|_\infty. \]

Proposición.

La aproximación minimax \(\hat{P}_n\) es única.

Los polinomios de Chebyschev \((T_k(x))_{k\geq 0}\) son los polinomios ortogonales correspondientes al intervalo \([-1,1]\) y a la función peso \(w(x)=\frac{1}{\sqrt{1-x^2}}\) y como veremos a continuación tienen un papel fundamental en la aproximación minimax.

En un ejemplo anterior, demostramos que dicha función peso cumplía las propiedades correspondientes.

Aunque los polinomios de Chebyschev se pueden definir a partir de la recurrencia siguiente tal como hemos visto: \[ \begin{align*} T_0(x)= & 1,\\ T_1(x)= & x-\alpha_1,\ \alpha_1 = \frac{\int_{-1}^1 \frac{x}{\sqrt{1-x^2}}\, dx}{\int_{-1}^1 \frac{1}{\sqrt{1-x^2}}\, dx}=0,\ \Rightarrow T_1(x)=x, \end{align*} \]

\[ T_k(x) = (x-\alpha_k)T_{k-1}(x)-\beta_k T_{k-2}(x), \] donde \[ \alpha_k =\frac{\int_{-1}^1 \frac{x}{\sqrt{1-x^2}}T_{k-1}(x)^2\, dx}{\int_{-1}^1 \frac{1}{\sqrt{1-x^2}} T_{k-1}(x)^2\, dx},\quad \beta_k =\frac{\int_{-1}^1 \frac{x}{\sqrt{1-x^2}}T_{k-1}(x) T_{k-2}(x)\, dx}{\int_{-1}^1 \frac{1}{\sqrt{1-x^2}} T_{k-2}(x)^2\, dx}, \] vamos a definirlos de la siguiente manera mucho más simple y más fácil de trabajar: \[ T_n(x)=\cos(n\arccos(x)),\ n\geq 0. \]

Teorema. Expresión de los polinomios de Chebyschev.

Los polinomios de Chebyschev cumplen la expresión siguiente: \[ T_n(x)=\cos(n\arccos(x)),\ n\geq 0, \] junto con la recurrencia: \[ T_{n+1}(x)=2xT_n(x)-T_{n-1}(x). \]

Observación.

Queremos realizar una aclaración. Cuando decimos que los polinomios de Chebyschev se pueden escribir como \(T_n(x)=\cos(n\arccos(x))\) estamos diciendo que salvo una constante, se obtienen los polinomios de Chebyschev dados por la recurrencia general en los polinomios ortogonales.

Recordemos que los polinomios ortogonales vistos por la recurrencia general son mónicos.

Sin embargo, veremos que los polinomios dados por la expresión \(T_n(x)=\cos(n\arccos(x))\) no lo son. De hecho, veremos que el coeficiente principal o el coeficiente de \(x^n\) del polinomio \(T_n(x)\) será \(2^{n-1}\).

Pasemos, pues, a demostrar que efectivamente los polinomios de Chebyschev salvo una constante se pueden definir de la forma anterior.

Observación.

Para \(n=m\), usando el mismo procedimiento anterior, tenemos que: \[ \begin{align*} \int_{-1}^1 \frac{T_n(x)^2}{\sqrt{1-x^2}}\, dx = & \int_{0}^{\pi}\cos(n\theta)^2\,d\theta =\int_{0}^{\pi}\left(\frac{1+\cos(2n\theta)}{2}\right)\,d\theta \\ = & \frac{1}{2}\left(\int_0^{\pi}1\,d\theta +\int_0^{\pi}\cos(2n\theta)\,d\theta\right)=\frac{1}{2}\left(\pi +\left[\frac{\sin(2n\theta)}{2n}\right]_0^\pi\right)=\frac{\pi}{2}, \end{align*} \] ya que \(\sin(2n\pi)=\sin(2n0)=0\), para todo \(n\geq 1\).

Para el caso \(n=0\), obtenemos: \[ \int_{-1}^1 \frac{T_0(x)^2}{\sqrt{1-x^2}}\, dx =\int_{-1}^1 \frac{1}{\sqrt{1-x^2}}\, dx = \int_{0}^{\pi}1\,d\theta =\pi. \]

Vamos a calcular los primeros polinomios de Chebyschev: \[ \begin{align*} T_0 (x)= & 1,\\ T_1 (x) = & x,\\ T_2 (x) = & 2xT_1(x)-T_0(x)=2x\cdot x-1=2x^2-1,\\ T_3 (x) = & 2xT_2(x)-T_1(x)=2x(2x^2-1)-x=4x^3-3x,\\ T_4 (x) = & 2x T_3(x)-T_2(x)=2x(4x^3-3x)-(2x^2-1)=8x^4-8x^2+1,\\ T_5 (x) = & 2x T_4(x)-T_3(x)=2x(8x^4-8x^2+1)-(4x^3-3x)=16 x^5-20x^3+5x, \end{align*} \] y así sucesivamente.

La gráfica siguiente muestra los polinomios de Chebyschev \(T_i(x)\), \(i=1,2,3,4,5\), donde \(T_1(x)\) está en negro, \(T_2(x)\), en morado, \(T_3(x)\), en verde, \(T_4(x)\), en azul y \(T_5(x)\) en rojo.

Veamos cuál es la relación de los polinomios de Chebyschev con la aproximación minimax. En primer lugar, necesitamos los resultados siguientes:

Proposición: ceros y extremos de los polinomios de Chebyschev.

El polinomio de Chebyschev \(T_n(x)\) de grado \(n\) tiene los \(n\) ceros siguientes en el intervalo \([-1,1]\): \[x_k = \cos\left(\frac{(2k+1)\pi}{2n}\right),\ \mbox{para }k=0,\ldots,n-1,\] y tiene los siguientes \(n+1\) extremos relativos en el intervalo \([-1,1]\): \[\hat{x}_k=\cos\left(\frac{k\pi}{n}\right),\ \mbox{para }k=0,1,\ldots,n,\] con valores extremos \(T_n(\hat{x}_k)=(-1)^k\).

Proposición. Coeficiente del término de mayor grado del polinomio de Chebyschev.

El coeficiente del término de \(x^n\) del polinomio de Chebyschev de grado \(n\), \(T_n(x)\) vale \(2^{n-1}\) para \(n\geq 1\).

Definimos los polinomios mónicos de Chebyschev \(\tilde{T}_n\) como los polinomios de Chebyschev divididos por su término principal de \(x^n\), \(2^{n-1}\):

Definición de los polinomios mónicos de Chebyschev.

El polinomio mónimo de Chebyschev \(\tilde{T}_n\) se define como: \(\tilde{T}_n=\frac{T_n(x)}{2^{n-1}}\).

Se llaman mónicos porque el coeficiente principal del polinomio mónico de Chebyschev de grado \(n\), \(\tilde{T}_n\) vale \(1\).

De esta forma podemos escribir: \(\tilde{T}_n(x) = x^n+\cdots\)

Usando que los polinomios de Chebyschev verifican la recurrencia \(T_{n}(x)=2xT_{n-1}(x)-T_{n-2}(x),\) podemos deducir la recurrencia que verifican los polinomios mónicos de Chebyschev: \[ \begin{align*} T_{n}(x)=& 2xT_{n-1}(x)-T_{n-2}(x),\ \Rightarrow \frac{T_n(x)}{2^{n-1}}=\frac{2xT_{n-1}(x)}{2^{n-1}}-\frac{T_{n-2}(x)}{2^{n-1}},\\ \frac{T_n(x)}{2^{n-1}}=&\frac{2x}{2}\cdot \frac{T_{n-1}(x)}{2^{n-2}}-\frac{1}{2^2}\cdot \frac{T_{n-2}(x)}{2^{n-3}},\ \Rightarrow \\ \tilde{T}_n = & x\cdot\tilde{T}_{n-1}-\frac{1}{4}\tilde{T}_{n-2}. \end{align*} \]

Los primeros polinomios mónicos de Chebyschev son los siguientes:

\[ \begin{align*} \tilde{T}_0(x) = & 1,\\ \tilde{T}_1(x) = & x,\\ \tilde{T}_2 (x)= & \frac{T_2(x)}{2} =\frac{2x^2-1}{2}=x^2-\frac{1}{2},\\ \tilde{T}_3 (x)= & x\tilde{T}_2(x) -\frac{1}{4}\tilde{T}_1(x) =x\left(x^2-\frac{1}{2}\right)-\frac{1}{4}x =x^3-\frac{3}{4}x,\\ \tilde{T}_4 (x)= & x\tilde{T}_3(x) -\frac{1}{4}\tilde{T}_2(x) =x\left(x^3-\frac{3}{4}x\right)-\frac{1}{4}\left(x^2-\frac{1}{2}\right)=x^4-x^2+\frac{1}{8}, \end{align*} \] y así sucesivamente.

Los primeros polinómicos mónicos de Chebyschev se muestran en la tabla siguiente:

El polinomio mónico de Chebyschev de grado \(n\), \(\tilde{T}_n\), minimiza la norma infinito en el intervalo \([-1,1]\) entre todos los polinomios mónicos de grado \(n\):

Teorema: minimización de la norma infinito.

Sea \(P_n(x)\) un polinomio mónico cualquiera de grado \(n\). Entonces, \[ \frac{1}{2^{n-1}}=\|\tilde{T}_n\|_\infty =\max_{x\in [-1,1]}|\tilde{T}_n(x)|\leq \|P_n\|_\infty =\max_{x\in [-1,1]}|P_n(x)|. \] Además la igualdad ocurre sólo si \(P_n\equiv\tilde{T}_n\).

En la demostración del Teorema anterior, hemos demostrado el corolario siguiente:

Corolario: aproximación minimax de la función \(x^n\) en el intervalo \([-1,1]\).

La aproximación minimax de la función \(f(x)=x^n\) en el intervalo \([-1,1]\) entre todos los polinomios de grado menor o igual que \(n-1\) es la función \(\hat{P}_{n-1}(x)=x^n-\tilde{T}_n(x)\), es decir: \[ \min_{P_{n-1}}\|x^n-P_{n-1}(x)\|_\infty =\|x^n-(x^n-\tilde{T}_n(x))\|_\infty. \]

En el tema dedicado a interpolación del curso Métodos numéricos con Python. Cálculo Numérico estudiamos cómo hallar el polinomio de interpolación de grado \(m\), \(P_m(x)\), dada una nube de puntos \((x_k,y_k),\ k=1,\ldots,m\), es decir, hallar un polinomio de grado menor o igual que \(m-1\) tal que \(P_m(x_k)=y_k\), \(k=1,\ldots,m\).

Además, en el caso en que la nube de puntos provenga de la gráfica de una función \(f(x)\), es decir \(y_k=f(x_k)\), donde dicha función \(f\in {\cal C}^{m+1}\) es de clase \(m+1\), el error de interpolación se puede escribir de la forma siguiente: \[ f(x)-P_m(x)=\frac{f^{(m+1)}(\xi(x))}{(m+1)!}(x-x_0)(x-x_1)\cdots (x-x_m). \]

Notemos que el error depende del polinomio mónico \((x-x_0)(x-x_1)\cdots (x-x_m)\) de grado \(m+1\).

Si suponemos que los nodos de interpolación están en el intervalo \([-1,1]\), nos podríamos plantear cómo elegir los nodos \(x_k,\ k=0,1,\ldots,m\) en dicho intervalo de tal forma que el máximo del polinomio anterior \(|(x-x_0)(x-x_1)\cdots (x-x_m)|\) sea mínimo. Es decir, nos planteamos el problema siguiente: \[ \begin{align*} & \min_{x_0,x_1,\ldots,x_m}\max_{x\in [-1,1]}|(x-x_0)(x-x_1)\cdots (x-x_m)|\\ & \quad = \min_{x_0,x_1,\ldots,x_m}\|(x-x_0)(x-x_1)\cdots (x-x_m)\|_\infty. \end{align*} \] El problema anterior es equivalente a hallar el polinomio mónico de grado \(m+1\) que minimice la norma infinito en el intervalo \([-1,1]\).

Dicho problema ya está resuelto eligiendo el polinomio mónico de Chebyschev de grado \(m+1\), \(\tilde{T}_{m+1}\).

Por tanto, los nodos \(x_k,\ k=0,1,\ldots,m+1\) se tienen que elegir de tal forma que se cumpla: \[ \tilde{T}_{m+1}(x)=(x-x_0)(x-x_1)\cdots (x-x_m), \] es decir, los nodos \(x_k,\ k=0,1,\ldots,m\) deben ser los ceros del polinomio de Chebyschev \(T_{m+1}(x)\): \[ x_k = \cos\left(\frac{(2k+1)\pi}{2(m+1)}\right),\ k=0,1,\ldots,m. \]

Hemos demostrado el corolario siguiente:

Corolario.

Sea \(f\in {\cal C}^{m+1}\) una función de clase \(m+1\) y sea \(P_m(x)\) el polinomio de interpolación en los ceros \(x_k=\cos\left(\frac{(2k+1)\pi}{2(m+1)}\right),\ k=0,1,\ldots,m\) del polinomio de Chebyschev de grado \(m+1\). Entonces, \[ \max_{x\in [-1,1]}|(x-x_0)(x-x_1)\ldots (x-x_m)|=\|\tilde{T}_{m+1}\|_\infty = \frac{1}{2^m}. \]

Corolario. (continuación)

Además el error de interpolación \(f(x)-P_m(x)\) en cualquier \(x\in [-1,1]\) se puede acotar por: \[ \begin{align*} \max_{x\in [-1,1]}|f(x)-P_m(x)|= & \max_{x\in [-1,1]}\left|\frac{f^{(m+1)}(\xi(x))}{(m+1)!}(x-x_0)(x-x_1)\cdots (x-x_m)\right|\\\leq &\frac{1}{2^m\cdot (m+1)!}\max_{x\in [-1,1]}\left|f^{(m+1)}(\xi(x))\right|\leq \frac{\|f^{(m+1)}\|_\infty}{2^m\cdot (m+1)!} \end{align*} \]

Observación.

El fenómeno de Runge que vimos en el curso Métodos Numéricos con Python. Cálculo Numérico que ocurría cuando interpolamos la función \(f(x)=\frac{1}{1+5x^2}\) en el intervalo \([-1,1]\) desaparece cuando consideramos como nodos de interpolación los ceros del polinomio de Chebyschev correspondiente.

Supongamos que tenemos una función \(f\in {\cal C}^{m+1}\) de clase \(m+1\) en un intervalo cualquiera \([a,b]\) y consideramos \(t_0,t_1,\ldots,t_m\), \(m+1\) puntos en dicho intervalo.

Sea \(P_m(t)\) el polinomio de interpolación de grado menor o igual que \(m\) en los puntos anteriores, es decir \(P_m(t_k)=f(t_k),\ k=0,1,\ldots,m\)

Sabemos que el error de interpolación (ver tema de interpolación del curso Métodos numéricos con Python. Cálculo Numérico) tiene la expresión siguiente: \[ f(t)-P_m(t)=\frac{f^{m+1}(\xi(t))}{(m+1)!}(t-t_0)(t-t_1)\cdots (t-t_m). \]

A continuación, nos planteamos la misma cuestión que hacíamos en el intervalo \([-1,1]\): ¿existe alguna manera de elegir los nodos \(t_k,\ k=0,1,\ldots,m\) que minimice la norma infinito de la parte del error de interpolación que depende de dichos nodos, \(\|(t-t_0)(t-t_1)\cdots (t-t_m)\|_\infty =\max_{t\in [a,b]}|(t-t_0)(t-t_1)\cdots (t-t_m)|\)?

Para responder a dicha pregunta, basta considerar los ceros del polinomio de Chebyschev de grado \(m\), \(T_m(x)\) pero considerados en el intervalo \([a,b]\).

Es decir, hay que considerar la siguiente aplicación que pasa del intervalo \([-1,1]\) al intervalo \([a,b]\), \(t=\frac{1}{2}((b-a)x+a+b)\) y transformar los ceros del polinomio de Chebyschev de grado \(m\) usando dicha transformación.

Como los ceros del polinomio de Chebyschev de grado \(m\) en el intervalo \([-1,1]\) son \(x_k=\cos\left(\frac{(2k+1)\pi}{2(m+1)}\right)\), los nodos que minimizan la norma infinito de la parte del error de interpolación que depende de dichos nodos serán: \[ t_k =\frac{1}{2}((b-a)x_k+a+b)=\frac{1}{2}\left((b-a)\cos\left(\frac{(2k+1)\pi}{2(m+1)}\right)+a+b\right), \] para \(k=0,1,\ldots,m\).

El valor de la norma infinito \(\|(t-t_0)(t-t_1)\cdots (t-t_m)\|_\infty\) en los nodos anteriores será: \[ \begin{align*} & \|(t-t_0)(t-t_1)\cdots (t-t_m)\|_\infty =\max_{t\in [a,b]}|(t-t_0)(t-t_1)\cdots (t-t_m)|\\ & \quad =\max_{t\in [a,b]}\prod_{k=0}^m |t-t_k|=\max_{x\in [-1,1]}\prod_{k=0}^m \left|\frac{1}{2}((b-a)x+a+b)-\frac{1}{2}((b-a)x_k+a+b)\right|\\ & \quad = \left(\frac{b-a}{2}\right)^{m+1}\max_{x\in [-1,1]}\prod_{k=0}^m |x-x_k|=\left(\frac{b-a}{2}\right)^{m+1}\|\tilde{T}_{m+1}\|_\infty \\ & \quad =\left(\frac{b-a}{2}\right)^{m+1}\cdot\frac{1}{2^{m}}=\frac{(b-a)^{m+1}}{2^{2m+1}}. \end{align*} \]