En este capítulo vamos a aproximar una función \(f\) cualquiera por un miembro de la familia más conocida y más sencilla de tratar: los polinomios.

Recordemos que un polinomio de grado \(n\) tiene una expresión de la forma siguiente: \[ P_n(x)=a_n x^n+a_{n-1}x^{n-1}+\cdots +a_1 x+a_0 = \sum_{i=0}^n a_i x^i, \] con \(a_n\neq 0\) para que tenga grado \(n\).

Los valores \(a_i\), \(i=0,1,\ldots,n\) se denominan coeficientes del polinomio \(P_n\).

Cuanto “mejor” sea la función \(f\) a aproximar, es decir, cuanto más alto sea el valor de \(k\) donde \(f\in {\cal C}^k\), mejor control sobre el error cometido en la aproximación tendremos.

Para justificar la aproximación de una función \(f\) por polinomios veamos el Teorema de Weierstrass que dice básicamente que cualquier función continua puede aproximarse por un polinomio con un error tan pequeño como se quiera:

Teorema de Weierstrass. Sea \(f\in {\cal C}^0[a,b]\) una función continua en un intervalo \([a,b]\). Entonces dado un valor \(\epsilon >0\), existe un polinomio \(P_n(x)\) de un grado determinado \(n\) tal que: \[ |f(x)-P_n(x)|<\epsilon, \] para todo valor \(x\in [a,b]\).

Observación. La aproximación anterior es uniforme en el sentido de que el “error” \(\epsilon\) no depende del valor \(x\in [a,b]\).

El problema que intentamos resolver es el siguiente:

Problema: Dados \(n\) valores \((x_0,y_0),(x_1,y_1),\ldots,(x_n,y_n)\), hallar el polinomio \(P_n\) de grado mínimo tal que \(P_n(x_i)=y_i\), \(i=0,\ldots,n\).

Es decir, dados \(n+1\) puntos en el plano, hallar un polinomio de grado mínimo, (más adelante veremos que dicho grado es \(n\)) tal que \(P_n(x_i)=y_i\), \(i=0,\ldots,n\).

Observaciones:

• Si los puntos son parte de la gráfica de una función \(f\), entonces \(y_i=f(x_i)\) y las condiciones que debe verificar el polinomio \(P_n(x)\) son \(P_n(x_i)=f(x_i)\), \(i=0,\ldots,n\).
• Tenemos en total \(n+1\) condiciones, por tanto, el número de incógnitas debe ser \(n+1\). Pensemos que un polinomio de grado \(n\) tiene en total \(n+1\) coeficientes.

Sea pues \(P_n(x)=a_0+a_1 x+\cdots +a_n x^n\) el polinomio que queremos hallar.

Las condiciones \(P_n(x_i)=y_i\) serían las siguientes en función de los coeficientes \(a_i\), \(i=0,\ldots, n\): \[ a_0+a_1 x_i+\cdots +a_n x_i^n=y_i, \ i=0,\ldots, n. \]

Los coeficientes \(a_i\) deben verificar el siguiente sistema de ecuaciones lineal: \[ \left. \begin{align*} a_0+a_1 x_0+\cdots +a_n x_0^n = & y_0,\\ a_0+a_1 x_1+\cdots +a_n x_1^n = & y_1,\\ \vdots & \\ a_0+a_1 x_n+\cdots +a_n x_n^n = & y_n. \end{align*} \right\} \]

Veamos que el sistema lineal anterior tiene solución única.

El determinante del sistema es el siguiente: \[ \left| \begin{array}{cccc} 1 & x_0 & \cdots & x_0^n\\ 1 & x_1 & \cdots & x_1^n\\ \vdots &\vdots & \vdots & \vdots\\ 1 & x_n & \cdots & x_n^n \end{array} \right|. \]

El determinante anterior se llama determinante de Vandermonde y su valor vale: \[ \prod_{0\leq i <j\leq n} (x_i-x_j). \] Por tanto, si \(x_i\neq x_j\) para \(i\neq j\), el determinante del sistema no será cero y tendremos solución única para nuestro problema.

En resumen, tenemos el teorema siguiente:

Teorema. Sean \((x_0,y_0),(x_1,y_1),\ldots,(x_n,y_n)\), \(n\) valores con \(x_i\neq x_j\), para \(i\neq j\), (es decir, las abscisas son todas diferentes). Entonces existe un único polinomio \(P_n(x)\) de grado \(n\) tal que \(P_n(x_i)=y_i\), \(i=0,\ldots,n\).

Para hallar el polinomio \(P_n(x)\) deberíamos resolver el sistema anterior para calcular los coeficientes del polinomio \(a_i\), \(i=0\ldots, n\).

Resolver el sistema lineal anterior es muy tedioso, sobre todo si el número de nodos \(n+1\) es grande. Por dicho motivo, vamos a ver técnicas para obtener el polinomio interpolador sin tener que resolver dicho sistema.

Una de dichas técnicas consiste en escribir el polinomio interpolador de forma explícita en función de unos polinomios especiales denominados polinomios interpoladores de Lagrange:

Polinomios de Lagrange. Sean \(x_0,\ldots, x_n\), \(n+1\) nodos donde suponemos que \(x_i\neq x_j\), para \(i\neq j\). Se define el polinomio de Lagrange \(L_{n,k}(x)\) de grado \(n\) asociado al nodo \(x_k\) de la forma siguiente: \[ \begin{array}{rl} L_{n,k}(x)= & \displaystyle\frac{(x-x_0)\cdots (x-x_{k-1})\cdot (x-x_{k+1})\cdots (x-x_n)}{(x_k-x_0)\cdots (x_k-x_{k-1})\cdot (x_k-x_{k+1})\cdots (x_k-x_n)},\\ = & \displaystyle\frac{\prod_{i\neq k} (x-x_i)}{\prod_{i\neq k} (x_k-x_i)}. \end{array} \]

Los polinomios de Lagrange verifican la proposición siguiente:

Proposición. Sean \(x_0,\ldots, x_n\), \(n+1\) nodos donde suponemos que \(x_i\neq x_j\), para \(i\neq j\). Sean \(L_{n,k}(x)\) el polinomio de Lagrange de grado \(n\) asociado al nodo \(x_k\). Entonces dicho polinomio verifica que: \[ L_{n,k}(x_i)=0,\mbox{ si }i\neq k,\ L_{n,k}(x_k)=1. \]

El Teorema siguiente nos dice cómo calcular el polinomio interpolador a partir de los polinomios de Lagrange:

Teorema.

Sean \(n+1\) puntos \((x_0,y_0),(x_1,y_1),\ldots,(x_n,y_n)\) con \(x_i\neq x_j\), con \(i\neq j\). Entonces el polinomio interpolador \(P_n(x)\) se puede expresar de la forma siguiente: \[ \begin{array}{rl} P_n(x)= & y_0\cdot L_{n,0}(x)+y_1\cdot L_{n,1}(x)+\cdots + y_n\cdot L_{n,n}(x)\\ = & \displaystyle\sum_{k=0}^n y_k\cdot L_{n,k}(x), \end{array} \] donde \(L_{n,k}(x)\) es el polinomio de Lagrange correspondiente al nodo \(x_k\), \(k=0,\ldots,n\).

