9.3: Regla de Simpson

Nuestro análisis de la regla del punto medio y la regla del trapecio mostró que ambos métodos tienen error\(O(1/N^2)\) numérico. En ambos casos, el error se puede rastrear a la misma fuente: el hecho de que la estimación integral sobre cada segmento difiera de la serie Taylor da como resultado el término de segundo orden, el término proporcional a\(f''(x_n)\). Esto sugiere una manera de mejorar el resultado de la integración numérica: podríamos tomar un promedio ponderado de la regla del punto medio y la regla del trapecio, ¡de tal manera que los errores numéricos de segundo orden de los dos esquemas se cancelen entre sí! Este es el método de integración numérica conocido como regla de Simpson.

Para ser precisos, volvamos a considerar un par de segmentos adyacentes, que se encuentran entre los puntos de discretización igualmente espaciados\(\left\{x_{n-1}, x_n, x_{n+1}\right\}\). Como se derivó anteriormente, la integral sobre estos segmentos se puede expandir Taylor como

\[\mathcal{I}_n = \; 2 f(x_n) \Delta x \;+\; \frac{f''(x_n)}{3} \Delta x^3 + O(\Delta x^5) \;+\; \cdots \]

En comparación, los estimadores de regla de punto medio y regla de trapeio para la integral son

\[\mathcal{I}_n^{\mathrm{mp}} = 2f(x_n) \Delta x\]

\[\mathcal{I}_n^{\mathrm{trapz}} = 2 f(x_n) \Delta x + \frac{f''(x_n)}{2} \Delta x^3 + O(\Delta x^5).\]

Por lo tanto, podríamos tomar el siguiente promedio ponderado:

\[\mathcal{I}_n^{\mathrm{simp}} = \frac{1}{3}\mathcal{I}_n^{\mathrm{mp}} + \frac{2}{3}\mathcal{I}_n^{\mathrm{trapz}} = 2 f(x_n) \Delta x + \frac{f''(x_n)}{3} \Delta x^3 + O(\Delta x^5).\]

Tal promedio ponderado coincidiría con el resultado de la serie Taylor hasta\(O(\Delta x^4)\)! (Puedes verificar por ti mismo que los\(O(\Delta x^5)\) términos difieren.) En resumen, la regla de Simpson para este conjunto de tres puntos puede escribirse como

\[\begin{align} \mathcal{I}_n^{\mathrm{simp}} &= \frac{1}{3} \Big[2f(x_n) \Delta x\Big] + \frac{2}{3} \, \Delta x \Big[ \frac{f(x_{n-1})}{2} + f(x_n) + \frac{f(x_{n+1})}{2}\Big]\\ &= \frac{\Delta x}{3} \Big[f(x_{n-1}) + 4f(x_n) + f(x_{n+1})\Big]. \end{align}\]

El error numérico total, sobre un conjunto de\(O(N)\) segmentos, es\(O(1/N^4)\). ¡Eso es una mejora de dos poderes de\(1/N\) sobre el trapzezio y reglas de punto medio! Lo que es aún mejor es que implica exactamente el mismo número de operaciones aritméticas que la regla del trapecio. Esto es lo más cercano a un almuerzo gratis como se puede conseguir en la ciencia computacional.

9.3.1 Implementación Python de la regla de Simpson

En Scipy, la regla de Simpson es implementada por la función scipy.integrate.simps, que se define en el submódulo scipy.integrate. Similar a la función trapz, esta puede llamarse como simps (y, x) o simps (y, dx=s) para estimar la integral\(\int y \;dx\), utilizando los elementos de x como puntos de discretización, especificando y el conjunto de valores para el integrando.

Debido a que la regla de Simpson requiere dividir los segmentos en pares, si especifica un número par de puntos de discretización en x (es decir, un número impar de segmentos), la función se ocupa de esto haciendo una estimación de regla de trapecio en el primer y último segmento. Por lo general, el error es insignificante, así que no te preocupes por este detalle

Aquí hay un ejemplo de simps en acción:

>>> from scipy import *
>>> from scipy.integrate import simps
>>> x = linspace(0,10,25)
>>> y = exp(-x)
>>> t = simps(y,x)
>>> print(t)
1.00011864276

Para el mismo número de puntos de discretización, la regla del trapecio da\(1.01438\); el resultado exacto es\(0.9999546\dots\) Claramente, la regla de Simpson es más precisa.