Las técnicas vistas en un curso clásico de cálculo diferencial e integral para hallar la derivada o la integral de una función, no son válidas en general cuando nos enfrentamos en un problema complicado de análisis numérico.

La mayoría piensa que una función \(f(x)\) es una expresión de la forma \(f(x)=\ldots\) donde \(\ldots\) es una expresión que contiene términos de funciones conocidas aplicadas a la variable \(x\). En estos casos, los problemas suelen ser fáciles de tratar ya que conocemos explícitamente la expresión de \(f(x)\).

Sin embargo, cuando lidiamos con un problema complejo de análisis numérico, la expresión de \(f(x)\) no es conocida. Podemos pensar que \(f(x)\) es un programa informático de un número determinado de líneas que tiene la variable \(x\) como input y nos da un valor \(f(x)\) como output.

En casos como los descritos anteriormente, no podemos hallar \(f'(x)\) ni \(\int f(x)\, dx\) usando las técnicas vistas en el curso de cálculo ya que dichas técnicas presuponen que conocemos la expresión explícita de \(f(x)\).

En estos casos, para hallar la derivada de \(f(x)\) en un valor determinado \(x\), \(f'(x)\) o para hallar la integral \(\int_a^b f(x)\, dx\) de \(f(x)\) entre dos valores concretos \(a\) y \(b\), necesitamos conocer otro tipo de técnicas que vamos a aprender en este capítulo.

Además, el problema es incluso más grave de lo que uno podría pensar:

• en primer lugar, el coste computacional de las técnicas de cálculo para hallar \(f'(x)\) o \(\int f(x)\, dx\) es elevado y
• en segundo lugar,conocer la expresión explícita de \(f(x)\) no garantiza que podamos hallar la derivada o la integral de \(f(x)\). No sólo eso, la mayoría de funciones que uno podría pensar no se pueden integrar usando técnicas de cálculo.

Como ejemplo, consideremos la función campana de Gauss \(f(x)=\frac{1}{\sqrt{2\pi}}\mathrm{e}^{-\frac{x^2}{2}}\). Para dicha función, no es posible hallar una expresión de una primitiva en términos de funciones conocidas y hay que integrarla usando técnicas numéricas.

Recordemos la definición de derivada de una función \(f(x)\) en un valor \(x_0\): \[ f'(x_0)=\lim_{x\to x_0}\frac{f(x)-f(x_0)}{x-x_0}=\lim_{h\to 0}\frac{f(x_0+h)-f(x_0)}{h}. \] Las dos expresiones anteriores son equivalentes, basta considerar \(h=x-x_0\) para pasar de la primera a la segunda.

Una manera sencilla de aproximar \(f'(x_0)\) sería considerar el cociente incremental \[ f'(x_0)\approx \frac{f(x_0+h)-f(x_0)}{h}, \] donde \(h\) sería un valor pequeño de cara a que \(x_0+h\) esté cerca de \(x_0\).

El problema de la fórmula anterior es que no tenemos ninguna expresión del error cometido.

De cara a solventar dicho problema, vamos a interpolar la función \(f(x)\) en un conjunto de puntos “cercanos” a \(x_0\) y derivar la expresión obtenida de cara a obtener una expresión del error.

Consideremos en primer lugar los puntos \(x_0\) y \(x_0+h\). Si \(f\in {\cal C}^2[a,b]\) es de clase \({\cal C}^2\) en un intervalo que contenga los puntos anteriores, podemos usar la fórmula de error de interpolación y escribir que: \[ f(x)=P_{0,1}(x)+\frac{(x-x_0)(x-(x_0+h))}{2}\cdot f''(\xi(x)), \] donde \(P_{0,1}(x)\) es el polinomio de interpolación en los puntos \((x_0,f(x_0))\) y \((x_0+h,f(x_0+h))\) y \(\xi(x)\in <x,x_0,x_0+h>\) (mínimo intervalo que contiene los puntos \(x\), \(x_0\) y \(x_0+h\))

El polinomio \(P_{0,1}(x)\) vale usando los polinomios de Lagrange: \[ P_{0,1}(x)=f(x_0)\cdot \frac{(x-x_0-h)}{(-h)}+f(x_0+h)\frac{(x-x_0)}{h}. \]

Por tanto, \[ \begin{align*} f(x)= & f(x_0)\cdot \frac{(x-x_0-h)}{(-h)}+f(x_0+h)\frac{(x-x_0)}{h}\\ & +\frac{(x-x_0)(x-(x_0+h))}{2}\cdot f''(\xi(x)). \end{align*} \]

Derivando la expresión anterior, obtenemos: \[ \begin{align*} f'(x)= &\frac{f(x_0+h)-f(x_0)}{h} +D_x\left[\frac{(x-x_0)(x-(x_0+h))}{2}\cdot f''(\xi(x))\right]\\ = & \frac{f(x_0+h)-f(x_0)}{h} + \frac{2(x-x_0)-h}{2}f''(\xi(x))\\ & + \frac{(x-x_0)(x-(x_0+h))}{2}\cdot D_x\left[f''(\xi(x))\right]. \end{align*} \] El problema de la expresión anterior es el término \(D_x\left[f''(\xi(x))\right]\) del que no sabemos calcular al no conocer el valor de \(\xi(x)\).

Sin embargo, como nos interesa \(f'(x_0)\), este término desaparece para \(x=x_0\): \[ f'(x_0)= \frac{f(x_0+h)-f(x_0)}{h} -\frac{h}{2}f''(\xi(x)). \]

Entonces podemos aproximar \(f'(x_0)\) por \(\displaystyle\frac{f(x_0+h)-f(x_0)}{h}\) con un error acotado por \(\frac{M|h|}{2}\) donde \[ M=\max_{x\in <x_0,x_0+h>}|f''(x)|. \]

Si \(h>0\), la fórmula anterior se conoce como fórmula de diferencias hacia adelante y si \(h<0\), fórmula de diferencias hacia atrás.

Vamos a generalizar el procedimiento anterior.

Sean \(x_0,x_1,\ldots, x_n\), \(n+1\) números distintos en algún intervalo \([a,b]\) y sea \(f\in {\cal C}^{n+1}[a,b]\) una función de clase \({\cal C}^{n+1}\) en dicho intervalo. Usando la fórmula del error en la interpolación podemos escribir: \[ f(x)=\sum_{k=0}^n f(x_k)L_k(x)+\frac{(x-x_0)\cdots (x-x_n)}{(n+1)!}f^{(n+1)}(\xi(x)), \] donde \(L_k(x)\) son los polinomios de Lagrange para \(k=0,1,\ldots,n\) y \(\xi(x)\in <x_0,\ldots,x_n,x>\).

Si derivamos la expresión anterior, obtenemos: \[ \begin{align*} f'(x)=& \sum_{k=0}^n f(x_k)L_k'(x)+D_x\left[\frac{(x-x_0)\cdots (x-x_n)}{(n+1)!}f^{(n+1)}(\xi(x))\right]\\ = & \sum_{k=0}^n f(x_k)L_k'(x)+D_x\left[\frac{(x-x_0)\cdots (x-x_n)}{(n+1)!}\right]f^{(n+1)}(\xi(x))\\ & + \frac{(x-x_0)\cdots (x-x_n)}{(n+1)!}\cdot D_x\left[f^{(n+1)}(\xi(x))\right]. \end{align*} \]

Si el punto donde aproximamos la derivada es uno de los nodos \(x_i\), nos “desaparece” el término que “más molesta” y nos queda: \[ \begin{align*} f'(x_i)=& \sum_{k=0}^n f(x_k)L_k'(x_i)\\ & +\frac{1}{(n+1)!}D_x\left[(x-x_0)\cdots (x-x_n)\right]_{|x=x_i}f^{(n+1)}(\xi(x)). \end{align*} \]

El valor de \(D_x\left[(x-x_0)\cdots (x-x_n)\right]|_{x=x_i}\) vale: \[ \begin{align*} D_x\left[\prod_{k=0}^n (x-x_k)\right]_{|x=x_i}= & D_x\left[(x-x_i)\prod_{k\neq i}(x-x_k)\right]_{|x=x_i}\\ = & \prod_{k\neq i}(x-x_k)_{|x=x_i}+(x_i-x_i)\cdot \prod_{k\neq i}(x-x_k)_{|x=x_i}\\ = & \prod_{k\neq i}(x_i-x_k) \end{align*} \]

La expresión de \(f'(x_i)\) queda de la forma siguiente: \[ f'(x_i)= \sum_{k=0}^n f(x_k)L_k'(x_i)+\frac{f^{(n+1)}(\xi_i)}{(n+1)!}\prod_{k\neq i}(x_i-x_k). \] A la expresión anterior se le conoce como fórmula de \(n+1\) puntos para aproximar \(f'(x_i)\).

Consideremos el caso particular en que \(n=2\) o tenemos tres puntos \(x_0\), \(x_1\) y \(x_2\).

Los polinomios de Lagrange y sus derivadas son los siguientes: \[ \begin{align*} L_0(x)=&\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)},\ \Rightarrow L_0'(x)=\frac{2x-x_1-x_2}{(x_0-x_1)(x_0-x_2)},\\ L_1(x)=&\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)},\ \Rightarrow L_1'(x)=\frac{2x-x_0-x_2}{(x_1-x_0)(x_1-x_2)},\\ L_2(x)=&\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)},\ \Rightarrow L_2'(x)=\frac{2x-x_0-x_1}{(x_2-x_0)(x_2-x_1)}. \end{align*} \]

La expresión de \(f'(x_i)\) para \(i=0,1,2\) es: \[ \begin{align*} f'(x_i)=&f(x_0)\left(\frac{2x_i-x_1-x_2}{(x_0-x_1)(x_0-x_2)}\right)+f(x_1)\left(\frac{2x_i-x_0-x_2}{(x_1-x_0)(x_1-x_2)}\right)\\ & +f(x_2)\left(\frac{2x_i-x_0-x_1}{(x_2-x_0)(x_2-x_1)}\right)+\frac{1}{6}f'''(\xi_i)\prod_{k=0,k\neq i}^2 (x_i-x_k). \end{align*} \]

Supongamos ahora que los tres puntos están equiespaciados, es decir, \(x_1=x_0+h\) y \(x_2=x_0+2h\) para un cierto \(h\).

Aplicando la fórmula anterior, tenemos las expresiones siguientes para \(f'(x_0)\), \(f'(x_1)\) y \(f'(x_2)\): \[ \begin{align*} f'(x_0)= & \frac{1}{h}\left(-\frac{3}{2}f(x_0)+2f(x_1)-\frac{1}{2}f(x_2)\right)+\frac{h^2}{3}f'''(\xi_0),\\ f'(x_1)= & \frac{1}{h}\left(-\frac{1}{2}f(x_0)+\frac{1}{2}f(x_2)\right)-\frac{h^2}{6}f'''(\xi_1),\\ f'(x_2)= & \frac{1}{h}\left(\frac{1}{2}f(x_0)-2 f(x_1)+\frac{3}{2}f(x_2)\right)+\frac{h^2}{3}f'''(\xi_2). \end{align*} \]

Aunque aparezcan tres fórmulas en realidad sólo tenemos dos ya que la primera y la última son la misma.

En la primera tenemos una aproximación de \(f'(x_0)\) usando los valores \(x_0\), \(x_0+h\) y \(x_0+2h\) y en la tercera, una aproximación de \(f'(x_0+2h)\) usando los mismos valores.

Si en la tercera “cambiamos” los papeles de \(x_0+2h\) por \(x_0\) y consideramos \(h<0\), nos sale la primera: \[ \begin{align*} \mbox{Primera fórmula} & \leftrightarrow \mbox{Tercera fórmula}\\ x_0 & \leftrightarrow x_0+2h\\ x_0+h & \leftrightarrow x_0+h\\ x_0+2h &\leftrightarrow x_0 \end{align*} \]

En resumen, tenemos las dos fórmulas siguientes de tres puntos para la aproximación de la derivada:

• Fórmula de los tres puntos respecto del punto medio: \[ f'(x_0)= \frac{1}{2h}(f(x_0+h)-f(x_0-h))-\frac{h^2}{6}f'''(\xi_1). \] Basta aplicar la segunda expresión anterior cambiando los papeles de \(x_0\), \(x_1=x_0+h\) y \(x_2=x_0+2h\) por \(x_0-h\), \(x_0\) y \(x_0+h\).

• Fórmula de los tres puntos respecto del punto extremo: \[ f'(x_0)=\frac{1}{2h}(-3 f(x_0)+4 f(x_0+h)-f(x_0+2h))+\frac{h^2}{3}f'''(\xi_0). \] Basta aplicar la primera expresión anterior.

Observación.

Siempre que sea posible, hay que aplicar la fórmula respecto del punto medio ya que el error queda reducido a la mitad.

Si en lugar de usar tres puntos equiespaciados, usamos cinco puntos equiespaciados de la forma \(x_0,x_0+h,x_0+2h,x_0+3h,x_0+4h\), y razonamos de la misma manera, es decir, calculamos el polinomio de interpolación en los puntos anteriores usando polinomios de Lagrange, usamos la fórmula de error de interpolación y la derivamos, obtenemos las siguientes fórmulas de cinco puntos:

• Fórmula de los cinco puntos respecto del punto medio: \[ \begin{align*} f'(x_0)=& \frac{1}{12h}(f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h))\\ & +\frac{h^4}{30}f^{(5)}(\xi), \end{align*} \] donde \(\xi\in (x_0-2h,x_0+2h)\).

• Formula de los cinco punto respecto del valor extremo: \[ \begin{align*} f'(x_0)=& \frac{1}{12h}(-25f(x_0)+48f(x_0+h)-36f(x_0+2h)\\ & +16f(x_0+3h) -3f(x_0+4h))+\frac{h^4}{5}f^{(5)}(\xi), \end{align*} \] donde \(\xi\in (x_0,x_0+4h)\).

Otra manera de deducir fórmulas de derivación numérica es usar la fórmula de Taylor.

Vamos a ver cómo se obtiene la fórmula de los cinco puntos respecto del punto medio usando desarrollos de Taylor.

Si desarrollamos por Taylor la función \(f\) alrededor de \(x_0\) en los puntos \(x_0-2h\), \(x_0-h\), \(x_0+h\) y \(x_0+2h\), obtenemos: \[ \begin{align*} f(x_0-2h)=& f(x_0)-f'(x_0)2h+f''(x_0) \frac{4h^2}{2}-f'''(x_0)\frac{8 h^3}{6}\\ & +f^{(4)}(x_0)\frac{16 h^4}{24}-f^{(5)}(\xi_1)\frac{32 h^5}{120},\\ f(x_0-h)=& f(x_0)-f'(x_0)h+f''(x_0) \frac{h^2}{2}-f'''(x_0)\frac{h^3}{6}\\ & +f^{(4)}(x_0)\frac{ h^4}{24}-f^{(5)}(\xi_2)\frac{h^5}{120}, \end{align*} \]

\[ \begin{align*} f(x_0+h)=& f(x_0)+f'(x_0)h+f''(x_0) \frac{h^2}{2}+f'''(x_0)\frac{h^3}{6}\\ & +f^{(4)}(x_0)\frac{ h^4}{24}+f^{(5)}(\xi_3)\frac{h^5}{120},\\ f(x_0+2h)=& f(x_0)+f'(x_0)2h+f''(x_0) \frac{4h^2}{2}+f'''(x_0)\frac{8 h^3}{6}\\ & +f^{(4)}(x_0)\frac{16 h^4}{24}+f^{(5)}(\xi_4)\frac{32 h^5}{120}, \end{align*} \]

donde \(\xi_1,\xi_2,\xi_3,\xi_4\in (x_0-2h,x_0+2h)\).

El siguiente paso es hallar unos coeficientes \(A_1\), \(A_2\), \(A_3\), \(A_4\) y \(A_5\) por los que hay que multiplicar las expresiones de \(f(x_0-2h), f(x_0-h), f(x_0), f(x_0+h)\) y \(f(x_0+2h)\), respectivamente, de cara a eliminar los términos en \(h^0\), \(h^2\), \(h^3\) y \(h^4\) de cara a que quede el término en \(h f'(x_0)\) para obtener una fórmula aproximada de \(f'(x_0)\):

\[ \begin{align*} \color{red}{A_1}\cdot \Biggl( f(x_0-2h) = & f(x_0)-f'(x_0)2h+f''(x_0) \frac{4h^2}{2}-\cdots,\Biggr)\\ \color{red}{A_2}\cdot \Biggl(f(x_0-h)=& f(x_0)-f'(x_0)h+f''(x_0) \frac{h^2}{2}-\cdots\Biggr),\\ \color{red}{A_3}\cdot (f(x_0)= & f(x_0)),\\ \color{red}{A_4}\cdot \Biggl(f(x_0+h)=& f(x_0)+f'(x_0)h+f''(x_0) \frac{h^2}{2}+\cdots \Biggr),\\ \color{red}{A_5}\cdot \Biggl( f(x_0+2h) = & f(x_0)+f'(x_0)2h+f''(x_0) \frac{4h^2}{2}+\cdots,\Biggr)\\ \end{align*} \]

• El coeficiente correspondiente a \(f(x_0)\) es: \(A_1+A_2+A_3+A_4+A_5\).
• El coeficiente correspondiente a \(f''(x_0)h^2\) es: \(2A_1+\frac{1}{2}A_2+\frac{1}{2}A_4+2A_5\).
• El coeficiente correspondiente a \(f'''(x_0)h^3\) es: \(-\frac{4}{3} A_1-\frac{1}{6}A_2+\frac{1}{6}A_4+\frac{4}{3}A_5\).
• El coeficiente correspondiente a \(f^{(4)}(x_0)h^4\) es: \(\frac{2}{3} A_1+\frac{1}{24}A_2+\frac{1}{24}A_4+\frac{2}{3}A_5\).

Entonces de cara a anular los coeficientes anteriores, tenemos que resolver el sistema siguiente: \[ \left. \begin{align*} A_1+A_2+A_3+A_4+A_5= & 0,\\ 2A_1+\frac{1}{2}A_2+\frac{1}{2}A_4+2A_5 =& 0,\\ -\frac{4}{3} A_1-\frac{1}{6}A_2+\frac{1}{6}A_4+\frac{4}{3}A_5 = & 0,\\ \frac{2}{3} A_1+\frac{1}{24}A_2+\frac{1}{24}A_4+\frac{2}{3}A_5 = &0. \end{align*} \right\} \]

Es un sistema compatible indeterminado ya que tiene 4 ecuaciones con 5 incógnitas. Podemos añadir una ecuación extra imponiendo que el coeficiente correspondiente a \(hf'(x_0)\) sea \(1\): \[ -2A_1-A_2+A_4+2A_5=1. \]

Las soluciones del sistema anterior son las siguientes: \[ A_1=\frac{1}{12},\ A_2=-\frac{2}{3},\ A_3=0,\ A_4=\frac{2}{3},\ A_5=-\frac{1}{12}. \] Entonces \(A_1\cdot f(x_0-2h)+A_2 f(x_0-h)+A_3 f(x_0)+A_4 f(x_0+h)+A_5 f(x_0+2h)\) vale: \[ \begin{align*} & \frac{1}{12} f(x_0-2h)-\frac{2}{3}f(x_0-h)+\frac{2}{3}f(x_0+h)-\frac{1}{12} f(x_0-2h)\\ & =f'(x_0)h+\frac{h^5}{120}\left(-\frac{32}{12}f^{(5)}(\xi_1)+\frac{2}{3}f^{(5)}(\xi_2)+\frac{2}{3}f^{(5)}(\xi_3)-\frac{32}{12}f^{(5)}(\xi_4)\right). \end{align*} \]

Aplicando el Teorema de Bolzano generalizado podemos reducir el término del error de la forma siguiente:

• Existe un \(\xi_{1,4}\) tal que \(-\frac{32}{12}f^{(5)}(\xi_1)-\frac{32}{12}f^{(5)}(\xi_4)=-\frac{16}{3}f^{(5)}(\xi_{1,4})\).
• Existe un \(\xi_{2,3}\) tal que \(\frac{2}{3}f^{(5)}(\xi_2)+\frac{2}{3}f^{(5)}(\xi_3)=\frac{4}{3}f^{(5)}(\xi_{2,3})\).

El término del error queda pues: \[ \frac{h^5}{120}\left(-\frac{16}{3}f^{(5)}(\xi_{1,4})+\frac{4}{3}f^{(5)}(\xi_{2,3})\right). \] Si hubiésemos usado la fórmula del error de interpolación y los polinomios de Lagrange, hubiésemos visto que el término del error se podría escribir como: \[ -\frac{16}{3}f^{(5)}(\xi_{1,4})+\frac{4}{3}f^{(5)}(\xi_{2,3}) =\left(-\frac{16}{3}+\frac{4}{3}\right)f^{(5)}(\xi)=-4f^{(5)}(\xi). \]

En resumen, \[ \begin{align*} & \frac{1}{12} f(x_0-2h)-\frac{2}{3}f(x_0-h)+\frac{2}{3}f(x_0+h)-\frac{1}{12} f(x_0+2h)\\ & =f'(x_0)h-\frac{h^5}{120}4f^{(5)}(\xi), \end{align*} \] de donde despejando \(f'(x_0)\) obtenemos la fórmula de derivación de los cinco puntos: \[ \begin{align*} f'(x_0)=& \frac{1}{h}\left(\frac{1}{12} f(x_0-2h)-\frac{2}{3}f(x_0-h)+\frac{2}{3}f(x_0+h)-\frac{1}{12} f(x_0+2h)\right)\\ & +\frac{h^4}{30}f^{(5)}(\xi),\\ = & \frac{f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h)}{12h}+\frac{h^4}{30}f^{(5)}(\xi) \end{align*} \]

La técnica anterior nos permite usar fórmulas numéricas aproximadas para hallar derivadas de orden superior.

Calculemos como ejemplo una aproximación de \(f''(x_0)\).

Considerando los mismos puntos anteriores, \(x_0\pm h\) y \(x_0\pm 2h\) y los mismos desarrollos de Taylor anteriores, hemos de calcular unos coeficientes \(B_1,B_2,B_3,B_4\) y \(B_5\) por los que hay que multiplicar las expresiones de \(f(x_0-2h), f(x_0-h), f(x_0), f(x_0+h)\) y \(f(x_0+2h)\), respectivamente, de cara a eliminar los términos en \(h^0\), \(h\), \(h^3\) y \(h^4\) de cara a que quede el término en \(h^2 f''(x_0)\) para obtener una fórmula aproximada de \(f''(x_0)\).

En primer lugar, desarrollamos por Taylor la función \(f\) alrededor de \(x_0\) en los puntos \(x_0-2h\), \(x_0-h\), \(x_0+h\) y \(x_0+2h\) hasta llegar a términos del orden \(h^6\): \[ \begin{align*} f(x_0-2h)=& f(x_0)-f'(x_0)2h+f''(x_0) \frac{4h^2}{2}-f'''(x_0)\frac{8 h^3}{6}\\ & +f^{(4)}(x_0)\frac{16 h^4}{24}-f^{(5)}(x_0)\frac{32 h^5}{120}+f^{(6)}(\xi_1) \frac{64 h^6}{720},\\ f(x_0-h)=& f(x_0)-f'(x_0)h+f''(x_0) \frac{h^2}{2}-f'''(x_0)\frac{h^3}{6}\\ & +f^{(4)}(x_0)\frac{ h^4}{24}-f^{(5)}(x_0)\frac{h^5}{120}+f^{(6)}(\xi_2) \frac{h^6}{720}, \end{align*} \]

\[ \begin{align*} f(x_0+h)=& f(x_0)+f'(x_0)h+f''(x_0) \frac{h^2}{2}+f'''(x_0)\frac{h^3}{6}\\ & +f^{(4)}(x_0)\frac{ h^4}{24}+f^{(5)}(x_0)\frac{h^5}{120}+f^{(6)}(\xi_3) \frac{h^6}{720},\\ f(x_0+2h)=& f(x_0)+f'(x_0)2h+f''(x_0) \frac{4h^2}{2}+f'''(x_0)\frac{8 h^3}{6}\\ & +f^{(4)}(x_0)\frac{16 h^4}{24}+f^{(5)}(x_0)\frac{32 h^5}{120}+f^{(6)}(\xi_4) \frac{64 h^6}{720}, \end{align*} \]

\[ \begin{align*} \color{red}{B_1}\cdot \Biggl( f(x_0-2h) = & f(x_0)-f'(x_0)2h+f''(x_0) \frac{4h^2}{2}-\cdots,\Biggr)\\ \color{red}{B_2}\cdot \Biggl(f(x_0-h)=& f(x_0)-f'(x_0)h+f''(x_0) \frac{h^2}{2}-\cdots\Biggr),\\ \color{red}{B_3}\cdot (f(x_0)= & f(x_0)),\\ \color{red}{B_4}\cdot \Biggl(f(x_0+h)=& f(x_0)+f'(x_0)h+f''(x_0) \frac{h^2}{2}+\cdots \Biggr),\\ \color{red}{B_5}\cdot \Biggl( f(x_0+2h) = & f(x_0)+f'(x_0)2h+f''(x_0) \frac{4h^2}{2}+\cdots,\Biggr)\\ \end{align*} \]

• El coeficiente correspondiente a \(f(x_0)\) es: \(B_1+B_2+B_3+B_4+B_5\).
• El coeficiente correspondiente a \(f'(x_0)h\) es: \(-2B_1-B_2+B_4+2B_5\).
• El coeficiente correspondiente a \(f'''(x_0)h^3\) es: \(-\frac{4}{3} B_1-\frac{1}{6}B_2+\frac{1}{6}B_4+\frac{4}{3}B_5\).
• El coeficiente correspondiente a \(f^{(4)}(x_0)h^4\) es: \(\frac{2}{3} B_1+\frac{1}{24}B_2+\frac{1}{24}B_4+\frac{2}{3}B_5\).

Entonces de cara a anular los coeficientes anteriores, tenemos que resolver el sistema siguiente: \[ \left. \begin{align*} B_1+B_2+B_3+B_4+B_5= & 0,\\ -2B_1-B_2+B_4+2B_5 =& 0,\\ -\frac{4}{3} B_1-\frac{1}{6}B_2+\frac{1}{6}B_4+\frac{4}{3}B_5 = & 0,\\ \frac{2}{3} B_1+\frac{1}{24}B_2+\frac{1}{24}B_4+\frac{2}{3}B_5 = &0. \end{align*} \right\} \]

Es un sistema compatible indeterminado ya que tiene 4 ecuaciones con 5 incógnitas. Podemos añadir una ecuación extra imponiendo que el coeficiente correspondiente a \(h^2f''(x_0)\) sea \(1\): \[ 2B_1+\frac{1}{2}B_2+\frac{1}{2}B_4+2B_5=1. \]

Las soluciones del sistema anterior son las siguientes: \[ B_1=-\frac{1}{12},\ B_2=\frac{4}{3},\ B_3=-\frac{5}{2},\ B_4=\frac{4}{3},\ B_5=-\frac{1}{12}. \] En este caso, observamos que el coeficiente de \(h^5\) también se anula: \[ -\frac{32B_1}{120}-\frac{B_2}{120}+\frac{B_4}{120}+\frac{32B_5}{120}=\frac{1}{1440}(32-48+48-32)=0. \]

Entonces \(B_1\cdot f(x_0-2h)+B_2 f(x_0-h)+B_3 f(x_0)+B_4 f(x_0+h)+B_5 f(x_0+2h)\) vale: \[ \begin{align*} & -\frac{1}{12} f(x_0-2h)+\frac{4}{3}f(x_0-h)-\frac{5}{2}f(x_0)+\frac{4}{3}f(x_0+h)-\frac{1}{12} f(x_0-2h)\\ & =f''(x_0)h^2+\frac{h^6}{720}\left(-\frac{64}{12}f^{(6)}(\xi_1)+\frac{4}{3}f^{(6)}(\xi_2)+\frac{4}{3}f^{(6)}(\xi_3)-\frac{64}{12}f^{(6)}(\xi_4)\right). \end{align*} \]

Aplicando el Teorema de Bolzano generalizado podemos reducir el término del error de la forma siguiente:

• Existe un \(\xi_{1,4}\) tal que \(\frac{64}{12}f^{(6)}(\xi_1)+\frac{64}{12}f^{(6)}(\xi_4)=\frac{32}{3} f^{(6)}(\xi_{1,4})\).
• Existe un \(\xi_{2,3}\) tal que \(\frac{4}{3}f^{(6)}(\xi_2)+\frac{4}{3}f^{(6)}(\xi_3)=\frac{8}{3} f^{(6)}(\xi_{2,3})\).

El término del error queda pues: \[ \frac{h^6}{720}\left(\frac{8}{3}f^{(6)}(\xi_{2,3})-\frac{32}{3}f^{(6)}(\xi_{1,4})\right)=\frac{h^6}{270}\left(f^{(6)}(\xi_{2,3})-4f^{(6)}(\xi_{1,4})\right) \]

En resumen, \[ \begin{align*} & -\frac{1}{12} f(x_0-2h)+\frac{4}{3}f(x_0-h)-\frac{5}{2}f(x_0)+\frac{4}{3}f(x_0+h)-\frac{1}{12} f(x_0+2h)\\ & =f''(x_0)h^2-\frac{h^6}{270}\left(4f^{(6)}(\xi_{1,4})-f^{(6)}(\xi_{2,3})\right), \end{align*} \]

de donde despejando \(f''(x_0)\) obtenemos la fórmula de derivación de los cinco puntos: \[ \begin{align*} f''(x_0)=& \frac{1}{h^2}\Biggl( -\frac{1}{12} f(x_0-2h)+\frac{4}{3}f(x_0-h)-\frac{5}{2}f(x_0)+\frac{4}{3}f(x_0+h)\\ &-\frac{1}{12} f(x_0+2h)\Biggr) +\frac{h^4}{270}\left(4f^{(6)}(\xi_{1,4})-f^{(6)}(\xi_{2,3})\right),\\ = & \frac{1}{12h^2}\Bigl(-f(x_0-2h)+16f(x_0-h)-30 f(x_0)+16f(x_0+h)\\ & -f(x_0+2h)\Bigr)+\frac{h^4}{270}\left(4f^{(6)}(\xi_{1,4})-f^{(6)}(\xi_{2,3})\right). \end{align*} \]

En todas las aproximaciones anteriores, sólo hemos tenido en cuenta el error de truncamiento o el error que cometemos cuando aproximamos \(f'(x_0)\), \(f''(x_0)\) o incluso derivadas superiores a partir de interpolar o de desarrollar \(f\) usando la fórmula de Taylor pero no hemos tenido en cuenta el error de redondeo.

Es decir, no hemos tenido en cuenta que cuando evaluamos la función \(f\) en un cierto punto \(x\), en realidad no obtenemos el valor exacto \(f(x)\), sino una aproximación \(\tilde{f}(x)\).

Por ejemplo, en la fórmula de los tres puntos respecto del punto medio, para aproximar \(f'(x_0)\): \[ f'(x_0)= \frac{1}{2h}(f(x_0+h)-f(x_0-h))-\frac{h^2}{6}f'''(\xi_1), \] donde \(\xi_1\in <x_0-h,x_0+h>\), en realidad, cuando evaluamos la aproximación, \(\frac{1}{2h}(f(x_0+h)-f(x_0-h))\), hacemos lo siguiente: \[ \frac{1}{2h}(\tilde{f}(x_0+h)-\tilde{f}(x_0-h)), \] donde \(\tilde{f}(x_0+h)\), \(\tilde{f}(x_0-h)\) son aproximaciones de \(f(x_0+h)\) y \(f(x_0-h)\), respectivamente.

Por tanto, podemos escribir, \[ \tilde{f}(x_0+h)=f(x_0+h)+e(x_0+h),\quad \tilde{f}(x_0-h)=f(x_0-h)+e(x_0-h), \] donde \(e(x_0+h)\) y \(e(x_0-h)\) son los errores de redondeo que cometemos al evaluar \(f(x_0+h)\) y \(f(x_0-h)\), respectivamente.

El error global cometido en la fórmula de los tres puntos respecto del punto medio seria el siguiente: \[ \begin{align*} & f'(x_0)-\frac{1}{2h}(\tilde{f}(x_0+h)-\tilde{f}(x_0-h))=\frac{1}{2h}(e(x_0+h)+e(x_0-h))\\ & -\frac{h^2}{6}f'''(\xi_1) \end{align*} \] Sea \(\epsilon >0\) el error máximo cometido al evaluar \(f(x_0\pm h)\) con \(|h|\leq h_0\) para un cierto \(h_0\).

Sea \(M=\max_{x\in (x_0-h_0,x_0+h_0)}|f'''(x)|\).

Podemos acotar el error global por: \[ E=\frac{1}{2h}\cdot 2\epsilon +\frac{h^2}{6}\cdot M=\frac{\epsilon}{h}+\frac{M h^2}{6}. \]

A partir de la expresión anterior, vemos que no tiene sentido considerar \(h\) arbitrariamente pequeños ya que, aunque el término \(\frac{M h^2}{6}\) tiende a cero, el término \(\frac{\epsilon}{h}\) se hace arbitrariamente grande.

Si tuviésemos conocimiento del comportamiento de la tercera derivada, podríamos calcular el \(h\) óptimo para el cual el error sería mínimo pero como trabajamos con cotas, dicho \(h\) óptimo sería aproximado y nunca lo podríamos hallar.

Vamos a aprender una técnica que nos permite transformar aproximaciones numéricas de algoritmos con errores de truncamiento relativamente bajos, es decir, de orden \(O(h^k)\) con \(k\) pequeños en algoritmos con errores de truncamiento altos, es decir, de orden \(O(h^\hat{k})\) con \(\hat{k}>k\).

Dicha técnica se denomina extrapolación de Richardson. Veamos cómo funciona.

Sea \(F(h)\) un algoritmo que nos da una aproximación numérica de cierta constante \(C\) con error de truncamiento \(h\). Es decir: \[ F(h)=C+k_1 h+k_2 h^2+k_3 h^3+\cdots, \] donde \(k_1,k_2,k_3,\ldots\) son constantes en principio desconocidas donde sí sabemos que \(k_1\neq 0\) ya que como hemos indicado el algoritmo tiene error de truncamiento \(O(h)\).

Si aplicamos el algoritmo a \(\frac{h}{2}\), obtenemos: \[ F\left(\frac{h}{2}\right)=C+k_1\cdot\frac{h}{2}+k_2\cdot\frac{h^2}{4}+k_3\cdot\frac{h^3}{8}+\cdots, \] Consideremos ahora \(2F\left(\frac{h}{2}\right)-F(h)\) como nueva aproximación de la constante \(C\), \[ 2F\left(\frac{h}{2}\right)-F(h)=C-k_2\cdot\frac{h^2}{2}-k_3\cdot \frac{3h^3}{4}+\cdots \] Vemos que el “nuevo” algoritmo \(2F\left(\frac{h}{2}\right)-F(h)\) tiene error de truncamiento \(O(h^2)\) y por tanto es mejor que el algoritmo inicial \(F(h)\) de cara a aproximar la constante \(C\).

Sea \[ \begin{align*} F_0(h)= &F(h)=C+k_1\cdot h+\cdots,\\ F_1(h)= & 2F\left(\frac{h}{2}\right)-F(h)=C-k_2\cdot \frac{h^2}{2}+\cdots \end{align*} \] En general, escribimos \(F_n(h)=C+k_{n+1}\frac{p_n}{q_n} h^{n+1}+\cdots\), donde \(k_{n+1}\frac{p_n}{q_n}\) representaría el coeficiente de término principal.

Para calcular \(F_{n+1}(h)\), tenemos que eliminar el término \(k_{n+1}\frac{p_n}{q_n} h^{n+1}\). Por tanto, hacemos lo siguiente: \[ 2^{n+1}F_n\left(\frac{h}{2}\right)-F_n(h)=(2^{n+1}-1)C+k_{n+2}\frac{\tilde{p}_{n+1}}{\tilde{q}_{n+1}}h^{n+2}+\cdots \]

Por tanto, \[ F_{n+1}(h)=\frac{2^{n+1}F_n\left(\frac{h}{2}\right)-F_n(h)}{2^{n+1}-1}=C+k_{n+2}\frac{p_{n+1}}{q_{n+1}}h^{n+2}+\cdots \]

Para calcular aproximaciones de \(C\) con errores de truncamiento cada vez más pequeños, hacemos una tabla como la que sigue:

Supongamos que la aproximación \(F(h)\) para calcular \(C\) tiene una expresión donde sólo aparecen exponentes pares en el término del error: \[ F_0(h)=F(h)=C+k_1 h^2+k_2 h^4+\cdots +h_n k^{2n}+\cdots \] En este caso, las aproximaciones sucesivas usando la extrapolación de Richardson son las siguientes: \[ F_{n+1}(h)=\frac{4^{n+1}F_n\left(\frac{h}{2}\right)-F_n(h)}{4^{n+1}-1}=C+k_{n+2}\frac{p_{n+1}}{q_{n+1}}h^{2(n+2)}+\cdots \]

La tabla de las aproximaciones sería en este caso:

En esta sección vamos a dar fórmulas aproximadas para calcular una integral definida \(\int_a^b f(x)\,dx\), donde \(f(x)\) es una función que pertenece a \({\cal C}^{k}[a,b]\), es decir de clase \({\cal C}^k\) para un cierto \(k\) y un intervalo \([a,b]\).

Dichas fórmulas se denominan fórmulas de cuadratura numérica. Dichas fórmulas son expresiones del tipo: \[ \int_a^b f(x)\, dx\approx\sum_{i=0}^n a_i\cdot f(x_i), \] donde \(a_i\) son unos coeficientes a determinar y \(x_i\) son unos nodos a determinar normalmente en el intervalo \([a,b]\).

Las fórmulas de cuadratura numérica están basadas en el polinomio de interpolación en un conjunto de nodos \(x_0,x_1,\ldots,x_n\in [a,b]\).

Por tanto, si escribimos el polinomio de interpolación en función de los polinomios de Lagrange y usando la fórmula del error en la interpolación, podemos escribir: \[ \begin{align*} & \int_a^b f(x)\, dx= \int_a^b \sum_{i=0}^n f(x_i)L_i(x)\, dx+\int_a^b\prod_{i=0}^n (x-x_i)\frac{f^{(n+1)}(\xi(x))}{(n+1)!}\, dx \\ & = \sum_{i=0}^n f(x_i)\int_a^b L_i(x)\,dx+\frac{1}{(n+1)!}\int_a^b \prod_{i=0}^n (x-x_i)f^{(n+1)}(\xi(x))\,dx. \end{align*} \]

Los coeficientes \(a_i\) serían: \[ a_i=\int_a^b L_i(x)\,dx, \] y el error cometido sería: \[ E(f,a,b)=\frac{1}{(n+1)!}\int_a^b \prod_{i=0}^n (x-x_i)f^{(n+1)}(\xi(x))\,dx. \]

Consideremos el caso más simple: \(n=1\), \(x_0=a\) y \(x_1=b\), en este caso aproximamos la integral \(\int_a^b f(x)\, dx\) por el área del trapecio de vértices \((a,0)\), \((a,f(a))\), \((b,0)\) y \((b,f(b))\):

En este caso, el polinomio de interpolación en los nodos \((x_0,f(x_0))\) y \((x_1,f(x_1))\) vale: \[ P_1(x)=f(x_0)\frac{x-x_1}{(x_0-x_1)}+f(x_1)\frac{x-x_0}{x_1-x_0}. \] La cuadratura sería:

\[ \begin{align*} \int_a^b f(x)\, dx= & f(x_0)\int_{x_0}^{x_1}\frac{x-x_1}{(x_0-x_1)}\, dx+f(x_1)\int_{x_0}^{x_1}\frac{x-x_0}{x_1-x_0}\, dx\\ & +\frac{1}{2}\int_{x_0}^{x_1}f''(\xi(x))(x-x_0)(x-x_1)\, dx\\ = & f(x_0)\left[\frac{(x-x_1)^2}{2(x_0-x_1)}\right]_{x_0}^{x_1}+f(x_1)\left[\frac{(x-x_0)^2}{2(x_1-x_0)}\right]_{x_0}^{x_1}+E(f,a,b)\\ = & -\frac{1}{2}f(x_0)(x_0-x_1)+\frac{1}{2}f(x_1)(x_1-x_0)+E(f,a,b)\\ = & \frac{1}{2}(x_1-x_0)(f(x_0)+f(x_1))+E(f,a,b)\\ = & \frac{1}{2}h(f(x_0)+f(x_1))+E(f,a,b), \end{align*} \] donde \(h=x_1-x_0\).

Analicemos el término del error: \[ E(f,a,b)=\frac{1}{2}\int_{x_0}^{x_1}f''(\xi(x))(x-x_0)(x-x_1)\, dx. \] Teniendo en cuenta que \((x-x_0)(x-x_1)\) no cambia de signo en el intervalo \([x_0,x_1]=[a,b]\) y usando el Teorema del valor medio para integrales, podemos escribir que existe un \(\xi\in (a,b)\) tal que:

\[ \begin{align*} E(f,a,b)=&\frac{f''(\xi)}{2}\int_{x_0}^{x_1}(x-x_0)(x-x_1)\, dx\\ = & \frac{f''(\xi)}{2}\int_{x_0}^{x_1}(x-x_0)(x-x_0-h)\, dx\\ = & \frac{f''(\xi)}{2}\int_{x_0}^{x_1}(x-x_0)^2-h (x-x_0)\, dx\\ = & \frac{f''(\xi)}{2}\left[\frac{(x-x_0)^3}{3}-h\frac{(x-x_0)^2}{2}\right]_{x_0}^{x_1}=\frac{f''(\xi)}{2}\left(\frac{h^3}{3}-\frac{h^3}{2}\right)\\ = & -\frac{h^3}{12}f''(\xi). \end{align*} \]

En resumen: \[ \int_a^b f(x)\, dx=\frac{h}{2}(f(a)+f(b))-\frac{h^3}{12}f''(\xi), \] donde \(h=b-a\) y \(\xi\in (a,b)\).

Consideremos ahora \(n=2\), \(x_0=a\), \(x_1=\frac{a+b}{2}\) y \(x_2=b\).

Si definimos \(h\) como \(h=\frac{b-a}{2}\), podemos escribir los nodos como \(x_0=a\), \(x_1=a+h\) y \(x_2=b\).

En este caso aproximamos la integral \(\int_a^b f(x)\, dx\) por el área de la parábola que pasa por los puntos \((a,f(a))\), \((x_1,f(x_1))\) y \((b,f(b))\):

Para deducir la regla de Simpson, en lugar de interpolar, vamos a considerar el desarrollo por Taylor de la función \(f\) alrededor de \(x=x_1\) ya que de esta forma se consigue una expresión del error de orden menor: \[ \begin{align*} f(x)=& f(x_1)+f'(x_1)(x-x_1)+\frac{1}{2}f''(x_1)(x-x_1)^2+\frac{1}{6}f'''(x_1)(x-x_1)^3\\ & +\frac{1}{24}f^{(4)}(\xi(x))(x-x_1)^4. \end{align*} \]

La cuadratura sería: \[ \begin{align*} \int_a^b f(x)\, dx= & f(x_1)(b-a)+f'(x_1)\int_a^b (x-x_1)\, dx\\ & +\frac{1}{2}f''(x_1)\int_a^b (x-x_1)^2\, dx+\frac{1}{6}f'''(x_1)\int_a^b (x-x_1)^3\, dx\\ & +\frac{1}{24} \int_a^b f^{(4)}(\xi(x))(x-x_1)^4\, dx \end{align*} \]

La cuadratura sería: \[ \begin{align*} = & f(x_1)(b-a)+f'(x_1)\left[\frac{(x-x_1)^2}{2}\right]_{a}^b\\ & + \frac{1}{2}f''(x_1)\left[\frac{(x-x_1)^3}{3}\right]_a^b +\frac{1}{6}f'''(x_1)\left[\frac{(x-x_1)^4}{4}\right]_a^b+E(f,a,b)\\ = & 2hf(x_1)+\frac{h^3}{3}f''(x_1)+E(f,a,b) \end{align*} \]

Teniendo en cuenta que \(x_0=x_1-h\) y \(x_2=x_1+h\), podemos usar la fórmula del ejercicio en la expresión de la cuadratura obteniendo: \[ \begin{align*} \int_a^b f(x)\, dx = & 2hf(x_1)\\ & +\frac{h^3}{3}\left(\frac{f(x_1+h)-2f(x_1)+f(x_1-h)}{h^2}-\frac{h^2}{12}f^{(4)}(\xi_1)\right) \\ & +E(f,a,b) \\ = & \frac{h}{3} (f(x_1+h)+4 f(x_1)+f(x_1-h))-\frac{h^5}{36}f^{(4)}(\xi_1)\\ & +E(f,a,b)\\ = & \frac{h}{3}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right)-\frac{h^5}{36}f^{(4)}(\xi_1)\\ & +E(f,a,b) \end{align*} \]

Analicemos el error cometido: \[ E(f,a,b)=\frac{1}{24} \int_a^b f^{(4)}(\xi(x))(x-x_1)^4\, dx. \]

Teniendo en cuenta que \((x-x_1)^4\) no cambia de signo en el intervalo \([a,b]\) y usando el Teorema del valor medio para integrales, podemos escribir que existe un \(\xi_2\in (a,b)\) tal que: \[ \begin{align*} E(f,a,b)=&\frac{f^{(4)}(\xi_2)}{24}\int_a^b (x-x_1)^4\, dx=\frac{f^{(4)}(\xi_2)}{24}\left[\frac{(x-x_1)^5}{5}\right]_a^b \\ = & \frac{f^{(4)}(\xi_2)}{60}h^5. \end{align*} \]

En resumen: \[ \begin{align*} \int_a^b f(x)\, dx= & \frac{h}{3} (f(x_1+h)+4 f(x_1)+f(x_1-h))\\ & +\frac{h^5}{12}\left(\frac{f^{(4)}(\xi_2)}{5}-\frac{f^{(4)}(\xi_1)}{3}\right), \end{align*} \] donde \(h=\frac{b-a}{2}\).

Intentando deducir la regla de Simpson como aquella regla para la cual la integral es exacta para polinomios de grado menor o igual que \(3\), puede deducirse que existe un valor \(\xi\) tal que: \[ \frac{f^{(4)}(\xi_2)}{5}-\frac{f^{(4)}(\xi_1)}{3}=-\frac{2 f^{(4)}(\xi)}{15}. \]

La regla de Simpson sería, pues: \[ \int_a^b f(x)\, dx= \frac{h}{3} (f(x_1+h)+4 f(x_1)+f(x_1-h))-\frac{h^5}{90}f^{(4)}(\xi). \]

Observación

Recordemos que el error en el método de los trapecios era de orden \(O(h^3)\).

Por tanto, en la regla de Simpson el orden del error esperable era \(O(h^4)\) ya que añadimos un nodo más en la interpolación pero el error nos ha salido de orden mucho menor, \(O(h^5)\). Por dicha razón, es uno de los métodos de integración numérica más usados.

Definición de precisión o grado de exactitud de una cuadratura.

Diremos que una cuadratura tiene precisión o grado de exactitud \(k\) si la fórmula es exacta para todos los polinomios de grado menor o igual que \(k\).

Observación.

Las fórmulas de cuadratura son lineales para las funciones, es decir, dadas dos funciones \(f\) y \(g\) y dados dos números \(\lambda\) y \(\mu\), entonces si llamamos \(\displaystyle C(f)=\sum_{i=0}^n a_i f(x_i)\) a la cuadratura para calcular \(\displaystyle\int_a^b f(x)\, dx\), se cumple: \[ C(\lambda f+\mu g)=\lambda\cdot C(f)+\mu\cdot C(g). \]

Usando la observación anterior, para comprobar que una cuadratura tiene precisión \(k\), basta ver que es exacta para los monomios \(x^i\), para \(i=0,1,\ldots,k\).

Proposición. Precisión de la fórmula de los trapecios y de la fórmula de Simpson.

La fórmula de los trapecios tiene precisión \(k=1\) y la fórmula de Simpson, precisión \(k=3\).

Vamos a generalizar el método de los trapecios o la fórmula de Simpson suponiendo que en lugar de considerar \(2\) puntos (trapecios) o \(3\) (Simpson), consideramos \(n+1\) puntos equiespaciados en el intervalo \([a,b]\): \[ x_i=a+i\cdot h, \ i=0,\ldots,n,\ h=\frac{b-a}{n}, \] e interpolamos para hallar la cuadratura correspondiente: \[ \int_a^b f(x)\, dx\approx \sum_{i=0}^n a_i f(x_i). \]

Dichas cuadraturas se llaman fórmulas de Newton-Cotes cerradas porque incluyen como nodos los extremos \(a\) y \(b\) ya que \(x_0=a\) y \(x_n=b\).

Teorema.

Consideramos la cuadratura \(\displaystyle\sum_{i=0}^n a_i f(x_i)\) como la fórmula de Newton-Cotes cerrada considerando \(n+1\) puntos equiespaciados en \([a,b]\), con \(x_0=a\), \(x_n=b\), \(h=\frac{b-a}{n}\). Entonces:

• Si \(n\) es par y \(f\in {\cal C}^{n+2}[a,b]\), existe un \(\xi\in (a,b)\) tal que \[ \int_a^b f(x)\, dx=\sum_{i=0}^n a_i f(x_i)+\frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!}\int_0^n t^2(t-1)\cdots (t-n)\,dt. \]

Teorema. (continuación)

• Si \(n\) es impar y \(f\in {\cal C}^{n+1}[a,b]\), existe un \(\xi\in (a,b)\) tal que \[ \int_a^b f(x)\, dx=\sum_{i=0}^n a_i f(x_i)+\frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!}\int_0^n t(t-1)\cdots (t-n)\,dt. \]

Observación. Los coeficientes \(a_i\) en las fórmulas de Newton-Cotes se pueden escribir de la siguiente manera: \[ a_i=h\int_0^n \prod_{j=0,j\neq i}^n \frac{t-j}{i-j}\, dt, \] ya que (\(L_i(x)\) son los polinomios de Lagrange) si hacemos el cambio de variable \(x=a+t h\), tenemos: \[ \begin{align*} a_i= & \int_a^b L_i(x)\, dx=\int_a^b \prod_{j=0,j\neq i}^n \frac{x-x_j}{x_i-x_j}\, dx\\ = & h\int_0^n\prod_{j=0,j\neq i}^n \frac{a+th-a-jh}{a+ih-a-jh}\, dt=h\int_0^n \prod_{j=0,j\neq i}^n \frac{t-j}{i-j}\, dt. \end{align*} \]

En resumen, llamando \(\displaystyle\alpha_{i,n}=\int_0^n \prod_{j=0,j\neq i}^n \frac{t-j}{i-j}\, dt\) las fórmulas de Newton-Cotes serían: \[ \int_a^b f(x)\, dx=h \sum_{i=0}^n \alpha_{i,n} f(x_i)+E(f,a,b), \] donde \(E(f,a,b)\) es el término del error indicado en el Teorema anterior.

Los valores \(\alpha_{i,n}\) para \(n=1,2,3,4,5\) son los siguientes:

Observaciones

• Las dos primeras filas anteriores corresponden al método de los trapecios y a la fórmula de Simpson, respectivamente.
• Es mejor siempre considerar \(n\) par ya que el término del error es de orden \(h^{n+3}\), en cambio si \(n\) es impar el término del error es de orden sólo \(h^{n+2}\). Es decir, “ganamos” un \(h\) extra tal como pasa con la fórmula de Simpson.

Vamos a hacer lo mismo que antes pero ahora no incluiremos los extremos del intervalo de integración \(a\) y \(b\) en los nodos equiespaciados de interpolación.

Los nodos serán pues \(x_i=x_0+ih\), donde \(i=0,\ldots,n\), \(h=\frac{b-a}{n+2}\) y \(x_0=a+h\): \[ \begin{align*} x_0 & =a+h,\ x_1=a+2h,\ldots, \\ x_n= & a+h+n h=a+h+\frac{(b-a)n}{n+2}=b-h. \end{align*} \]

En este caso, etiquetamos los extremos del intervalo como \(x_{-1}=a\), \(x_{n+1}=b\) y la cuadratura será: \[ \int_a^b f(x)\, dx=\int_{x_{-1}}^{x_{n+1}} f(x)\, dx \approx \sum_{i=0}^n a_i f(x_i), \] donde recordemos que los \(a_i\) valen \(\displaystyle a_i=\int_a^b L_i(x)\, dx\), donde \(L_i(x)\) es el \(i-èsimo\) polinomio de Lagrange respecto los nodos \(x_0,x_1,\ldots, x_n\).

Teorema.

Consideramos la cuadratura \(\displaystyle\sum_{i=0}^n a_i f(x_i)\) como la fórmula de Newton-Cotes abierta considerando \(n+1\) puntos equiespaciados en \([a,b]\), con \(x_{-1}=a\), \(x_{n+1}=b\), \(h=\frac{b-a}{n+2}\). Entonces:

• Si \(n\) es par y \(f\in {\cal C}^{n+2}[a,b]\), existe un \(\xi\in (a,b)\) tal que \[ \int_a^b f(x)\, dx=\sum_{i=0}^n a_i f(x_i)+\frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!}\int_{-1}^{n+1} t^2(t-1)\cdots (t-n)\,dt. \]

Teorema. (continuación)

• Si \(n\) es impar y \(f\in {\cal C}^{n+1}[a,b]\), existe un \(\xi\in (a,b)\) tal que \[ \int_a^b f(x)\, dx=\sum_{i=0}^n a_i f(x_i)+\frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!}\int_{-1}^{n+1} t(t-1)\cdots (t-n)\,dt. \]

Observaciones.

• Tenemos una expresión del error muy parecida a las fórmulas de Newton-Cotes cerradas, lo único que cambia en el término del error es que ahora integramos entre \(-1\) y \(n+1\) en vez de entre \(0\) y \(n\).
• Vemos que si \(n\) es par la precisión vuelve a ser mejor que si \(n\) es impar tal como pasaba en las fórmulas de Newton-Cotes cerradas.

• Caso \(n=0\), fórmula del punto medio. En este caso sólo tenemos un punto \(x_0=\frac{a+b}{2}\) y \(h=\frac{b-a}{2}\): \[ \begin{align*} \int_a^b f(x)\, dx= & 2hf(x_0)+\frac{h^3}{2}\int_{-1}^1 t^2\, dt= 2hf(x_0)+\frac{h^3}{2}\cdot\frac{2}{3}\cdot f''(\xi)\\ = & 2hf(x_0)+\frac{h^3}{3}\cdot f''(\xi) \end{align*} \] con \(\xi\in (a,b)\).

• \(n=1\). En este caso, \(h=\frac{b-a}{3}\), \(x_0=a+h=\frac{2a+b}{3}\), \(x_1=a+2h=\frac{a+2b}{3}\): \[ \begin{align*} & \int_a^b f(x)\, dx= \frac{3h}{2}(f(x_0)+f(x_1))+\frac{h^3 f''(\xi)}{2}\int_{-1}^2 t(t-1)\, dt\\ = & \frac{3h}{2}(f(x_0)+f(x_1))+\frac{h^3 f''(\xi)}{2}\cdot \frac{3}{2}\\ = & \frac{3h}{2}(f(x_0)+f(x_1))+\frac{3h^3 f''(\xi)}{4}. \end{align*} \] con \(\xi\in (a,b)\).

• \(n=2\). En este caso, \(h=\frac{b-a}{4}\), \(x_0=a+h=\frac{3a+b}{4}\), \(x_1=a+2h=\frac{a+b}{2}\), \(x_2=a+3h=\frac{a+3b}{4}\): \[ \begin{align*} & \int_a^b f(x)\, dx= \frac{4h}{3}(2f(x_0)-f(x_1)+2f(x_2))\\ & +\frac{h^5 f^{(4)}(\xi)}{4!}\int_{-1}^3 t^2(t-1)(t-2)\, dt\\ = & \frac{4h}{3}(2f(x_0)-f(x_1)+2f(x_2)) +\frac{h^5 f^{(4)}(\xi)}{24}\cdot \frac{112}{15}\\ = & \frac{4h}{3}(2f(x_0)-f(x_1)+2f(x_2)) +\frac{14h^5 f^{(4)}(\xi)}{45}. \end{align*} \] con \(\xi\in (a,b)\).

Las fórmulas de Newton-Cotes introducidas en la sección anterior tienen los siguientes inconvenientes:

• Se requiere mucho cálculo para obtener los coeficientes \(a_i\) para valores de \(n\) grandes.
• Están basadas en interpolar polinomios de grado cada vez más alto en nodos equiespaciados. Ya comentamos en el capítulo de interpolación que aproximar una función con polinomios de grado cada vez más alto en nodos equiespaciados no siempre produce los efectos deseados, es decir, un polinomio que tenga poco error en todos los valores del intervalo en cuestión. Recordar el fenómeno de Runge, por ejemplo.

Por los motivos anteriormente expuestos, la fórmulas de Newton-Cotes se suelen usar para valores de \(n\) bajos.

Ahora bien, si el intervalo de integración tiene una longitud grande, tendremos un error grande ya que el valor de \(h\) será también grande.

Para resolver este problema, dividiremos el intervalo de integración en pequeños subintervalos y aplicaremos la fórmula de Newton-Cotes a cada uno de los subintervalos usando que la integral sobre el intervalo grande es la suma de las integrales sobre cada uno de los subintervalos.

Dicha técnica se denomina integración compuesta.

Vamos a generalizar el ejemplo anterior.

Sea \(f\in {\cal C}^2[a,b]\) de clase \({\cal C}^2\) en un cierto intervalo \([a,b]\). Sea \(n\) un número natural (número de subintervalos en los que subdividiremos el intervalo \([a,b]\)) y \(h=\frac{b-a}{n}\). Sean \(x_i=a+ih\), \(i=0,1,\ldots,n\), los extremos de dichos subintervalos, con \(x_0=a\) y \(x_n=b\). Entonces podemos aproximar \(\displaystyle\int_a^b f(x)\,dx\) de la forma siguiente: \[ \int_a^b f(x)\, dx=\sum_{i=0}^{n-1}\int_{x_{i}}^{x_{i+1}}f(x)\, dx=\sum_{i=0}^{n-1}\left(\frac{h}{2}\left(f(x_{i})+f(x_{i+1})\right)-\frac{h^3}{12}f''(\xi_i)\right), \] con \(\xi_i\in (x_{i},x_{i+1})\).

La expresión anterior puede simplificarse de la forma siguiente: \[ \int_a^b f(x)\, dx =\frac{h}{2}\left(f(a)+2\sum_{i=1}^{n-1}f(x_i)+f(b)\right)-\frac{h^3}{12}\sum_{i=0}^{n-1}f''(\xi_i). \] Usando el Teorema de Bolzano generalizado, podemos afirmar que existe un \(\xi\in (a,b)\) tal que \[ \sum_{i=0}^{n-1}f''(\xi_i)=n f''(\xi). \]

La fórmula de trapecios compuesta queda de la forma siguiente: \[ \begin{align*} \int_a^b f(x)\, dx = & \frac{h}{2}\left(f(a)+2\sum_{i=1}^{n-1}f(x_i)+f(b)\right)-\frac{h^3}{12} n f''(\xi)\\ = & \frac{h}{2}\left(f(a)+2\sum_{i=1}^{n-1}f(x_i)+f(b)\right)-\frac{(b-a)}{12}h^2 f''(\xi), \end{align*} \] donde hemos usado que \(h\cdot n=b-a\).

La aproximación de la fórmula de trapecios compuesta es: \[ T(h)=\frac{h}{2}\left(f(a)+2\sum_{i=1}^{n-1}f(x_i)+f(b)\right), \] con error: \[ E(f,a,b)=-\frac{(b-a)}{12}h^2 f''(\xi). \]

Observación.

La fórmula de trapecios compuesta tiene error \(O(h^2)\).

• INPUT a,b,n. (damos los extremos del intervalo de integración y el número de subintervalos)
• Set h=(b-a)/n. (calculamos el valor de \(h\))
• Set S0=f(a)+f(b). (calculamos \(f(a)+f(b)\))
• Set S=0 (en la variable \(S\) vamos a guardar la suma \(\displaystyle\sum_{i=1}^{n-1}f(x_i)\))
• For i=1,...,n-1 (calculamos la suma anterior)
  • Set X=a+i*h.
  • Set S=S+f(X).
• Set Iaprox = (h/2)*(S0+2*S). (calculamos el valor de la integral aproximada)
• Print (Iaprox). (damos dicho valor)
• STOP.

Vamos a aplicar la técnica anterior pero en lugar de considerar la fórmula de los trapecios, consideraremos la fórmula de Simpson.

Recordemos que para aplicar la fórmula de Simpson en un intervalo \([c,d]\) necesitábamos tres puntos: \(c\), \(d\) y \(\frac{c+d}{2}\): \[ \int_c^d f(x)\, dx=\frac{h}{3}\left(f(c)+4f\left(\frac{c+d}{2}\right)+f(d)\right)-\frac{h^5}{90}f^{(4)}(\xi), \] donde \(\xi\in (c,d)\) y \(h=\frac{d-c}{2}\).

Sea \(f\in {\cal C}^4[a,b]\) de clase \({\cal C}^4\) en un cierto intervalo \([a,b]\). Sea \(n\) un número par, \(h=\frac{b-a}{n}\) y \(x_i=a+ih\), con \(i=0,1,\ldots,n\). Dividimos el intervalo \([a,b]\) en los subintervalos siguientes: \[ [a,b]=[x_0,x_2]\cup [x_2,x_4]\cup\cdots\cup [x_{n-2},x_n], \] es decir, aplicaremos la fórmula de Simpson a cada subintervalo de la forma \([x_{2i},x_{2i+2}]\), \(i=0,1,\ldots,\frac{n-2}{2}\): \[ \int_{x_{2i}}^{x_{2i+2}} f(x)\, dx=\frac{h_i}{3}\left(f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})\right)-\frac{h_i^5}{90}f^{(4)}(\xi_i), \] donde \(h_i=\frac{x_{2i+2}-x_{2i}}{2}=\frac{a+(2i+2)h-a-2ih}{2}=\frac{2h}{2}=h\) y \(\xi_i\in (x_{2i},x_{2i+2})\).