Observación:

Si los \(n+1\) puntos forman parte de la gráfica de una función \(f\), es decir \(y_i=f(x_i)\), para \(i=0.\ldots,n\), entonces el polinomio interpolador se escribiría de la forma siguiente: \[ \begin{array}{rl} P_n(x)= & f(x_0)\cdot L_{n,0}(x)+f(x_1)\cdot L_{n,1}(x)+\cdots + f(x_n)\cdot L_{n,n}(x)\\ = & \displaystyle\sum_{k=0}^n f(x_k)\cdot L_{n,k}(x), \end{array} \]

Interpolar una función \(f\) en unos nodos determinados puede interpretarse como una manera de aproximar la función \(f\) en un entorno de los nodos, es decir, en un dominio que esté relativamente cerca de dichos nodos.

Es fundamental estimar de alguna manera el error cometido en un valor cualquiera cuando se intenta aproximar una función. El Teorema siguiente nos da una expresión del error cometido cuando aproximamos una función \(f\) por un polinomio interpolador:

Teorema. Error en la interpolación.

Sea \(f\in {\cal C}^{n+1}[a,b]\) una función de clase \(n+1\) en un intervalo \([a,b]\). Consideramos \(n+1\) nodos \(x_0, x_1,\ldots,x_n\) en dicho intervalo \([a,b]\) con \(x_i\neq x_j\), con \(i\neq j\). Sea \(P_n(x)\) el polinomio interpolador de la función \(f\) en los puntos \((x_0,f(x_0)),(x_1,f(x_1)),\ldots,(x_n,f(x_n))\).

Sea \(x\) un valor cualquiera dentro del intervalo \([a,b]\). Entonces existe un valor \(\xi(x)\in <x_0,\ldots,x_n,x>=(\min\{x_0,\ldots,x_n,x\},\max\{x_0,\ldots,x_n,x\})\) tal que: \[ f(x)-P_n(x)=\frac{f^{(n+1)}(\xi(x))}{(n+1)!}\cdot (x-x_0)\cdot (x-x_1)\cdots (x-x_n). \]

Observación:

El error cometido \(f(x)-P_n(x)\) cuando interpolamos una función \(f\) tiene tres partes bien diferenciadas:

• El producto \((x-x_0)\cdot (x-x_1)\cdots (x-x_n)\) nos dice que cuándo más cerca esté el valor \(x\) de los nodos, mejor será la aproximación.
• El denominador \((n+1)!\) nos dice que cuantos más nodos tengamos en cuenta, mejor será la aproximación y menor el error cometido.
• La “caja negra” \(f^{(n+1)}(\xi(x))\) que depende de la función \(f\) nos dice que tenemos que tener “controladas” las derivadas de la función \(f\) ya que si éstas van aumentando sin control, no tendrá ningun sentido interpolar ya que el error cometido puede aumentar indefinidamente.

Muchas aplicaciones derivadas de la interpolación consisten en la interpolación de datos tabulados, es decir, datos que tenemos en una tabla de la forma:

En estos casos, no necesitamos conocer explícitamente la expresión del polinomio interpolador.

En el caso de datos tabulados, al no conocer la función \(f\) que interpolamos, no podemos aplicar la fórmula del error en la interpolación y no podemos saber cuál es el error cometido en un valor cualquiera \(x\).

El método de Lagrange para hallar el polinomio interpolador tiene el siguiente problema:

• Imaginemos que hemos hallado el polinomio interpolador para los puntos \((x_0,f(x_0)),(x_1,f(x_1)),\ldots,(x_n,f(x_n))\) que forman parte de la gráfica de una cierta función \(f\).

• De cara a hallar una mejor aproximación de la función \(f\), añadimos un nuevo punto \((x_{n+1},f(x_{n+1}))\). Entonces para hallar el nuevo polinomio interpolador en los puntos \((x_0,f(x_0)),(x_1,f(x_1)),\ldots,(x_n,f(x_n)),(x_{n+1},f(x_{n+1}))\) usando el método de Lagrange, tenemos que “empezar de nuevo” ya que los polinomios de Lagrange hallados previamente para los puntos \((x_0,f(x_0)),(x_1,f(x_1)),\ldots,(x_n,f(x_n))\) no nos sirven y hay que volver a hallar los “nuevos polinomios de Lagrange”.

Lo dicho anteriormente hace que el método de Lagrange no sea apropiado para hallar el polinomio interpolador para datos tabulados.

Vamos a ver un método nuevo que aprovecha el trabajo realizado en interpolar los puntos \((x_0,f(x_0)),(x_1,f(x_1)),\ldots,(x_n,f(x_n))\) de cara a interpolar los puntos anteriores añadiendo un nuevo punto \((x_{n+1},f(x_{n+1}))\). Dicho método se llama método de Neville.

En primer lugar necesitamos unas definiciones previas:

Definición.

Sea \(f\) una función definida en los valores \(x_0,x_1,\ldots,x_n\). Sean \(m_1,m_2,\ldots,m_k\), \(k\) enteros diferentes, con \(m_i\in\{0,1,\ldots,n\}\), \(i=1,\ldots,k\). El polinomio interpolador de Lagrange que interpola los puntos \((x_{m_1},f(x_{m_1})),(x_{m_2},f(x_{m_2})),\ldots,(x_{m_k},f(x_{m_k}))\) se denota por \(P_{m_1,m_2,\ldots,m_k}(x)\).

El siguiente resultado nos dice cómo hallar el valor del polinomio interpolador en un punto \(x\) a partir de los Polinomios interpoladores de Lagrange introducidos anteriormente:

Proposición.

Sean \(x_0,\ldots,x_n\), \(n+1\) nodos, sea \(f\) una función y sea \(P_n(x)\) el polinomio interpolador de la función \(f\) en los nodos anteriores.

Sean \(x_i\) y \(x_j\) dos nodos cualquiera con \(i\neq j\). Sean \(P_{0,1,\ldots,i-1,i+1,\ldots,n}(x)\) y \(P_{0,1,\ldots,j-1,j+1,\ldots,n}(x)\) los Polinomios interpoladores de Lagrange en los nodos \(x_0,x_1,\ldots,x_{i-1},x_{i+1},\ldots,x_n\) y \(x_0,x_1,\ldots,x_{j-1},x_{j+1},\ldots,x_n\), respectivamente.

Proposición.(continuación)

Entonces el polinomio interpolador \(P_n(x)\) que interpola todos los nodos \(x_0,x_1,\ldots,x_n\) puede escribirse de la forma siguiente: \[ P_n(x)=\frac{(x-x_j)P_{0,1,\ldots,j-1,j+1,\ldots,n}(x)-(x-x_i)P_{0,1,\ldots,i-1,i+1,\ldots,n}(x)}{x_i-x_j}. \]

El método de Neville nos permite hallar los distintos polinomios de interpolación de forma recursiva.

