Regra de Simpson

Na regra do trapézio o integrando é aproximado por um polinômio de grau 1. Vimos na aula passada que essa aproximação leva a um método efetivo para a estimativa de integrais definidas, porém com a eventual necessidade de repartir o intervalo de integração em muitos subintervalos. O caminho natural para conseguir uma fórmula de integração mais acurada é melhorar a aproximação do integrando. O método de Simpson faz isso aproximando o integrando por um polinômio de grau 2. Para a determinação de um polinômio interpolador de grau 2, devemos prescrever 3 pontos de interpolação. Na figura a seguir ilustra o polinômio utilizado na regra de Simpson, interpolando a função nos extremos do intervalo e também no ponto médio.

A determinação do polinômio $p_2$ pode ser feita através da representação em termos de polinômios de Lagrange, diretamente como $$p_2(x) = f_a {(x-m)(x-b)\over(a-m)(a-b)}+ f_m {(x-a)(x-b)\over(m-a)(m-b)}+ f_b{(x-a)(x-m)\over(b-a)(b-m)}, $$ onde $f_a = f(a)$, $f_m = f(m)$ e $f_b=f(b)$. Definindo $h = (m-a) = (b-m)$, $p_2$ pode ser reescrito como \begin{equation}\label{p2} p_2(x) = f_a {(x-m)(x-b)\over 2h^2}- f_m {(x-a)(x-b)\over h^2}+ f_b{(x-a)(x-m)\over 2h^2}. \end{equation} A regra de Simpson surge de aproximar a integral de $f$ pela integral de $p_2$. Com efeito, $$ I \equiv \int_a^b f(x) \, dx \approx \int_a^b p_2(x)\, dx = 2h\int_0^{1} p_2(2hu+a)\, du, $$ onde usamos a mudança de variável $u = (x-a)/(2h)$. Ao utilizar essa mudança, o integrando se escreve em termos apenas de $u$ e $h$ como \begin{align*} p_2(2hu+a) & = f_a {(2hu-h)(2hu-2h)\over 2h^2} - f_m {2hu(2hu-2h)\over h^2} + f_b{2hu(2hu-h)\over 2h^2} \\ & ={f_a\over 2} (2u-1)(2u-2) -f_m 2u(2u-2) + {f_b\over 2}2u(2u-1). \end{align*} Logo \begin{align*} \int_a^b f(x)\, dx & \approx 2h \left[f_a\int_0^1 (2u-1)(u-1)\, du -4f_m \int_0^1 u(u-1)\,du + f_b \int_0^1 u(2u-1)\, du \right]\\ & = 2h\left[{1\over 6}f_a + {4\over 6}f_m + {1\over 6}f_b\right]\\ & = {h\over 3} \left(f_a + 4f_m + f_b\right). \end{align*} Desta forma, definimos a fórmula de quadratura de Simpson como $$ Q_S[f] = {h\over 3} \left(f_a + 4f_m + f_b\right). $$

Assim como fizemos para a regra do trapézio composta, uma maneira de melhorar a qualidade da integração numérica é aplicar repetidamente a regra acima a subintervalos do intervalo de integração. Como a regra de Simpson simples já trabalha com dois subintervalos (de $a$ a $m$ e de $m$ a $b$), na regra composta precisamos obrigatoriamente trabalhar com uma quantidade par de subintervalos, para poder usá-lo dois a dois na aplicação da regra simples. Considere a situação ilustrada na figura a seguir.

Nesse exemplo, o intervalo original foi fracionado em 4 subintervalos regulares. Logo, a aplicação repetida da regra de Simpson simples leva a $$ I \approx {h\over 3} \big[ (f_0 + 4f_1 + f_2) + (f_2+4f_3+f_4)\big], $$ onde $f_j = f(x_j)$ e $h$ é comprimento dos subintervalos. Perceba que os pontos internos a cada dupla de subintervalos ficam com peso $4$, enquanto que pontos fronteiriços a duplas de subintervalos ficam com peso $2$. De modo geral, a regra de Simpson composta é dada por \begin{equation} Q_{SC}[f] = {h\over 3} \big[ f_0 + 4f_1 + 2f_2 +4f_3+ \cdots + 2f_{n-2} + 4f_{n-1}+ f_n \big ] \end{equation} ou ainda \begin{equation} Q_{SC}[f] = {h\over 3} \left[ f_0 + 4 \sum_{j=1}^{n/2} f_{2j-1} + 2 \sum_{j=1}^{n/2 - 1} f_{2j} + f_n\right], \end{equation} onde $h = (b-a)/n$ e $n$ é par.

Utilizar essa fórmula é tão simples quanto usar a fórmula do trapézio. Vejamos como fazer isto no Octave, usando novamente como exemplo $\int_1^4 (1/x)\exp(x/2)\, dx.$

> f = @(x) 1./x .* exp(x/2);
> n = 20;                     # quantidade par de subintervalos
> x = linspace (1,4,n+1);     # malha de pontos
> h = x(2)-x(1)               # espaçamento entre pontos consecutivos
h =  0.15000
> y = f(x);
>  Q = (h/3) * (y(1) + 4*sum(y(2:2:end-1)) + 2*sum(y(3:2:end-2)) + y(end))
Q =  4.5000
> R = quad(f,1,4);
> abs(R-Q)
ans =  0.000016145

Usando apenas 21 pontos conseguimos um erro absoluto, quando comparado à estimativa da rotina quad do Octave, da ordem de $10^{-5}$. Esta estimativa é melhor que a obtida com a regra do trapézio quando foram usados 101 pontos (veja aula sobre trapézios).

Com este exemplo parece que podemos afirmar que conseguimos atender à expectativa inicial de construir uma fórmula mais acurada que a regra do trapézio composta.

Claro que um único exemplo não é suficiente e uma análise do erro de integração é necessária. Por ora, fique tranquilo, a regra de Simpson é melhor mesmo! Isso significa que devemos sempre usar Simpson ao invés de usar a regra do trapézio? Bom, nem sempre. A regra de Simpson tem duas restrições importantes que sua rival não tem. Em primeiro lugar, os nós de integração devem ser igualmente espaçados, pois a fórmula \eqref{p2} foi derivada com essa hipótese. Isso nem sempre é possível, principalmente quando apenas amostras de $f$ são conhecidas, possivelmente vindas de algum experimento. Em segundo lugar, a quantidade de subintervalos, $n$, deve ser par. Se alguma dessas hipóteses falhar, resta-nos usar a regra do trapézio.

Na próxima aula estudaremos o erro de integração numérica na regra de Simpson.

1. Utilize a regra de Simpson composta para estimar $$ \int_0^{0.5} {1\over x}\mbox{arcsinh} x\, dx. $$ Experimente com diferentes valores de $n$ (quantidade de subintervalos) e determine qual o valor mínimo para o qual a diferença entre sua estimativa e a estimativa produzida pelo comando quad do Octave fica abaixo de $10^{-6}.$ Compare seu resultado de agora com o obtido nesse mesmo exercício quando feito com a regra do trapézio composta.

Dica: $\displaystyle{\lim_{x\to 0} {1\over x}\mbox{arcsinh} x =1}.$