La fórmula resultante será la siguiente: \[ \begin{align*} \int_a^b f(x)\, dx=& \sum_{i=0}^{\frac{n-2}{2}}\int_{x_{2i}}^{x_{2i+2}}f(x)\, dx \\ =& \sum_{i=0}^{\frac{n-2}{2}}\left(\frac{h}{3}\left(f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})\right)-\frac{h^5}{90}f^{(4)}(\xi_i)\right)\\ = & \frac{h}{3}\left(f(a)+2\sum_{i=1}^{\frac{n-2}{2}}f(x_{2i})+4\sum_{i=0}^{\frac{n-2}{2}}f(x_{2i+1})+f(b)\right)\\ & -\frac{h^5}{90}\sum_{i=0}^{\frac{n-2}{2}}f^{(4)}(\xi_i). \end{align*} \]

Usando el Teorema de Bolzano generalizado, podemos afirmar que existe un \(\xi\in (a,b)\) tal que: \[ \sum_{i=0}^{\frac{n-2}{2}}f^{(4)}(\xi_i)=\left(1+\frac{n-2}{2}\right)f^{(4)}(\xi)=\frac{n}{2}f^{(4)}(\xi). \]

La fórmula de Simpson compuesta será, pues: \[ \begin{align*} \int_a^b f(x)\, dx=& \frac{h}{3}\left(f(a)+2\sum_{i=1}^{\frac{n-2}{2}}f(x_{2i})+4\sum_{i=0}^{\frac{n-2}{2}}f(x_{2i+1})+f(b)\right)\\ & -\frac{h^5}{90}\cdot\frac{n}{2}f^{(4)}(\xi) \\ =& \frac{h}{3}\left(f(a)+2\sum_{i=1}^{\frac{n-2}{2}}f(x_{2i})+4\sum_{i=0}^{\frac{n-2}{2}}f(x_{2i+1})+f(b)\right)\\ & -\frac{(b-a)}{180}\cdot h^4 f^{(4)}(\xi). \end{align*} \]

