Integral usando o método de Romberg e Simpson

---

Algum conhecimento de integração numérica é necessário para esta postagem, sugerimos a esta leitura aqui antes de continuar.

Em alguma situações desejamos calcular a integral de uma função de $a$ até $b$:
\[ \int_a^b f(x) dx = F(a) - F(b) \]
Entretanto, nem sempre é possível ou viável calcular $F(x)$ ou não conhecemos a função $f(x)$, mas temos somente um conjunto de pontos $(x_i,y_i)$. Nesses casos devemos lançar mão de métodos numéricos para calcular a integral. Um dos mais simples é a chamada regra dos trapézios:
\[ \int_a^b f(x) dx \cong \frac{h}{2}(f(a) + f(b)) + E_O(h^2)\]
sendo $h = (b-a)$ e $E_O(.)$ é o erro que é função de $h^2$ e da derivada segunda da função $f(x)$. Se dividirmos o intervalo na metade, o erro tente a ser 4 vezes menor. Vejamos um exemplo simples: qual a integral da função $x^3$ no intervalo $(0,2)$? Numericamente, usando a regra dos trapézios com passo $h=2$:
\[ I_{t2} \cong f(0) + f(2) = 8\]
que apresenta um erro considerável. Se fizermos $h=1$:
\[ I_{t1} \cong \frac{1}{2}(f(0) + 2f(1) + f(2)) = 5\]
o erro é bastante reduzido. Podemos, aproximadamente, cancelar esse erro com a seguinte observação:
\[ 4 \times I_{t1} - I_{t2} \cong 3\int_0^2 x^3 dx \]
ou melhor:
\[ \int_0^2 x^3 dx  \cong \frac{4 \times I_{t1} - I_{t2}}{3}\]
Para os valores calculados acima:
\[ I = \frac{4 \times 5 - 8}{3} = 4\]
que é o valor exato da integral da função $x^3$ no intervalo $(0,2)$.

Em geral, o método de Romberg* utiliza a regra dos trapézios em seu algoritmo.  De forma geral, o algoritmo de Romberg usa a seguinte expressão recursiva:
\[     I_{j,k} \cong \frac{4^{k-1}I_{j+1,k-1} - I_{j,k-1}}{4^{k-1}-1}  \]
onde $I_{j,k}$ é a integral melhorada, $I_{j+1,k-1}$ e $I_{j,k-1}$ são as versões menos acuradas. A integração usando o algoritmo de Romberg exige o conhecimento da função $f(x)$, pois dados tabulados dificilmente teriam os pontos necessários para os refinamentos.

Podemos melhorar ainda mais o método de Romberg pelo uso da primeira regra de Simpson** no lugar da regra do trapézio. O esforço computacional aumenta muito pouco e o erro cai duas ordens de grandeza. A primeira regra de Simpson para o cálculo de uma integral exige 3 pontos:
\[ \int_a^b f(x) dx \cong \frac{h}{3}(f(a) + 4f(a+h) + f(a+2h)) + E_O(h^4)\]
com $a+2h = b$, $h = (b-a)/2$. Notamos que o erro cai 16 vezes quando o passo cai pela metade. Se $I_{s1}$ é uma aproximação da integral usando a primeira regra de Simpson com passo $h/2$ e $I_{s2}$ é uma aproximação da integral usando a primeira regra de Simpson com passo $h$, então:
\[ \int_a^b f(x) dx  \cong \frac{16 \times I_{s1} - I_{s2}}{15}\]
E o algoritmo de Romberg pode usar a seguinte expressão recursiva quando a primeira regra de Simpson é usada para calcular as integrais numéricas:
\[     I_{j,k} \cong \frac{16^{k-1}I_{j+1,k-1} - I_{j,k-1}}{16^{k-1}-1} \]

O código Scilab abaixo implementa o algoritmo Romberg e calcula a integral da função $f(x) = 4.5 + 4\cos(x) - 8 e^{-4x}$ no intervalo $(0,4)$. O passo inicial é $h=1$.

Código:

/////// Romberg
clear; clc; close;

function is=itsp1(xi, yi);
    h = xi(2) - xi(1);
    t = max(size(xi));
    p = ones(1,t);
    p(2:2:t) = 4;
    p(3:2:t-1) = 2;
    is = h*sum(p.*yi)/3;
endfunction

function f=fx(x)   //função
    f = 4.5 + 4*cos(x) - 8*exp(-x*4);
endfunction

function g=fxi(x)   //integral da função
    g = 4.5*x + 4*sin(x) + 2*exp(-x*4);
endfunction

xg = 0:0.001:4;
fg = fx(xg);
title('Função de teste');
plot(xg,fg); xgrid(4);
valor_teorico = fxi(4) - fxi(0);
mprintf('valor_teorico = %1.6f \n',valor_teorico);
xgi = 0:1:4;
fgi = fx(xgi);
plot([0, xgi, 4],[0, fgi, 0],'m');

MR = zeros(5,5);
h = 1;
for k=1:5
    x = 0:h:4;
    y = fx(x);
    MR(k,1) = itsp1(x,y);
    h = h/2;
end;
tol = 1e-9;
df = 1;
C = 1;
while (df>tol)&(C<6)
    C=C+1;
    k4 = 16^(C-1);
    L = 0;
    while (df>tol)&(L<6-C)
        L = L + 1;
        MR(L,C) = (k4*MR(L+1,C-1) - MR(L,C-1))/(k4-1);
        df = abs(MR(L,C) - MR(L+1,C-1));
        m_est = MR(L,C);
    end
end

for L=1:5
    mprintf('%1.6f & %1.6f & %1.6f & %1.6f & %1.6f\n',MR(L,1:5));
end

erro = m_est - valor_teorico;
mprintf('Erro = %e,  *** %1.6f',erro,valor_teorico);

---

Resultados numéricos:

valor_teorico = 12.972790

12.089847 & 12.904267 & 12.970364 & 12.972743 & 12.972789
12.853366 & 12.970106 & 12.972742 & 12.972789 & 0.000000
12.962810 & 12.972732 & 12.972789 & 0.000000 & 0.000000
12.972112 & 12.972789 & 0.000000 & 0.000000 & 0.000000
12.972747 & 0.000000 & 0.000000 & 0.000000 & 0.000000

Erro = -7.631405e-07,  *** 12.972790

---

* Werner Romberg (1909–2003) foi criado em Berlim. Entrou na Universidade de Heidelberg em 1928, onde estudou matemática e ciências físicas. Também estudou em Munique. Devido o nazismo, ele foi obrigado a fazer carreira fora de sua terra natal. Seus interesses incluíram o uso de computadores para computação numérica.

** Thomas Simpson nasceu em 20 de agosto de 1710 em Market Bosworth, Leicestershire, Inglaterra. Morreu em 14 maio de 1761 in Market Bosworth. Thomas recebeu pouca educação formal. Ele frequentou a escola no Market Bosworth por um tempo, mas seu primeiro trabalho foi como um tecelão. Foi autodidata em matemática. Ele trabalhou com probabilidade, mas ficou mais conhecido pelos seus trabalho sobre interpolação e integração numérica.