• Los polinomios interpoladores de Lagrange de grado \(0\) serían: \[ P_0 =f(x_0),\ P_1=f(x_1),\ldots, P_n=f(x_n), \] en general \(P_i=f(x_i)\), \(i=0,\ldots,n\).

• La proposición anterior nos da los polinomios interpoladores de Lagrange de grado \(1\): \[ \begin{align*} P_{0,1}(x)= & \frac{(x-x_0)P_1-(x-x_1)P_0}{x_1-x_0},\\ P_{1,2}(x)= & \frac{(x-x_1)P_2-(x-x_2)P_1}{x_2-x_1},\\ \vdots &\\ P_{n-1,n}(x) = & \frac{(x-x_{n-1})P_n-(x-x_n)P_{n-1}}{x_n-x_{n-1}}. \end{align*} \]

• Usando la misma técnica podemos obtener los polinomios interpoladores de Lagrange de grado \(2\):

\[ \begin{align*} P_{0,1,2}(x)= & \frac{(x-x_0)P_{1,2}(x)-(x-x_2)P_{0,1}(x)}{x_2-x_0},\\ P_{1,2,3}(x)= & \frac{(x-x_1)P_{2,3}(x)-(x-x_3)P_{1,2}(x)}{x_3-x_1},\\ \vdots &\\ P_{n-2,n-1,n}(x) = & \frac{(x-x_{n-2})P_{n-1,n}(x)-(x-x_n)P_{n-2,n-1}(x)}{x_n-x_{n-2}}. \end{align*} \]

• Y así sucesivamente.

Los polinomios generados anteriormente se pueden escribir en forma de tabla de la siguiente manera:

La manera de indicar los polinomios de Lagrange interpoladores es un poco tediosa ya que el número de subíndices va aumentando a medida que aumenta el número de nodos a interpolar.

Como los subíndices son consecutivos, basta tener en cuenta el subíndice del nodo en que empieza la interpolación y el número de nodos que se interpola.

Por ejemplo, para indicar el polinomio interpolador de Lagrange \(P_{2,3,4}\) bastaría “guardar” el valor \(4\), es decir, el último nodo que se interpola \(x_4\) y el número \(3\) que es el número de nodos que se interpola ya que con estos dos valores ya quedaría claro que los nodos a interpolar serían \(x_2,x_3\) y \(x_4\).

Por dicho motivo, introduciremos una nueva notación para indicar los polinomios interpoladores de Lagrange: para indicar el polinomio \(P_{i-j,i-j+1,\ldots,i}\) definimos: \[ Q_{i,j}:=P_{i-j,i-j+1,\ldots,i}. \]

La tabla anterior sería en la nueva notación:

Los polinomios interpoladores de Lagrange rara vez se calculan.

Lo que habitualmente se hace es usar el algoritmo anterior para evaluar el valor de \(P_n(x)\) para un valor determinado \(x\).

De esta manera, la tabla de los \(Q_{i,j}\) no son polinomios sino valores.

De hecho, el valor de \(Q_{i,i}\) es el valor de polinomio interpolador en los valores \((x_0,y_0),\ldots (x_i,y_i)\) en el punto \(x\).

Por tanto, \(Q_{n,n}(x)\) es el valor de \(P_n(x)\).

• INPUT x, x_0, x_1,..., x_n; (damos los valores \(x\) y los nodos \(x_0,x_1,\ldots,x_n\))
• INPUT f(x_0), f(x_1),..., f(x_n) (damos los valores \(f(x_0),f(x_1),\ldots, f(x_n)\))
• For i=0,...,n
  • Set Q_{i,0}=f(x_i). (guardamos los valores \(f(x_i)\) en \(Q_{i,0}\) para \(i=0,1,\ldots,n\))

• For i=1,...,n
  • For j=1,...,i
    • Set Q_{i,j}=
    ((x-x_{i-j})*Q_{i,j-1}-(x-x_i)*Q_{i-1,j-1})/(x_i-x_{i-j}). (calculamos los valores \(Q_{i,j}\) de la tabla)
• Print Q_{n,n} (damos el valor de \(P_n(x)\))

El método de Neville anteriormente descrito se usa básicamente para hallar el valor del polinomio interpolador \(P_n(x)\) en un valor concreto \(x\).

El método de las diferencias divididas permite hallar el polinomio interpolador \(P_n(x)\) de la función \(f\) en los nodos \(x_0,x_1,\ldots,x_n\) expresando \(P_n(x)\) de la forma siguiente: \[ \begin{align*} P_n(x)= & a_0+a_1 (x-x_0)+a_2 (x-x_0)(x-x_1)+\cdots \\ & + a_n (x-x_0)\cdots (x-x_{n-1}), \end{align*} \] es decir escribimos el polinomio interpolador como combinación lineal de los elementos de la base siguiente del espacio vectorial de polinomios de grado \(n\): \[ 1,\ (x-x_0),\ (x-x_0)(x-x_1),\ \ldots, (x-x_0)\cdots (x-x_{n-1}). \]

Calculemos los primeros coeficientes \(a_i\).

• Como \(P_n(x_0)=f(x_0)\), deducimos que \(a_0=f(x_0)\),
• Como \(P_n(x_1)=a_0+a_1 (x_1-x_0)=f(x_1)\), deducimos \(a_1=\frac{f(x_1)-f(x_0)}{x_1-x_0}\).

Para calcular los demás coeficientes \(a_i\), necesitamos introducir las diferencias divididas de cierto orden:

Definición de diferencias divididas.

Las diferencias divididas de orden \(0\) valen \(f[x_i]=f(x_i)\).

Las diferencias divididas de orden \(1\) valen \(f[x_i,x_{i+1}]=\frac{f[x_{i+1}]-f[x_i]}{x_{i+1}-x_i}\).

Las diferencias divididas de orden \(2\) valen \(f[x_i,x_{i+1},x_{i+2}]=\frac{f[x_{i+1},x_{i+2}]-f[x_i,x_{i+1}]}{x_{i+2}-x_i}\).

En general, las diferencias divididas de orden \(k\) se definen en función de las diferencias divididas de orden \(k-1\): \[ f[x_i,\ldots,x_{i+k}]=\frac{f[x_{i+1},\ldots,x_{i+k}]-f[x_i,\ldots, x_{i+k-1}]}{x_{i+k}-x_i}. \]

Entonces tenemos el siguiente Teorema:

Teorema.

El polinomio interpolador \(P_n(x)\) de la función \(f\) en los nodos \(x_0,x_1\ldots,x_n\) vale en función de las diferencias divididas: \[ \begin{align*} P_n(x)= & f[x_0]+f[x_0,x_1](x - x_0)+\cdots \\ & + f[x_0,\ldots,x_n](x-x_0)\cdots (x-x_{n-1}),\\ = & f[x_0]+\sum_{k=1}^n f[x_0,\ldots,x_k](x-x_0)\cdots (x-x_{k-1}). \end{align*} \]

Observación.

