alta precisión inesperada de la regla de Simpson

Question

He intentado utilizar la regla de Simpson para calcular la siguiente integral

$$ \pi = 4 \int_0^1 \frac{dx}{1+x^2} .$$

El código Matlab es el siguiente. Esperaba una ley de escala $error \propto h^4 $

Nlist = 2.^(2:6);
elist = 0 * Nlist;
for s1 = 1 : length(Nlist)
    N = Nlist(s1);
    h = 1/N;
    xlist = 0 : h/2 : 1;
    ylist = 1 ./ (1 + xlist.^2 );

    sum1 =4*h * ( (ylist(1) + ylist(end))*(1/6)...
        + sum(ylist(2:2: end-1))*(4/6)...
        + sum(ylist(3:2 : end-2))* (2/6));

    error = abs(sum1 - pi);
    elist(s1) = error;
end

plot(log2(Nlist), log2(elist),'*')
xlabel('log2(N)')
ylabel('log2(error)')

La cifra generada es

La pendiente es de $-6$ . Esperaba que fuera $-4$ . La regla de Simpson funciona con más precisión de lo esperado. Estoy totalmente desconcertado. ¿Se debe a alguna propiedad oculta del integrante?

Answer 1

Consideremos el problema de calcular la integral $$T = \int_a^b f(x)dx $$ utilizando tanto la regla trapezoidal compuesta $T_h$ así como la regla de Simpson $S_h$ . Si $f \in C^{(2k+1)}([a,b],\mathbb{R})$ entonces el error satisface una expansión asintótica del error de la forma $$T - T_h = \sum_{j=1}^k \alpha_j h^{2j} + O(h^{2k+1}), \quad h \rightarrow 0, \quad h > 0.$$ El valor exacto de los coeficientes de expansión viene dado por $$\alpha_j = \frac{B_{2j}}{(2j)!}(f^{(2j-1)}(b) - f^{(2j-1)}(a))$$ donde $B_{j}$ es el j-ésimo número de Bernoulli. Sólo necesitamos los primeros términos, a saber $$B_2 = \frac{1}{6}, \quad B_4 = \frac{1}{30}, \quad B_6 = \frac{1}{42}.$$ En el caso que nos ocupa, tenemos $$f(x) = (1+x^2)^{-1}.$$ Es tedioso pero sencillo comprobar que $$\alpha_1 \not = 0, \quad \alpha_2 = 0, \quad \alpha_3 \not = 0.$$ Por lo tanto, el error de la regla trapezoidal es $O(h^2)$ . La regla de Simpson puede obtenerse a partir de la regla trapezoidal utilizando la extrapolación de Richardson. En concreto, tenemos $$S_h = T_h + \frac{T_h - T_{2h}}{3}.$$ La extrapolación de Richardson elimina el término de error primario de la expansión de error y expone el siguiente término. En este caso es $O(h^6)$ porque nunca hubo un término que fuera $O(h^4)$ .

Answer 2

El término de error de orden principal de una cuadratura de integración $Q$ en un intervalo $[a, b]$ que integra monomios hasta el grado $n$ es exactamente de la forma $C h^{n + 1}(f^{(n)}(b) - f^{(n)}(a))$ donde $C$ es una constante independiente de $f$ y $h$ . Esto se demuestra fácilmente de la siguiente manera.

Consideremos primero la aplicación de la regla en un intervalo pequeño $[-h, h]$ . Defina $E(f) = \int_{-h}^{h}f(x)\,dx - Q(f)$ . Obsérvese que, dado que la integración y $Q$ son lineales, $E$ es lineal. Ampliar $f$ como una serie de Taylor en torno a $0$ : $$f(x) = \sum_{j = 0}^{n}\frac{f^{(j)}(0)}{j!}x^j + \frac{f^{(n + 1)}(0)}{(n + 1)!}x^{n + 1} + \frac{f^{(n + 2)}(\theta(x)x)}{(n + 2)!}x^{n + 2}.$$ Solicitar $E$ a ambas partes, señalando que los términos $E(x^j) = 0$ para $j \leq n$ para obtener $$E(f) = \frac{f^{(n + 1)}(0)}{(n + 1)!}E(x^{n + 1}) + E(\frac{f^{(n + 2)}(\theta(x)x)}{(n + 2)!}x^{n + 2}).$$ Tenga en cuenta que $E(x^j) = O(h^{j + 1})$ . Del mismo modo, mientras $f^{(n + 2)}$ está acotado (lo que ocurre si $f \in C^{n + 2}([a, b])$ ), entonces $E(\frac{f^{(n + 2)}(\theta(x)x)}{(n + 2)!}x^{n + 2}) = O(h^{n + 3})$ . Por lo tanto, para encontrar el término de error principal, podemos eliminar el término $E(\frac{f^{(n + 2)}(\theta(x)x)}{(n + 2)!}x^{n + 2})$ . Así $$E(f) \approx \frac{f^{(n + 1)}(0)}{(n + 1)!}E(x^{n + 1}) = Cf^{(n + 1)}(0)h^{n + 2}.$$ Ahora el error $e(h)$ sobre el método compuesto en $[a, b]$ con subintervalos de longitud $2h$ sumamos los errores en cada subintervalo $[c_j - h, c_j + h]$ hasta conseguir \begin{align} e(h) &\approx \sum_{j = 0}^{N - 1}Cf^{(n + 1)}(c_j)(2h)^{n + 2} \\ &= C(2h)^{n + 1}\sum_{j = 0}^{N - 1}f^{(n + 1)}(c_j)(2h) \\ &\approx C(2h)^{n + 1}\int_{a}^{b}f^{(n + 1)}(x)\,dx \\ &= C(2h)^{n + 1}(f^{(n)}(b) - f^{(n)}(a)). \end{align} Nótese que en la tercera igualdad reconocimos la regla del punto medio como una aproximación de $\int_{a}^{b}f^{(n + 1)}(x)\,dx$ .

En su caso, tenemos la regla de Simpson, que tiene $n = 3$ . Con su función $f(x) = \frac{1}{1 + x^2}$ , $f^{(3)}(1) - f^{(3)}(0) = 0$ . La explicación de por qué su pedido fue $6$ es la siguiente. Realizando una demostración similar a la anterior sobre la regla del punto medio (digamos para los dos primeros términos), se puede ver que el error de la regla del punto medio $e_m(h)$ tiene una "expansión asintótica" como $$e_m(h) = C_1h^2(f'(b) - f'(a)) + C_2h^4(f^{(3)}(b) - f^{(3)}(a)) + \dots.$$ Con esto en la mano, puede llevar a cabo la prueba anterior para la regla de Simpson más (digamos que a los dos primeros términos. Tenga en cuenta que $E(x^j) = 0$ para impar $j$ debido a la simetría), y verá que el error asintótico es de la forma ( $C_1$ y $C_2$ aquí son diferentes a los de la regla del punto medio) $$e_s(h) = C_1h^4(f^{(3)}(b) - f^{(3)}(a)) + C_2h^6(f^{(5)}(b) - f^{(5)}(a)) + \dots$$ Los coeficientes exactos $C_j$ en la expansión del punto medio puede escribirse en términos de números de Bernoulli (fórmula de Euler Maclaurin).