3: Integración

Queremos construir algoritmos numéricos que puedan realizar integrales definidas de la forma

$$I=\int_{a}^{b} f(x) d x . \nonumber$$

El cálculo numérico de estas integrales definidas se denomina integración numérica, cuadratura numérica o, más simplemente, cuadratura.

Fórmulas elementales

Primero consideramos la integración de 0 a$h$, con$h$ pequeños, para servir como los bloques de construcción para la integración en dominios más grandes. Aquí definimos$I_{h}$ como la siguiente integral:

$$I_{h}=\int_{0}^{h} f(x) d x . \nonumber$$

Para realizar esta integral, consideramos una expansión de la serie Taylor de$f(x)$ aproximadamente el valor$x=h / 2$:

$$\begin{aligned} f(x)=f(h / 2)+(x-h / 2) f^{\prime}(h / 2) &+\frac{(x-h / 2)^{2}}{2} f^{\prime \prime}(h / 2) \\ +& \frac{(x-h / 2)^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{(x-h / 2)^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots \end{aligned} \nonumber$$

Regla de punto medio

La regla del punto medio hace uso solo del primer término en la expansión de la serie Taylor. Aquí, determinaremos el error en esta aproximación. Integrando,

$$\begin{aligned} I_{h}=h f(h / 2)+\int_{0}^{h}((x-&h / 2) f^{\prime}(h / 2)+\frac{(x-h / 2)^{2}}{2} f^{\prime \prime}(h / 2) \\ &\left.+\frac{(x-h / 2)^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{(x-h / 2)^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d x . \end{aligned} \nonumber$$

Cambiando variables dejando$y=x-h / 2$ y$d y=d x$, y simplificando la integral dependiendo de si el integrando es par o impar, tenemos

$$\begin{aligned} I_{h}=& h f(h / 2) \\ &+\int_{-h / 2}^{h / 2}\left(y f^{\prime}(h / 2)+\frac{y^{2}}{2} f^{\prime \prime}(h / 2)+\frac{y^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{y^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d y \\ =& h f(h / 2)+\int_{0}^{h / 2}\left(y^{2} f^{\prime \prime}(h / 2)+\frac{y^{4}}{12} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d y . \end{aligned} \nonumber$$

Las integrales que necesitamos aquí son

$$\int_{0}^{\frac{h}{2}} y^{2} d y=\frac{h^{3}}{24}, \quad \int_{0}^{\frac{h}{2}} y^{4} d y=\frac{h^{5}}{160} . \nonumber$$

Por lo tanto,

$$I_{h}=h f(h / 2)+\frac{h^{3}}{24} f^{\prime \prime}(h / 2)+\frac{h^{5}}{1920} f^{\prime \prime \prime \prime}(h / 2)+\ldots \nonumber$$

Regla trapezoidal

De la expansión de la serie Taylor de$f(x)$ aproximadamente$x=h / 2$, tenemos

$$f(0)=f(h / 2)-\frac{h}{2} f^{\prime}(h / 2)+\frac{h^{2}}{8} f^{\prime \prime}(h / 2)-\frac{h^{3}}{48} f^{\prime \prime \prime}(h / 2)+\frac{h^{4}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots, \nonumber$$

y

$$f(h)=f(h / 2)+\frac{h}{2} f^{\prime}(h / 2)+\frac{h^{2}}{8} f^{\prime \prime}(h / 2)+\frac{h^{3}}{48} f^{\prime \prime \prime}(h / 2)+\frac{h^{4}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber$$

Sumando y multiplicando por$h / 2$ obtenemos

$$\frac{h}{2}(f(0)+f(h))=h f(h / 2)+\frac{h^{3}}{8} f^{\prime \prime}(h / 2)+\frac{h^{5}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber$$

Ahora sustituimos el primer término en el lado derecho usando la fórmula de regla de punto medio:

$$\begin{aligned} \frac{h}{2}(f(0)+f(h))=\left(I_{h}-\frac{h^{3}}{24} f^{\prime \prime}(h / 2)-\frac{h^{5}}{1920} f^{\prime \prime \prime \prime}(h / 2)\right) &\, \\ &+\frac{h^{3}}{8} f^{\prime \prime}(h / 2)+\frac{h^{5}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots, \end{aligned} \nonumber$$

y resolviendo para$I_{h}$, encontramos

$$I_{h}=\frac{h}{2}(f(0)+f(h))-\frac{h^{3}}{12} f^{\prime \prime}(h / 2)-\frac{h^{5}}{480} f^{\prime \prime \prime \prime}(h / 2)+\ldots \nonumber$$

La regla de Simpson

Para obtener la regla de Simpson, combinamos la regla de punto medio y trapezoidal para eliminar el término de error proporcional a$h^{3}$. Multiplicando (3.4) por dos y sumando a (3.8), obtenemos

$$3 I_{h}=h\left(2 f(h / 2)+\frac{1}{2}(f(0)+f(h))\right)+h^{5}\left(\frac{2}{1920}-\frac{1}{480}\right) f^{\prime \prime \prime \prime}(h / 2)+\ldots, \nonumber$$

o

$$I_{h}=\frac{h}{6}(f(0)+4 f(h / 2)+f(h))-\frac{h^{5}}{2880} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber$$

Por lo general, la regla de Simpson se escribe considerando los tres puntos consecutivos$0, h$ y$2 h$. Sustituyendo$h \rightarrow 2 h$, obtenemos el resultado estándar

$$I_{2 h}=\frac{h}{3}(f(0)+4 f(h)+f(2 h))-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(h)+\ldots \nonumber$$

Reglas compuestas

Ahora utilizamos nuestras fórmulas elementales obtenidas para (3.2) para realizar la integral dada por (3.1).

Regla trapezoidal

Suponemos que la función$f(x)$ es conocida en los$n+1$ puntos etiquetados como$x_{0}, x_{1}, \ldots, x_{n}$, con los puntos finales dados por$x_{0}=a$ y$x_{n}=b$. Definir

$$f_{i}=f\left(x_{i}\right), \quad h_{i}=x_{i+1}-x_{i} . \nonumber$$

Entonces la integral de (3.1) puede descomponerse como

$$\begin{aligned} \int_{a}^{b} f(x) d x &=\sum_{i=0}^{n-1} \int_{x_{i}}^{x_{i+1}} f(x) d x \\ &=\sum_{i=0}^{n-1} \int_{0}^{h_{i}} f\left(x_{i}+s\right) d s, \end{aligned} \nonumber$$

donde la última igualdad surge del cambio de variables$s=x-x_{i}$. Aplicando la regla trapezoidal a la integral, tenemos

$$\int_{a}^{b} f(x) d x=\frac{1}{2} \sum_{i=0}^{n-1} h_{i}\left(f_{i}+f_{i+1}\right) . \nonumber$$

Si los puntos no están espaciados uniformemente, digamos porque los datos son valores experimentales, entonces el$h_{i}$ puede diferir para cada valor de$i$ y (3.13) se va a usar directamente.

No obstante, si los puntos están espaciados uniformemente, digamos porque se$f(x)$ pueden computar, tenemos$h_{i}=h$, independientemente de$i$. Entonces podemos definir

$$x_{i}=a+i h, \quad i=0,1, \ldots, n ; \nonumber$$

y como el punto final$b$ satisface$b=a+n h$, tenemos

$$h=\frac{b-a}{n} . \nonumber$$

La regla trapezoidal compuesta para puntos de espacio uniformemente luego se convierte en

$$\begin{align} \nonumber \int_{a}^{b} f(x) d x &=\frac{h}{2} \sum_{i=0}^{n-1}\left(f_{i}+f_{i+1}\right) \\ &=\frac{h}{2}\left(f_{0}+2 f_{1}+\cdots+2 f_{n-1}+f_{n}\right) . \end{align} \nonumber$$

El primer y último término tienen un múltiplo de uno; todos los demás términos tienen un múltiplo de dos; y la suma completa se multiplica por$h / 2$

Figura 3.1: Cuadratura adaptativa Simpson: Nivel$1 .$

La regla de Simpson

Aquí consideramos la regla compuesta de Simpson para puntos espaciados uniformemente. Aplicamos la regla de Simpson a intervalos de$2 h$, comenzando$a$ y terminando en$b$:

$$\begin{aligned} \int_{a}^{b} f(x) d x=\frac{h}{3}\left(f_{0}+4 f_{1}+f_{2}\right)+\frac{h}{3}\left(f_{2}+4 f_{3}+f_{4}\right)+\ldots & \\ &+\frac{h}{3}\left(f_{n-2}+4 f_{n-1}+f_{n}\right) . \end{aligned} \nonumber$$

Tenga en cuenta que$n$ debe ser parejo para que este esquema funcione. Combinando términos, tenemos

$$\int_{a}^{b} f(x) d x=\frac{h}{3}\left(f_{0}+4 f_{1}+2 f_{2}+4 f_{3}+2 f_{4}+\cdots+4 f_{n-1}+f_{n}\right) . \nonumber$$

El primer y último término tienen un múltiplo de uno; los términos indexados pares tienen un múltiplo de 2; los términos indexados impares tienen un múltiplo de 4; y la suma completa se multiplica por$h / 3$.

Integración adaptativa

La útil función de MATLAB quad.m realiza integración numérica utilizando cuadratura adaptativa Simpson. La idea es dejar que el cálculo mismo decida sobre el tamaño de cuadrícula requerido para lograr un cierto nivel de precisión. Además, el tamaño de la cuadrícula no necesita ser el mismo en toda la región de integración.

Comenzamos la integración adaptativa en lo que se denomina Nivel 1. Los puntos uniformemente espaciados en los que$f(x)$ se va a evaluar la función se muestran en la Fig. 3.1. La distancia entre los puntos$a$ y$b$ se toma para ser$2 h$, de manera que

$$h=\frac{b-a}{2} . \nonumber$$

Integración usando la regla de Simpson (3.11) con$h$ rendimientos de tamaño de cuadrícula para la integral$I$,

$$I=\frac{h}{3}(f(a)+4 f(c)+f(b))-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(\xi), \nonumber$$

donde$\xi$ hay algún valor satisfactorio$a \leq \xi \leq b$. Integración usando la regla de Simpson dos veces con$h / 2$ rendimientos de tamaño de cuadrícula

$$I=\frac{h}{6}(f(a)+4 f(d)+2 f(c)+4 f(e)+f(b))-\frac{(h / 2)^{5}}{90} f^{\prime \prime \prime \prime}\left(\xi_{l}\right)-\frac{(h / 2)^{5}}{90} f^{\prime \prime \prime \prime}\left(\xi_{r}\right), \nonumber$$

con$\xi_{l}$ y$\xi_{r}$ algunos valores satisfactorios$a \leq \xi_{l} \leq c$ y$c \leq \xi_{r} \leq b$.

Ahora definimos las dos aproximaciones a la integral por

$$\begin{aligned} &S_{1}=\frac{h}{3}(f(a)+4 f(c)+f(b)), \\ &S_{2}=\frac{h}{6}(f(a)+4 f(d)+2 f(c)+4 f(e)+f(b)), \end{aligned} \nonumber$$

y los dos errores asociados por

$$\begin{aligned} &E_{1}=-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(\xi), \\ &E_{2}=-\frac{h^{5}}{2^{5} \cdot 90}\left(f^{\prime \prime \prime \prime}\left(\xi_{l}\right)+f^{\prime \prime \prime \prime}\left(\xi_{r}\right)\right) . \end{aligned} \nonumber$$

Ahora nos preguntamos si el valor de$S_{2}$ para la integral es lo suficientemente exacto, o ¿debemos afinar aún más el cálculo e ir al Nivel 2? Para responder a esta pregunta, hacemos la aproximación simplificadora de que todas las derivadas de cuarto orden$f(x)$ en términos de error son iguales; es decir,

$$f^{\prime \prime \prime \prime}(\xi)=f^{\prime \prime \prime \prime}\left(\xi_{l}\right)=f^{\prime \prime \prime \prime}\left(\xi_{r}\right)=C . \nonumber$$

Entonces

$$\begin{aligned} &E_{1}=-\frac{h^{5}}{90} C, \\ &E_{2}=-\frac{h^{5}}{2^{4} \cdot 90} C=\frac{1}{16} E_{1} . \end{aligned} \nonumber$$

Ahora como la integral es igual a la aproximación más su error asociado,

$$S_{1}+E_{1}=S_{2}+E_{2} \text {, } \nonumber$$

y desde

$$E_{1}=16 E_{2} \text {, } \nonumber$$

podemos derivar una estimación para el término de error$E_{2}$:

$$E_{2}=\frac{1}{15}\left(S_{2}-S_{1}\right) \text {. } \nonumber$$

Por lo tanto, dado algún valor específico de la tolerancia tol, si

$$\left|\frac{1}{15}\left(S_{2}-S_{1}\right)\right|<\text { tol, } \nonumber$$

entonces podemos aceptar$S_{2}$ como$I$. Si la estimación del error es mayor en magnitud que tol, entonces procedemos al Nivel 2. El cálculo en el Nivel 2 divide aún más el intervalo de integración de$a$ a$b$ en los dos intervalos de integración$a$ a$c$ y$c$ a$b$, y continúa con el procedimiento anterior independientemente en ambas mitades. La integración se puede detener en cualquiera de las dos mitades siempre que la tolerancia sea menor que tol/2 (ya que la suma de ambos errores debe ser menor que tol). De lo contrario, cualquiera de las dos mitades puede pasar al Nivel 3, y así sucesivamente.

Como nota al margen, los dos valores$I$ dados anteriormente (para la integración con el tamaño del paso$h$ y se$h / 2)$ pueden combinar para dar un valor más preciso para I dado por

$$I=\frac{16 S_{2}-S_{1}}{15}+\mathrm{O}\left(h^{7}\right), \nonumber$$

donde los términos de error de$\mathrm{O}\left(h^{5}\right)$ aproximadamente cancelan. Este almuerzo gratuito, por así decirlo, se llama extrapolación de Richardson.