El polinomio interpolador también puede escribirse como: \[ \begin{align*} P_n(x)= & f[x_n]+f[x_n,x_{n-1}](x - x_n)+f[x_n,x_{n-1},x_{n-2}](x-x_n)(x-x_{n-1})\cdots \\ & + f[x_n,\ldots,x_0](x-x_n)(x-x_{n-1})\cdots (x-x_1),\\ = & f[x_n]+\sum_{k=n-1}^{0} f[x_n,\ldots,x_{k}](x-x_n)(x-x_{n-1})\cdots (x-x_{k+1})\\ = & f[x_n]+\sum_{k=0}^{n-1} f[x_k,\ldots,x_{n}](x-x_{k+1})\cdots (x-x_{n-1})(x-x_n). \end{align*} \]

• INPUT x_0,x_1,...,x_n, f(x_0), f(x_1),...,f(x_n) (damos los valores de los nodos \(x_i\) y sus valores a interpolar \(f(x_i)\), \(i=0,1,\ldots,n\))
• For i=0,...,n
  • Set F_{i,0}=f(x_i). (Definimos \(F_{i,j}=f[x_{i-j},\ldots,x_i]\). Los coeficientes \(a_i=f[x_0,\ldots,x_n]\) serán \(F_{i,i}\). Por tanto, \(F_{i,0}=f[x_i]=f(x_i)\))
• For i=1,...,n
  • For j=1,...,i
    • Set F_{i,j}=(F_{i,j-1}-F_{i-1,j-1})/(x_i-x_{i-j}). (Calculamos \(F_{i,j}=f[x_{i-j},\ldots,x_i]\) en función de \(F_{i,j-1}=f[x_{i-j+1},\ldots,x_i]\) y \(F_{i-1,j-1}=f[x_{i-j},\ldots,x_{i-1}]\))
• Print F_{0,0},...,F_{n,n}. (damos los coeficientes \(a_0\ldots, a_n\))

El teorema siguiente nos da la relación entre las diferencias divididas y la función que interpolamos \(f\):

Teorema.

Sea \(f\in{\cal C}^n\) una función de classe \({\cal C}^n\) en un intervalo \([a,b]\) que contenga los nodos a interpolar \(x_0,x_1,\ldots, x_n\). Entonces existe un número \(\xi\in (a,b)\) tal que \[ f[x_0,x_1,\ldots,x_n]=\frac{f^{(n)}(\xi)}{n!}. \]

Supongamos que los nodos son equiespaciados, esto es, que la diferencia entre dos nodos consecutivos es constante, o \(x_{i}-x_{i-1}=h\), donde \(h\) es independiente del nodo considerado.

En este caso, podemos escribir los nodos de la siguiente manera, suponiendo que el primero \(x_0\) es el menor de todos ellos: \[ x_i=x_0+i\cdot h,\ i=0,1,\ldots,n. \] Así por ejemplo, \(x_1=x_0+h\), \(x_2=x_0+2h\), y así sucesivamente hasta llegar a \(x_n=x_0+n\cdot h.\)

Vamos a hallar una expresión del polinomio interpolador en los nodos anteriores para una determinada función \(f\) en una forma simplificada.

Recordemos que el polinomio interpolador se escribía como: \[ P_n(x)= f[x_0]+\sum_{k=1}^n f[x_0,\ldots,x_k](x-x_0)\cdots (x-x_{k-1}). \]

Sea \(x\) un valor cualquiera. Asociamos a este valor \(x\) un valor \(s\) tal que \(x=x_0+s\cdot h\), o, si se quiere, \(s=\frac{x-x_0}{h}\). El polinomio interpolador anterior valdrá en función de \(s\): \[ \begin{align*} P_n(x)=& P_n(x_0+sh) \\ = & f[x_0]+\sum_{k=1}^n f[x_0,\ldots,x_k](x_0+sh-x_0)\cdots (x_0+sh-(x_0+(k-1)h)),\\ = & f[x_0]+\sum_{k=1}^n f[x_0,\ldots,x_k] s(s-1)\cdots (s-k+1)\cdot h^k. \end{align*} \]

Si definimos el número binomial \(\binom{s}{k}\) como \(\binom{s}{k}=\frac{s(s-1)\cdots (s-k+1)}{k!}\), podemos escribir el polinomio interpolador en los nodos equiespaciados como: \[ P_n(x)= f[x_0]+\sum_{k=1}^n \binom{s}{k} k! h^k f[x_0,\ldots,x_k], \] donde \(s=\frac{x-x_0}{h}\).

A continuación, vamos a escribir la parte \(f[x_0,\ldots,x_k]\) de una forma más compacta usando las diferencias hacia adelante introducidas en el capítulo de ceros cuando se explicó el Método de Aitken.

Definición de diferencias hacia delante.

Sean \(x_0,x_1\ldots, x_n\), \(n+1\) nodos equiespaciados y \(f\) una función que queremos interpolar. Definimos las diferencias \(\Delta^k f(x_i)\) de forma recurrente de la siguiente manera: \[ \Delta^0 f(x_i)=f(x_i),\ i=0,\ldots, n,\ \Delta^k f(x_i)=\Delta^{k-1}f(x_{i+1})-\Delta^{k-1}f(x_{i}). \]

Así por ejemplo, \[ \begin{align*} \Delta^0 f(x_0)= & f(x_0),\\ \Delta^1 f(x_0)= & \Delta f(x_0)=f(x_1)-f(x_0),\\ \Delta^2 f(x_0)= & \Delta f(x_1)-\Delta f(x_0)=f(x_2)-f(x_1)-(f(x_1)-f(x_0))\\ = & f(x_2)-2f(x_1)+f(x_0). \end{align*} \]

Los valores de \(f[x_0,x_1]\) y \(f[x_0,x_1,x_2]\) serán en función de las diferencias hacia adelante \(\Delta f(x_0)\) y \(\Delta^2 f(x_0)\): \[ \begin{align*} f[x_0,x_1]= & \frac{f(x_1)-f(x_0)}{x_1-x_0}=\frac{1}{h}\Delta f(x_0),\\ f[x_0,x_1,x_2]= & \frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}=\frac{\frac{1}{h}\left(\Delta f(x_1)-\Delta f(x_0)\right)}{2h}\\ = & \frac{1}{2h^2}\Delta^2 f(x_0). \end{align*} \] En general, puede demostrarse por inducción que: \[ f[x_0,x_1,\ldots,x_k]=\frac{1}{k!h^k}\Delta^k f(x_0). \]

Usando la expresión anterior, el polinomio de interpolación de la función \(f\) en los nodos equiespaciados se puede escribir como: \[ P_n(x)= f[x_0]+\sum_{k=1}^n \binom{s}{k} \Delta^k f(x_0). \]

Vamos a hacer los mismo pero en lugar de introducir diferencias hacia adelante, introduciremos las diferencias hacia atrás:

Definición de diferencias hacia atrás

Sean \(x_0,x_1\ldots, x_n\), \(n+1\) nodos equiespaciados y \(f\) una función que queremos interpolar. Definimos las diferencias \(\nabla^k f(x_i)\) de forma recurrente de la siguiente manera: \[ \nabla^0 f(x_i)=f(x_i),\ i=0,\ldots, n,\ \nabla^k f(x_i)=\nabla^{k-1}f(x_{i})-\nabla^{k-1}f(x_{i-1}). \]