Observación.

La fórmula de Simpson compuesta tiene error \(O(h^4)\).

• INPUT a,b,n. (damos los extremos del intervalo de integración y el valor \(n\))
• Set h=(b-a)/n. (calculamos el valor de \(h\))
• Set S0=f(a)+f(b). (calculamos \(f(a)+f(b)\))
• Set S1=0. (en la variable \(S1\) vamos a guardar la suma \(\displaystyle\sum_{i=1}^{\frac{n-2}{2}}f(x_{2i})\)
• Set S2=0. (y en \(S2\), \(\displaystyle \sum_{i=0}^{\frac{n-2}{2}}f(x_{2i+1})\))

• For i=1,...,n-1 (calculamos las sumas anteriores)
  • Set X=a+i*h.
  • If i par then
    • Set S1=S1+f(X)
  • Else
    • Set S2=S2+f(X).
• Set Iaprox = (h/3)*(S0+2*S1+4*S2). (calculamos el valor de la integral aproximada)
• Print (Iaprox). (damos dicho valor)
• STOP.

Los algoritmos para calcular aproximaciones de derivadas son inestables como ya estudiamos en su momento, ver la parte de inestabilidad del error de redondeo.

En pocas palabras, vimos que no tenía sentido considerar valores de \(h\) muy pequeños ya que el error de redondeo debido a la evaluación de la función \(f\) tenía una parte que tendía a infinito o “explotaba” cuando \(h\) tendía a cero.

Nos preguntamos lo mismo con los algoritmos de integración numérica. ¿Son estables o inestables?

Es decir, si consideramos una fórmula de integración compuesta y aumentamos el número de subintervalos considerados, ¿cómo afectará este hecho al error global cometido, teniendo en cuenta aparte del error de truncamiento del método, los errores de redondeo cometidos al evaluar la función \(f\)?

Estudiemos la estabilidad para el método de trapecios compuesto: \[ T(h)=\frac{h}{2}\left(f(a)+2\sum_{i=1}^{n-1}f(x_i)+f(b)\right). \] Sean \(e_i\), \(i=0,1,\ldots,n\) los errores de redondeo cometidos al evaluar \(f(x_i)\), es decir: \[ f(x_i)=\tilde{f}(x_i)+e_i,\ i=0,1,\ldots,n, \] donde \(\tilde{f}(x_i)\) es el valor que calculamos realmente y \(f(x_i)\) su valor exacto.

Entonces el error de redondeo cometido por el método para un \(h\) determinado está acotado por: \[ e(h)\leq \frac{h}{2}\left(|e_0|+2\sum_{i=1}^{n-1} |e_i|+|e_n|\right). \] Si todos los errores \(e_i\) están acotados por \(e\), tenemos que: \[ e(h)\leq \frac{h}{2}(1+2(n-1)+1)e=nhe=(b-a)e. \] El error de redondeo \(e(h)\) está acotado por una constante independiente de \(h\) y, por tanto, de \(n\). Esto significa que si aumentamos \(n\) o disminuimos \(h\), el número de operaciones requeridas no aumentará el error de redondeo global del método lo que significa que el método es estable si \(h\) tiende a cero.

\[ S(h)=\frac{h}{3}\left(f(a)+2\sum_{i=1}^{\frac{n-2}{2}}f(x_{2i})+4\sum_{i=0}^{\frac{n-2}{2}}f(x_{2i+1})+f(b)\right). \] En este caso, el error global de redondeo del método estará acotado por: \[ \begin{align*} e(h)\leq & \frac{h}{3}\left(e+2\frac{(n-2)}{2} e+4\left(\frac{(n-2)}{2}+1\right)e+e\right) \\ = & \frac{h}{3}\cdot 3ne=hne=(b-a)e. \end{align*} \]

Nos pasa lo mismo que en el método de trapecios compuesto, lo que deducimos que el método de Simpson compuesto también es estable cuando \(h\) tiende a cero o si aumentamos \(n\).

Hemos dado un procedimiento para estudiar la estabilidad de un método numérico de integración.

En general, como los métodos de integración compuesta tendrán expresiones similares a la fórmula de los trapecios compuesta o a la fórmula de Simpson serán estables. De todas formas, para estar seguros, podemos aplicar el procedimiento anterior y estudiar su estabilidad.

El método de integración de Romberg consiste en aplicar la extrapolación de Richardson a un método de integración numérica.

En esta sección, para ilustrar su funcionamiento, aplicaremos la extrapolación de Richardson al método de los trapecios compuesto.

En primer lugar, se puede demostrar el resultado siguiente:

Teorema. Expresión del error del método de los trapecios compuesto.

Sea \(f\in {\cal C}^2[a,b]\) de clase \({\cal C}^2\) en un cierto intervalo \([a,b]\). Sea \(n\) un número natural y \(h=\frac{b-a}{n}\). Sean \(x_i=a+ih\), \(i=0,1,\ldots,n\), los extremos de dichos subintervalos, con \(x_0=a\) y \(x_n=b\). Entonces el método de los trapecios compuesto admite una expresión como la que sigue: \[ \begin{align*} \int_a^b f(x)\, dx= & \frac{h}{2}\left(f(a)+2\sum_{i=1}^{n-1}f(x_i)+f(b)\right)\\ & +k_1 h^2+k_2 h^4+\cdots + k_n h^{2n}+\cdots \end{align*} \]

Entonces sólo aparecen exponentes pares en los exponentes de la \(h\) en la expresión del error.

Llamamos \(T_0=T(h)=\frac{h}{2}\left(f(a)+2\sum_{i=1}^{n-1}f(x_i)+f(b)\right)\). Para hallar las aproximaciones de orden superior, usábamos la expresión siguiente: \[ T_{n+1}(h)=\frac{4^{n+1}T_n\left(\frac{h}{2}\right)-T_n(h)}{4^{n+1}-1}=C+k_{n+2}\frac{p_{n+1}}{q_{n+1}}h^{2(n+2)}+\cdots \]

• INPUT a,b,n0,n. (damos los extremos del intervalo de integración, el valor \(n_0\) del número de subintervalos inicial a considerar y el orden \(h^{2n}\) al que queremos llegar)
• For i=0,...,n do (calculamos primero \(T_0\left(\frac{h_0}{2^i}\right)\))
  • Set ni=n0*2^i. (calculamos el número de subintervalos para el nivel \(i\)-ésimo)
  • Set h=(b-a)/(ni). (calculamos el \(h\) para el nivel \(i\)-ésimo)
  • Set T_0[i]=(h/2)*(f(a)+f(b)+Sum(f(a+k*h),k=1,...,ni-1)). (calculamos \(T_0\left(\frac{h_0}{2^i}\right)\))

• For k=1,...,n do (el índice \(k\) nos indicará el nivel donde estamos. Guardaremos los valores \(T_k\left(\frac{h_0}{2^i}\right)\), \(i=0,\ldots,n-k\) en el vector \(T_k\))
  • For i=0,...,n-k do (bucle para el índice \(i\))
    • Set ni=n0*2^i.
    • Set h=(b-a)/(ni).
    • Set T_k[i]=(4^k*T_{k-1}(h/2)-T_{k-1}(h))/(4^k-1).
• Print T_n[0]. (damos la aproximación de orden \(h^{2n}\))
• STOP.

Las fórmulas de cuadratura numérica vistas hasta el momento tienen la expresión siguiente: \[ \int_a^b f(x)\, dx\approx \sum_{i=0}^n a_i\cdot f(x_i), \] donde \(a_i\), \(i=0,\ldots,n\) son unos coeficientes a determinar de cara a tener el menor error posible o la mayor precisión posible y \(x_i\), \(i=0,1,\ldots, n\) son los nodos que consideramos equiespaciados y que vienen dados en función del intervalo de integración \([a,b]\).

Es decir dado \([a,b]\) y valor \(n\), se determinan los nodos equiespaciados \(x_i\), \(i=0,1,\ldots, n\).

En esta sección, vamos a introducir unas fórmulas de cuadratura donde se admite más libertad en el sentido de que los nodos de integración \(x_i\) no están predeterminados y forman parte de los parámetros a determinar en la fórmula de cuadratura numérica correspondiente: \[ \int_a^b f(x)\, dx\approx \sum_{i=0}^n a_i\cdot f(x_i). \] Es decir, no se presupone que los nodos \(x_i\) sean equiespaciados sino que se pueden escoger de tal forma que la expresión anterior tenga la precisión más alta posible.

Supongamos que \(n=0\) y el intervalo \([a,b]\) es \([-1,1]\). Nos planteamos hallar una fórmula de cuadratura de la forma: \(\displaystyle\int_{-1}^1 f(x)\, dx \approx a_0 f(x_0)\) que tenga la precisión más alta posible.

Como tenemos dos parámetros, \(a_0\) y \(x_0\), podemos exigir precisión igual a \(1\): \[ \begin{align*} \int_{-1}^1 1\, dx = & 2=a_0,\\ \int_{-1}^1 x\, dx= & \left[\frac{x^2}{2}\right]_{-1}^1 = 0=a_0\cdot x_0. \end{align*} \] Los valores de los parámetros \(a_0\) y el nodo \(x_0\) son los siguientes: \(a_0=2,\ x_0=0\).

Entonces, la fórmula de cuadratura: \[ \int_{-1}^1 f(x)\, dx\approx 2\cdot f\left(0\right), \] tiene precisión \(1\).

Consideremos \(n=1\). En este caso, buscamos coeficientes \(a_0\) y \(a_1\) y nodos \(x_0\) y \(x_1\) tal que la fórmula de cuadratura \(\displaystyle\int_{-1}^1 f(x)\, dx\approx a_0 f(x_0)+a_1 f(x_1)\) tenga precisión lo más alta posible.

Como tenemos \(4\) parámetros, \(a_0,a_1,x_0\) y \(x_1\) podemos exigir precisión \(3\):

\[ \begin{align*} \int_{-1}^1 1\, dx = & 2=a_0+a_1,\\ \int_{-1}^1 x\, dx= & 0=a_0\cdot x_0+a_1\cdot x_1,\\ \int_{-1}^1 x^2\, dx= & \frac{2}{3}=a_0\cdot x_0^2+a_1\cdot x_1^2,\\ \int_{-1}^1 x^3\, dx= & 0=a_0\cdot x_0^3+a_1\cdot x_1^3. \end{align*} \] Las soluciones del sistema anterior son: \(a_0=a_1=1\), \(x_0=\frac{1}{\sqrt{3}}\) y \(x_1=-\frac{1}{\sqrt{3}}\).

Entonces, la fórmula de cuadratura: \[ \int_{-1}^1 f(x)\, dx\approx f\left(-\frac{1}{\sqrt{3}}\right)+f\left(\frac{1}{\sqrt{3}}\right), \] tiene precisión \(3\).

Por último, consideremos \(n=2\). En este caso, buscamos coeficientes \(a_0\), \(a_1\) y \(a_2\) y nodos \(x_0\), \(x_1\) y \(x_2\) tal que la fórmula de cuadratura \(\displaystyle\int_{-1}^1 f(x)\, dx\approx a_0 f(x_0)+a_1 f(x_1)+a_2 f(x_2)\) tenga precisión lo más alta posible.

Como tenemos \(6\) parámetros, \(a_0,a_1,a_2,x_0,x_1\) y \(x_2\) podemos exigir precisión \(5\):

\[ \begin{align*} \int_{-1}^1 1\, dx = & 2=a_0+a_1+a_2,\\ \int_{-1}^1 x\, dx= & 0=a_0\cdot x_0+a_1\cdot x_1+a_2 x_2,\\ \int_{-1}^1 x^2\, dx= & \frac{2}{3}=a_0\cdot x_0^2+a_1\cdot x_1^2+a_2\cdot x_2^2,\\ \int_{-1}^1 x^3\, dx= & 0=a_0\cdot x_0^3+a_1\cdot x_1^3+a_2\cdot x_2^3,\\ \int_{-1}^1 x^4\, dx= & \frac{2}{5}=a_0\cdot x_0^4+a_1\cdot x_1^4+a_2\cdot x_2^4,\\ \int_{-1}^1 x^5\, dx= & 0=a_0\cdot x_0^5+a_1\cdot x_1^5+a_2\cdot x_2^5. \end{align*} \]

Las soluciones del sistema anterior son: \[ a_0=\frac{5}{9},\ a_1=\frac{5}{9},\ a_2=\frac{8}{9},\ x_0=\sqrt{\frac{3}{5}},\ x_1=-\sqrt{\frac{3}{5}},\ x_2=0. \]

Entonces, la fórmula de cuadratura: \[ \int_{-1}^1 f(x)\, dx\approx \frac{5}{9}\cdot f\left(-\sqrt{\frac{3}{5}}\right)+\frac{8}{9}\cdot f(0)+\frac{5}{9}\cdot f\left(\sqrt{\frac{3}{5}}\right), \] tiene precisión \(5\).

¿Cómo podemos obtener fórmulas para valores de \(n\) más grandes y generalizar dichas fórmulas para cualquier intervalo \([a,b]\)?

Para poder hacerlo, necesitamos conocer los polinomios de Legendre:

Definición de los polinomios de Legendre.

Las sucesión de polinomios \(P_0(x),P_1(x),\ldots,P_n(x),\ldots\) forma la sucesión de los polinomios de Legendre si se cumplen las condiciones siguientes:

• el grado de \(P_n(x)\) es \(n\),
• \(\displaystyle\int_{-1}^1 P(x)\cdot P_n(x)\, dx=0\), para cualquier polinomio \(P(x)\) de grado menor o igual que \(n-1\).

Como los polinomios de Legendre están definidos salvo una constante multiplicativa, consideraremos polinomios de Legendre mónicos, es decir, polinomios \(P_n(x)\) cuyo coeficiente de \(x^n\) es \(1\).

Los primeros polinomios de Legendre son los siguientes: \[ \begin{align*} P_0(x)= & 1,\ P_1(x)=x,\ P_2(x)=x^2-\frac{1}{3},\ P_3(x)=x^3-\frac{3}{5}x,\\ P_4(x)= & x^4-\frac{6}{7}x^2+\frac{3}{35}. \end{align*} \]

Proposición. Expresiones de los polinomios de Legendre.

• Los polinomios de Legendre mónicos verifican la siguiente relación de recurrencia (fórmula de Bonet): \[ P_{n+1}(x)=x\cdot P_n(x)-\frac{n^2}{(4n^2-1)}P_{n-1}(x),\ n=1,2,\ldots \]
• Los polinomios de Legendre también pueden calcularse de la forma siguiente (fórmula de Rodrigues): \[ P_n(x)=\frac{1}{2^n\cdot (2n-1)!!}\frac{d^n}{dx^n}(x^2-1)^n. \]

Las expresiones anteriores nos dan una manera para ir calculando todos los polinomios de Legendre.

Las raíces de los polinomios de Legendre cumplen la condición siguiente:

Proposición. Raíces de los polinomios de Legendre. Sea \(P_n(x)\) el polinomio de Legendre de grado \(n\). Entonces sus raíces \(x_1,\ldots,x_n\) son reales, todas están en el intervalo \((-1,1)\) y son simétricas respecto al origen, es decir, si \(x_i\) es una raíz de \(P_n(x)\), \(P_n(x_i)=0\), entonces \(-x_i\) también lo es, es decir, \(P_n(-x_i)=0\).

Observación. Se puede demostrar a partir de la fórmula de Bonet que los polinomios de Legendre \(P_n(x)\) son pares si \(n\) es par e impares si \(n\) es impar. Es decir: \[ P_{2k}(-x)=P_{2k}(x),\quad P_{2k+1}(-x)=-P_{2k+1}(x), \] para todo valor de \(k\).

El teorema siguiente nos dice que si escogemos como nodos los ceros de los polinomios de Legengre, obtenemos una cuadratura de precisión \(2n-1\) en el intervalo \([-1,1]\), donde \(n\) es el grado del polinomio de Legendre considerado.

Teorema. Cuadratura gaussiana Sea \(P_n(x)\) el polinomio de Legendre de grado \(n\). Sean \(x_1,x_2,\ldots,x_n\) las raíces del mismo en el intervalo \((-1,1)\). Consideramos los números \(c_i\) definidos de la forma siguiente: \[ c_i=\int_{-1}^1 \prod_{j=1,\\ j\neq i}^n \frac{x-x_j}{x_i-x_j}\, dx. \] Entonces si \(P(x)\) es un polinomio de grado menor o igual que \(2n-1\), \[ \int_{-1}^1 P(x)\,dx=\sum_{i=1}^n c_i P(x_i). \]

Los nodos \(x_i\), \(i=1,\ldots,n\) de los primeros polinomios de la cuadratura gaussiana son los siguientes:

Los coeficientes \(c_i\), \(i=1,\ldots,n\) de las primeras expresiones de la cuadratura gaussiana son los siguientes: (están asociados en el mismo orden que los nodos \(x_i\))

Error en la cuadratura gaussiana.

Sea \(n\) un entero positivo y \(f\in {\cal C}^{2n}[-1,1]\) una función de clase \(2n\) en el intervalo \([-1,1]\). Entonces existe un valor \(\xi\in (-1,1)\) tal que el error de la cuadratura gaussiana viene dado por: \[ \int_{-1}^1 f(x)\, dx =\sum_{i=1}^n c_i f(x_i)+\frac{f^{(2n)}(\xi)}{(2n)!}\int_{-1}^1 P_n(x)^2\, dx, \] donde \(P_n(x)\) es el polinomio de Legendre de grado \(n\).

Si en lugar de integrar en el intervalo \([-1,1]\), tenemos que integrar en un intervalo cualquiera \([a,b]\), podemos considerar el cambio de variable siguiente \(t=\varphi(x)=\left(\frac{b-a}{2}\right)x+\frac{a+b}{2}\) que nos pasa del intervalo \([-1,1]\) al intervalo \([a,b]\): \[ \begin{align*} \varphi: [-1,1]\longrightarrow & [a,b]\\ x \longrightarrow & t=\varphi(x)=\left(\frac{b-a}{2}\right)x+\frac{a+b}{2}. \end{align*} \]

Usando el cambio anterior, transformamos la integral de aproximar una función \(f(x)\) en el intervalo \([a,b]\) en una integral en el intervalo \([-1,1]\) de la que sabemos aproximar usando la cuadratura gaussiana: \[ \begin{align*} \int_a^b f(t)\, dt =& \frac{(b-a)}{2}\int_{-1}^1 f\left(\left(\frac{b-a}{2}\right)x+\frac{a+b}{2}\right)\, dx\\ \approx & \frac{(b-a)}{2}\sum_{i=1}^n c_i f\left(\left(\frac{b-a}{2}\right)x_i+\frac{a+b}{2}\right), \end{align*} \] donde \(x_i\), \(i=1,\ldots,n\) son los \(n\) ceros del \(n\)-ésimo polinomio de Legendre \(P_n(x)\).