Los valores de \(f[x_n,x_{n-1}]\) y \(f[x_n,x_{n-1},x_{n-2}]\) serán en función de las diferencias hacia atrás \(\nabla f(x_n)\) y \(\nabla^2 f(x_n)\): \[ \begin{align*} f[x_n,x_{n-1}]= & \frac{f(x_{n-1})-f(x_n)}{x_{n-1}-x_n}=\frac{1}{h}\nabla f(x_n),\\ f[x_n,x_{n-1},x_{n-2}]= & \frac{f[x_{n-1},x_{n-2}]-f[x_n,x_{n-1}]}{x_{n-2}-x_n}=\frac{\frac{1}{h}\left(\nabla f(x_{n-1})-\nabla f(x_n)\right)}{-2h}\\ = & \frac{1}{2h^2}\nabla^2 f(x_n). \end{align*} \]

En general, puede demostrarse por inducción que: \[ f[x_n,x_{n-1},\ldots,x_k]=\frac{1}{(n-k)!h^{n-k}}\nabla^{n-k} f(x_n). \]

Recordemos que el polinomio de interpolación de la función \(f\) en los nodos equiespaciados se podía expresar como: \[ P_n(x)= f[x_n]+\sum_{k=0}^{n-1} f[x_k,\ldots,x_{n}](x-x_{k+1})\cdots (x-x_{n-1})(x-x_n). \]

Dado un valor \(x\), sea \(s\) tal que \(x=x_n-sh\), donde \(h=x_i-x_{i-1}\).

Entonces \(x_i=x_n-(n-i)h\) y \(x-x_i=x_n-sh-(x_n-(n-i)h)=(n-i-s)h\).

El productorio \(\pi_{k,n}=(x-x_{k+1})\cdots (x-x_{n-1})(x-x_n)\) se puede expresar como: \[ \begin{align*} \pi_{k,n}= & (n-(k+1)-s)\cdots (n-(n-1)-s)(n-n-s)h^{n-k}\\ = & (-1)^{n-k}h^{n-k}s (s-1)\cdots (s-n+k+1)\\ = & (-1)^{n-k}h^{n-k}(n-k)!\binom{s}{n-k}. \end{align*} \]

El polinomio interpolador será: \[ \begin{align*} P_n(x)= & f[x_n]+\sum_{k=0}^{n-1}f[x_k,\ldots,x_{n}](-1)^{n-k}h^{n-k}(n-k)!\binom{s}{n-k}\\ = & f[x_n]+\sum_{k=0}^{n-1}(-1)^{n-k}\binom{s}{n-k}\nabla^{n-k} f(x_n)\\ = & f[x_n]+\sum_{k=1}^{n}(-1)^{k}\binom{s}{k}\nabla^{k} f(x_n). \end{align*} \]

Vamos a generalizar el concepto de interpolación.

Hasta ahora, hemos estudiado cómo dada una función \(f\) y unos nodos \(x_0,x_1,\ldots, x_n\) distintos, hallar el único polinomio interpolador \(P_n(x)\) en los nodos anteriores, es decir, \(P_n(x_i)=f(x_i)\), \(i=0,1,\ldots, n\).

Imaginemos que queremos hacer lo mismo pero no nos “conformamos” en que en el nodo \(x_i\), la función \(f(x)\) y el polinomio \(P_n(x)\) coincidan, sino que exigimos que dichas funciones coincidan en sus derivadas en el nodo \(x_i\) hasta cierto orden \(m_i\), donde para cada nodo \(x_i\), tendremos un orden de coincidencia \(m_i\), para \(i=0,1,\ldots, n\).

A dicho polinomio se le denomina polinomio osculador de aproximación de \(f\) en los nodos \(x_0,x_1,\ldots,x_n\) con órdenes de coincidencia \(m_0,m_1,\ldots, m_n\):

Definición de polinomio osculador.

Sea \(x_0,x_1,\ldots,x_n,\) \(n+1\) nodos diferentes y sean \(m_0,m_1,\ldots,m_n,\) \(n+1\) enteros positivos. Sea \(N=\displaystyle\sum_{i=0}^n m_i\). Consideremos una función \(f\in {\cal C}^m\) de clase \({\cal C}^m\), con \(m=\max\{m_0,m_1,\ldots,m_n\}\). Entonces el polinomio osculador \(P_{N+n}(x)\) de grado menor o igual que \(N+n\) es el que verifica: \(P_{N+n}^{(j)}(x_i)=f^{(j)}(x_i)\), \(j=0,1,\ldots, m_i\), \(i=0,1,\ldots,n\), con \(f^{(0)}=f\).

Observación.

Para calcular el polinomio osculador tenemos en total \(N=(m_0+1)+(m_1+1)+\cdots + (m_n+1)=N+n+1\) condiciones. Por tanto, el polinomio osculador tendrá grado \(N+n\) al tener en total \(N+n+1\) coeficientes.

En este curso vamos estudiar en detalle la interpolación de Hermite para el caso en que quisiéramos interpolar la función \(f\) y su derivada \(f'\).

Es decir, estudiaremos el caso en que \(m_0=m_1=\cdots =m_n=1\).

El grado del polinomio interpolador será \(N+n=1+\cdots +1+n=n+1+n=2n+1\).

El problema que nos planteamos, pues, es el siguiente:

Problema:

Sean \(x_0,x_1,\ldots, x_n\), \(n+1\) nodos distintos. Sea \(f\in {\cal C}^1\) una función de clase \({\cal C}^1\), es decir, la función y su primera derivada son funciones continuas.

El problema es hallar un polinomio \(P_{2n+1}(x)\) de grado menor o igual que \(2n+1\) tal que la función \(f\) y el polinomio anterior \(P_{2n+1}\) coincidan en los puntos \(x_0,x_1,\ldots, x_n\) y tengan la misma recta tangente dichos puntos, o que las pendientes de las rectas tangentes coincidan: \[ P_{2n+1}(x_i)=f(x_i),\ P'_{2n+1}(x_i)=f'(x_i),\ i=0,1,\ldots,n. \]

El teorema siguiente nos da la expresión del polinomio \(P_{2n+1}(x)\):

Teorema. Sea \(f\in {\cal C}^1[a,b]\) una función de clase \({\cal C}^1\) en un intervalo \([a,b]\) y sean \(x_0,x_1,\ldots,x_n\), \(n+1\) valores distintos en \([a,b]\), \(x_i\in [a,b]\), \(i=0,1,\ldots,n\). El único polinomio \(P_{2n+1}(x)\) de grado menor o igual que \(2n+1\) que coincide con \(f\) y \(f'\) en los nodos anteriores es el siguiente: \[ P_{2n+1}(x)=\sum_{j=0}^n f(x_j)H_{n,j}(x)+\sum_{j=0}^n f'(x_j)\hat{H}_{n,j}(x), \] donde los polinomios \(H_{n,j}(x)\) y \(\hat{H}_{n,j}(x)\) son los siguientes (\(L_{n,j}(x)\) es el \(j\)-ésimo polinomio de Lagrange de grado \(n\)): \[ \begin{align*} H_{n,j}(x)= & (1-2(x-x_j)L'_{n,j}(x_j)) L_{n,j}^2(x),\\ \hat{H}_{n,j}(x) = & (x-x_j)L_{n,j}^2 (x). \end{align*} \]

Teorema (continuación).

Además, si \(f\in {\cal C}^{2n+2}[a,b]\), entonces el error cometido \(f(x)-P_{2n+1}(x)\) en un punto \(x\in [a,b]\) vale: \[ f(x)-P_{2n+1}(x)=\frac{(x-x_0)^2\cdots (x-x_n)^2}{(2n+2)!}f^{(2n+2)}(\xi(x)), \] donde \(\xi(x)\) es un valor que está en el intervalo \((a,b)\): \(\xi(x)\in (a,b)\).

Observación.

El polinomio del teorema anterior tiene como máximo grado \(2n+1\) ya que como el grado de \(L_{n,j}(x)=n\), el grado de \(H_{n,j}(x)\) será \(1+2n\) y el de \(\hat{H}_{n,j}(x)\), también. Por tanto, como \(P_{2n+1}(x)\) es una combinación lineal de los polinomios anteriores, tendrá grado \(2n+1\).

El método anterior para el cálculo del polinomio interpolador de Hermite es muy tedioso ya que hay que calcular en total

• \(n+1\) polinomios de Lagrange \(L_{n,j}(x)\),
• \(n+1\) polinomios \(H_{n,j}(x)\),
• \(n+1\) polinomios \(\hat{H}_{n,j}(x)\),

lo que hace un total de \(3(n+1)\) polinomios!

Otra forma de calcular el polinomio interpolador de Hermite es usar una fórmula parecida a la que teníamos para calcular el polinomio interpolador usando el método de Newton: \[ P_n(x)= f[x_0]+\sum_{k=1}^n f[x_0,\ldots,x_k](x-x_0)\cdots (x-x_{k-1}). \]

Para hallar una expresión parecida a la anterior para la interpolación de Hermite dados los nodos \(x_0,x_1,\ldots, x_n\), introducimos los nodos \(z_i\), \(i=0,1,\ldots,2n+1\) de la forma siguiente: \[ z_{2i}=z_{2i+1}=x_i,\ i=0,1,\ldots, n. \] Es decir, \[ z_0=x_0,\ z_1=x_0,\ z_2=x_1,\ z_3=x_1,\ldots, z_{2n}=x_n,\ z_{2n+1}=x_{n}. \] La idea es repetir cada nodo \(x_i\) dos veces.

Ahora vamos a aplicar la fórmula de las diferencias divididas a los nuevos nodos \(z_i\).

El problema se plantea en hallar las diferencias divididas de primer orden \(f[z_{2i},z_{2i+1}]\) ya que nos queda \[ f[z_{2i},z_{2i+1}]=\frac{f[z_{2i+1}]-f[z_{2i}]}{z_{2i+1}-z_{2i}}=\frac{f[x_{i}]-f[x_{i}]}{x_{i}-x_{i}}=\frac{0}{0}! \] Entonces es cuando surge usar el papel que juega la derivada. Definimos \[ f[z_{2i},z_{2i+1}]=f'(z_{2i})=f'(x_i),\ i=0,1,\ldots,n. \]

Entonces el polinomio de interpolación de Hermite se escribe de la forma siguiente: \[ P_{2n+1}(x)=f[z_0]+\sum_{k=1}^{2n+1}f[z_0,\ldots,z_k](x-z_0)(x-z_1)\cdots (x-z_{k-1}). \]

• INPUT x_0,...x_n, f(x_0),...,f(x_n), f'(x_0),...,f'(x_n). (entramos los nodos \(x_i\), \(f(x_i)\) y \(f'(x_i)\), Guardaremos en la matriz
• For i=0,1,...,n (Definimos \(Q[i,j]=F[z_{i-j},\ldots, z_i]\). Los coeficientes \(a_i=f[z_0,\ldots,z_i]\) serán \(Q_{i,i}\). Por tanto, \(Q_{2i,0}=Q_{2i+1,0}=f[x_i]=f(x_i)\) y \(Q_{2i+1,1}=f'(x_i)\))
  • Set z[2i]=x_i;
  • Set z[2i+1]=x_i;
  • Set Q[2i+1,0]=f(x_i);
  • Set Q[2i+1,1]=f'(x_i);
  • If i!=0 then
    • Set Q[2i,1]=(Q[2i,0]-Q[2i-1,0])/(z[2i]-z[2i-1]).

• • For i=2,...,2n+1
    • For j=2,...,i
      • Set Q[i,j]=(Q[i,j-1]-Q[i-1,j-1])/(z[i]-z[i-j]).
• Print Q[0,0],...,Q[2n+1,2n+1].

La interpolación usando polinomios de grados grandes tiene el inconveniente que pueden oscilar mucho.

En la figura siguiente, podemos observar los polinomios de interpolación de grados \(10\) (en negro), \(20\) (en azul) y \(30\) (en verde) de la función \(f(x)=\frac{1}{1+25 x^2}\) en nodos equiespaciados en el intervalo \([-1,1]\).

Observamos que los polinomios de interpolación oscilan muchísimo en los extremos lo que hace que, en este caso, la interpolación polinomial tenga un error muy alto para valores de \(x\) cercanos a \(-1\) y a \(1\) y no sea viable.

A dicho fenómeno se le denomina Fenómeno de Runge.

Una alternativa para evitar este tipo de problemas es usar interpolación por splines o interpolación a trozos.

Dados \(n+1\) nodos \(x_0,x_1,\ldots,x_n\), dicha interpolación consiste en hallar en cada subintervalo \([x_i,x_{i+1}]\) un polinomio de grado bajo de tal manera que la función a interpolar tenga grado de derivabilidad lo más alto posible.

En esta sección veremos la interpolación cúbica ya que en este caso, se puede conseguir que la función interpoladora tenga grado de derivabilidad \(2\) o que sea de clase \({\cal C}^2\) en el intervalo considerado.

En muchas aplicaciones, como por ejemplo en informática gráfica, es suficiente que la función interpoladora sea de clase \({\cal C}^2\) ya que el ojo humano no distingue grados de derivabilidad superior.

La interpolación spline más sencilla es la lineal a trozos.

Dicha interpolación consiste en juntar de forma lineal los puntos \((x_0,f(x_0)), (x_1,f(x_1)),\ldots, (x_n,f(x_n))\). En la figura siguiente se muestra la interpolación lineal para \(n=6\) donde la función que se interpola está en color rojo.

El problema de la interpolación lineal es que la función interpoladora sólo es continua y no es ni siquiera derivable en los nodos \(x_i\).

Esto hace que la interpolación lineal sea inviable en aplicaciones donde se requiera diferenciabilidad de la función interpoladora.

Una manera de solventar el problema anterior es usar interpolación cúbica en lugar de interpolación lineal.

Es decir, dados \(f(x_i)\) y \(f'(x_i)\) para \(i=0,1,\ldots,n\) hallar un polinomio de interpolación de grado \(3\) para cada subintervalo usando interpolación de Hermite.

El problema es que muchas veces no se conoce el valor de \(f'(x_i)\). Por dicho motivo, vamos a explicar cómo realizar interpolación por polinomios de grado menor o igual que \(3\) sin tener conocimiento explícito de la derivada \(f'(x)\).

Dicho tipo de interpolación se denomina interpolación por splines cúbicos.

Definición de spline cúbico interpolador.

Sea \(f\) una función definida en un intervalo \([a,b]\) y sean \(a=x_0<x_1<\cdots <x_n=b\), \(n+1\) nodos en dicho intervalo. Una función spline cúbico interpolador \(S\) para \(f\) es una función que satisface las condiciones siguientes:

1. \(S(x)\) es un polinomio de grado menor o igual que \(3\), al que denotaremos por \(S_j(x)\) en el subintervalo \(x_i,x_{i+1}\), para \(i=0,1,\ldots,n-1\).
2. \(S_i(x_i)=f(x_i)\) y \(S_i(x_{i+1})=f(x_{i+1})\) para \(i=0,1,\ldots, n-1\).
3. \(S_{i+1}'(x_{i+1})=S_i'(x_{i+1})\) para \(i=0,1,\ldots, n-2\).
4. \(S_{i+1}''(x_{i+1})=S_i''(x_{i+1})\) para \(i=0,1,\ldots, n-2\).
5. Una de las condiciones siguientes debe cumplirse:
   1. \(S''(x_0)=S''(x_n)=0\), (frontera natural o libre).
   2. \(S'(x_0)=f'(x_0)\) y \(S'(x_n)=f'(x_n)\) (frontera fijada).

Observación.

A partir de la condición 2. se deduce que el spline cúbico interpolador es una función continua: \[ S_{i+1}(x_{i+1})=S_i(x_{i+1}) \mbox{ para }i=0,1,\ldots, n-2, \] y las condiciones 3. y 4. nos dicen que dicho spline es de clase \({\cal C}^2[a,b]\).

Observación.

Calcular el spline cúbico interpolador equivale a calcular en total \(n\) polinomios de grado menor o igual que \(3\), \(S_i(x)\), \(i=0,1,\ldots,n-1\).

Por tanto, tenemos en total \(4n\) incógnitas a calcular, \(4\) para cada polinomio \(S_i(x)\), \(i=0,1,\ldots,n-1\).

Veamos cuantas condiciones tenemos:

• Condición 2: \(2n\) condiciones.
• Condición 3: \(n-1\) condiciones.
• Condición 4: \(n-1\) condiciones.

En total tenemos \(2n+n-1+n-1=4n-2\) condiciones. Quedan, por tanto, \(2\) incógnitas “libres”. Las condiciones 5.1 o 5.2 nos fijan dichas incógnitas.

Vamos a dar un algoritmo para calcular los \(n\) polinomios de grado menor o igual que \(3\), \(S_j(x)\).

Recordemos que los nodos son los siguientes: \(a=x_0<x_1<\cdots < x_n=b\) y la función a interpolar era \(f\). Por tanto, nuestra función spline interpoladora debe pasar por los puntos \((x_i,f(x_i))\), para \(i=0,1,\ldots,n\).

En primer lugar, escribimos dichos polinomios de la forma siguiente: \[ S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i (x-x_i)^3, \] donde \(a_i,b_i,c_i\) y \(d_i\) son las incógnitas que debemos calcular para \(i=0,1,\ldots,n-1\).

• Como \(S_i(x_i)=f(x_i)\), tenemos que \(a_i=f(x_i)\) para \(i=0,1,\ldots,n-1\). Hemos reducido el número de incógnitas a \(3n\) ya que \(n\) de ellas ya están calculadas (las \(a_i\))
• Como \(a_{i+1}=f(x_{i+1})=S_{i+1}(x_{i+1})=S_i(x_{i+1})\), tenemos que: \[ a_i+b_i(x_{i+1}-x_i)+c_i(x_{i+1}-x_i)^2+d_i (x_{i+1}-x_i)^3 =a_{i+1}, \] para \(i=0,1,\ldots,n-2\).

Sea \(h_i=x_{i+1}-x_i\). La condición anterior será: \[ a_i+b_i h_i+c_i h_i^2+d_i h_i^3=a_{i+1},\ i=0,1,\ldots,n-2. \] Si definimos \(a_n=f(x_n)\), la condición anterior será cierta para \(i=0,1,\ldots,n-1\) ya que la condición para \(i=n-1\) equivale a decir que \(S_{n-1}(x_n)=f(x_n)\).

Definimos \(b_n=S'(x_n)=S'_{n-1}(x_n)\).

Los valores de \(S_i'(x)\) son: \[ S_i'(x)=b_i+2c_i(x-x_i)+3d_i(x-x_i)^2. \] Por tanto \(S_i'(x_i)=b_i\), para \(i=0,1,\ldots,n-1\).

Si expresamos la condición \(S_{i}'(x_{i+1})=S_{i+1}'(x_{i+1})=b_{i+1}\), tenemos que: \[ b_i+2c_i h_i+3 d_i h_i^2= b_{i+1},\ i=0,1,\ldots,n-1. \] La condición para \(i=n-1\) es consecuencia de la definición de \(b_n\) o de afirmar que \(S'_{n-1}(x_n)=b_n\).

Definimos \(c_n=\frac{S''(x_n)}{2}=\frac{S''_{n-1}(x_n)}{2}.\)

Los valores de \(S_i''(x)\) son: \[ S_i''(x)=2c_i+6d_i(x-x_i). \] Por tanto \(S_i''(x_i)=2c_i\), para \(i=0,1,\ldots,n-1\).

Si expresamos la condición \(S_{i}''(x_{i+1})=S_{i+1}''(x_{i+1})=2c_{i+1}\), tenemos que: \[ 2c_i +6 d_i h_i= 2c_{i+1},\ \Rightarrow c_i+3d_ih_i=c_{i+1},\ i=0,1,\ldots,n-1. \] La condición para \(i=n-1\) es consecuencia de la definición de \(c_n\) o de afirmar que \(\frac{S''_{n-1}(x_n)}{2}=c_n\).

Si despejamos los valores \(d_i\) de las ecuaciones anteriores obtenemos que: \[ d_i=\frac{c_{i+1}-c_i}{3h_i},\ i=0,1,\ldots,n-1. \] Sustituyendo estos valores en los conjuntos de ecuaciones anteriores para \(b_i\) y \(c_i\), obtenemos: \[ \begin{align*} a_i+b_i h_i+\frac{(c_{i+1}+2c_i)}{3} h_i^2= & a_{i+1},\\ b_i+(c_{i+1}+c_i) h_i= & b_{i+1}, \end{align*} \] para \(i=0,1,\ldots,n-1\).

Del primer conjunto de ecuaciones, podemos obtener los valores \(b_i\): \[ b_i=\frac{1}{h_i}(a_{i+1}-a_i)-\frac{h_i}{3}(c_{i+1}+2c_i),\ i=0,1,\ldots,n-1. \] Sustituyendo en el segundo sistema de ecuaciones obtenemos: \[ h_{i+1}c_{i+2}+2(h_{i+1}+h_i)c_{i+1}+h_i c_i=\frac{3}{h_{i+1}}(a_{i+2}-a_{i+1})-\frac{3}{h_i}(a_{i+1}-a_i), \] para \(i=0,1,\ldots,n-2\).

El último conjunto de ecuaciones es un sistema lineal con \(n-1\) ecuaciones y \(n+1\) incógnitas: \(c_0,c_1,\ldots,c_n\) ya que recordemos que los valores \(a_i=f(x_i)\) son conocidos.

Las dos incógnitas libres de las \(c_i\) quedarán fijadas usando las condiciones 5.1 (splines naturales) o 5.2 (splines con frontera fija).

Una vez hallados los coeficientes \(c_i\), podemos hallar los coeficientes \(b_i\) y \(d_i\) a partir de las relaciones anteriores.

Teorema de existencia y unicidad de splines naturales.

Sea \(f\) una función definida en un intervalo \([a,b]\) y sean \(a=x_0<x_1<\cdots <x_n=b\), \(n+1\) nodos en dicho intervalo.

Entonces existe un único spline natural interpolador \(S\) de \(f\) en los nodos anteriores. Es decir, \(S\) interpola los nodos anteriores, es de clase \({\cal C}^2[a,b]\) en \([a,b]\) y además, \(S''(a)=S''(b)=0\).

• INPUT n; x_0,x_1,...,x_n; a_0=f(x_0), a_1=f(x_1),...,a_n=f(x_n). (Damos los valores de \(n\), \(x_0,x_1,\ldots,x_n\) y los valores de \(a_0,a_1,\ldots,a_n\) con \(a_i=f(x_i)\), \(i=0,1,\ldots,n\))
• For i=0,...,n-1
  • Set h_i=x_{i+1}-x_i (Definimos los valores \(h_i\), \(i=0,1,\ldots,n-1\))
• For i=1,...,n-1
  • Set y_i=(3/h_i)*(a_{i+1}-a_i)-(3/(h_{i-1})*(a_i-a_{i-1})). (Definimos el vector de términos independientes \(\mathbf{y}\))
• Set A=0. (Vamos a definir la matriz del sistema \(\mathbf{A}\). Primero la inicializamos a \(0\) y seguidamente definiremos todas sus componentes)
• Set A[1,1]=2*(h_0+h_1).
• Set A[1,2]=h_1.

• For i=2,...,n-2
  • Set A[i,i-1]=h_{i-1};
  • Set A[i,i]=2*(h_{i-1}+h_i);
  • Set A[i,i+1]=h_i.
• Set A[n-1,n-2]=h_{n-2}.
• Set A[n-1,n-1]=2*(h_{n-2}+h_{n-1}).
• Resolvemos el sistema A*c=y usando un método de Álgebra lineal numérica directo como por ejemplo LU, o iterativo como Jacobi, Gauss-Seidel o SOR.
• Obtenemos c.

• For i=0,...,n-1 (Vamos a calcular los coeficientes \(b_i\) y \(d_i\), \(i=0,1,\ldots,n-1\))
  • Set b_i=(1/h_i)*(a_{i+1}-a_i)-(h_i/3)*(c_{i+1}+2*c_i);
  • Set d_i=(1/(3*h_i))*(c_{i+1}-c_i).
• Print a_i,b_i,c_i,d_i (Damos los coeficientes de los polinomios \(S_i(x)\), \(i=0,1,\ldots,n-1\))

Teorema de existencia y unicidad de splines de frontera fija

Sea \(f\) una función definida en un intervalo \([a,b]\) y sean \(a=x_0<x_1<\cdots <x_n=b\), \(n+1\) nodos en dicho intervalo.

Entonces existe un único spline de frontera fija interpolador \(S\) de \(f\) en los nodos anteriores. Es decir, \(S\) interpola los nodos anteriores, es de clase \({\cal C}^2[a,b]\) en \([a,b]\) y además, \(S'(a)=f'(a)\) y \(S'(b)=f'(b)\).

• INPUT n; x_0,x_1,...,x_n; a_0=f(x_0), a_1=f(x_1),...,a_n=f(x_n), fp0, fpn. (Damos los valores de \(n\), \(x_0,x_1,\ldots,x_n\) y los valores de \(a_0,a_1,\ldots,a_n\) con \(a_i=f(x_i)\), \(i=0,1,\ldots,n\), fp0\(=f'(a)\) y fpn\(=f'(b)\))
• For i=0,...,n-1
  • Set h_i=x_{i+1}-x_i. (Definimos los valores \(h_i\), \(i=0,1,\ldots,n-1\))
• Set y_0=3*(a_1-a_0)/h_0-3*fp0. (Vamos a definir el vector de términos independientes \(\mathbf{y}\))
• For i=1,...,n-1
  • Set y_i=(3/h_i)*(a_{i+1}-a_i)-(3/(h_{i-1})*(a_i-a_{i-1})).
• Set y_n=3*fpn-3*(a_n-a_{n-1})/h_{n-1}.

• Set A=0. (Vamos a definir la matriz del sistema \(\mathbf{A}\). Primero la inicializamos a \(0\) y seguidamente definiremos todas sus componentes)
• Set A[1,1]=2*h_0.
• Set A[1,2]=h_0.
• For i=2,...,n
  • Set A[i,i-1]=h_{i-2};
  • Set A[i,i]=2*(h_{i-2}+h_{i-1});
  • Set A[i,i+1]=h_{i-1}.
• Set A[n+1, n]=h_{n-1}.
• Set A[n+1,n+1]=2*h_{n-1}.

• Resolvemos el sistema A*c=y usando un método de Álgebra lineal numérica directo como por ejemplo LU, o iterativo como Jacobi, Gauss-Seidel o SOR.
• Obtenemos c.
• For i=0,...,n-1 (Vamos a calcular los coeficientes \(b_i\) y \(d_i\), \(i=0,1,\ldots,n-1\))
  • Set b_i=(1/h_i)*(a_{i+1}-a_i)-(h_i/3)*(c_{i+1}+2*c_i);
  • Set d_i=(1/(3*h_i))*(c_{i+1}-c_i).
• Print a_i,b_i,c_i,d_i (Damos los coeficientes de los polinomios \(S_i(x)\), \(i=0,1,\ldots,n-1\))

Desde el punto de vista visual, la aproximación usando splines de frontera fija es mejor que la aproximación usando splines naturales.

El siguiente resultado nos da una expresión del error cometido cuando usamos splines de frontera fija:

Teorema. Sea \(f\in {\cal C}^4[a,b]\) una función de clase \({\cal C}^4\) en el intervalo \([a,b]\). Sea \(M=\max_{x\in [a,b]}|f^{(4)}(x)|\). Si \(S\) es el único spline cúbico de frontera fija que interpola \(f\) en los nodos \(a=x_0<x_1<\cdots < x_n=b\). Entonces, para todo valor \(x\in [a,b]\),

\[ |f(x)-S(x)|\leq \frac{5M}{384}\max_{i=0,\ldots,n-1}(x_{i+1}-x_i)^4. \]