TÉCNICAS DE INTEGRAÇÃO NUMÉRICA

Prof. Luiz Eduardo Ourique
Pontifícia Universidade Católica do Rio Grande do Sul – PUCRS
ourique@pucrs.br
Rodrigo Mello Nardi
Pontifícia Universidade Católica do Rio Grande do Sul – PUCRS
rodrigomellonardi@gmail.com 

Introdução
Nos cursos de graduação da área de Ciências Exatas das universidades brasileiras, a formação geral dos alunos inclui as disciplinas de Cálculo Diferencial e Integral, onde são estudadas diversas aplicações de derivadas, conforme a apresentação clássica dos livros-textos: estudo de crescimento e decrescimento de funções, localização de máximos e mínimos locais e globais de funções de uma ou mais variáveis, problemas de taxas de variação, entre outros. As aplicações das integrais incluem cálculo de áreas de regiões planas, volume de sólidos de revolução, momentos de inércia, entre outros; em geral, a formação inclui o estudo de propriedades, como a linearidade da derivada e da integral, até chegar a modelagem de problemas simples e a solução analítica dos mesmos. Os livros mais antigos apresentam problemas resolvidos de modo exato, como podemos ver em Courant (1970). Os livros mais modernos apresentam os problemas clássicos do Cálculo Diferencial e Integral citados acima, além de formularem questões a serem resolvidas através do uso de recursos computacionais, em particular, os CAS , abreviatura em inglês para “sistemas algébricos computacionais”, ou programas de computador que efetuam cálculos matemáticos numérica e simbolicamente, como podemos ver em Anton (2007)  e Stewart (2006).  Como exemplos de sistemas algébricos computacionais, podemos citar o Derive, o Maple e o Mathematica. 
Alguns cursos de graduação da área de Ciências Exatas incluem a disciplina de Cálculo Numérico como obrigatória:  os alunos então estudam algumas técnicas para a resolução numérica de problemas, entre os quais o cálculo aproximado de raízes de equações algébricas e transcendentes, de derivadas e integrais. Em geral, as técnicas de integração numérica tradicionais incluem as chamadas fórmulas de Newton-Cotes, como a fórmula dos trapézios e a fórmula de Simpson. Num universo mais restrito, alguns cursos incluem as quadraturas de Gauss-Legendre em seus programas. Estas quadraturas apresentam, em alguns casos, vantagens em relação às fórmulas de Newton-Cotes, não obstante, a dedução das fórmulas envolvidas e algumas de suas propriedades serem mais complexas, uma vez que são obtidas a partir de propriedades dos polinômios ortogonais de Legendre. A utilização de técnicas numéricas de aproximação envolve desde a dedução das fórmulas a partir de alguns princípios, até a geração de aproximações sucessivas, levando a questões como convergência de seqüências (a seqüência converge? Em caso afirmativo, converge para a solução exata do problema?) , determinação do erro, existência e unicidade de solução, que são questões discutidas tanto no Cálculo como na Análise Matemática. 
O objetivo deste trabalho é apresentar os resultados de uma comparação entre as técnicas de integração numérica tradicionais  e as quadraturas de Gauss-Legendre, numa aplicação prática, qual seja de calcular a eficiência luminosa, conforme o problema proposto por Dorn e McCracken (1981). A eficiência luminosa, escrita na forma de uma integral definida, foi calculada por diferentes quadraturas (fórmula dos trapézios, Simpson e de Gauss-Legendre) e foi feita a comparação do resultado, em termos do  esforço computacional. Os cálculos foram efetuados com o software Maple. Os resultados foram apresentados no Salão de Iniciação Científica da PUCRS, em outubro de 2006, na forma de pôster e comunicação oral, por parte do aluno orientando. A pesquisa foi realizada ao longo de todo o ano de 2006, o que incluiu levantamento bibliográfico dos títulos disponíveis na biblioteca da universidade, dedução das fórmulas de quadratura, estudo do erro, extrapolação pelo limite e implementação computacional, usando o Maple, na resolução de uma aplicação encontrada na literatura. Neste trabalho, vamos apresentar os resultados desta pesquisa, conforme o seu andamento. Inicialmente, são deduzidas as fórmulas de integração baseadas em pontos uniformemente distribuídos, como a fórmula dos trapézios e a fórmula de Simpson; depois, são apresentadas  as quadraturas de Gauss-Legendre. Na parte final, foi feita uma aplicação das fórmulas no cálculo da eficiência luminosa e a comparação dos resultados e uma discussão das práticas de ensino de Cálculo Numérico.
Quadraturas ou fórmulas de integração numérica 
Consideremos o problema de calcular a integral definida 
                                                                                                 (1) onde a e b são finitos e f (x) é uma função contínua de x  para   a  x  b. Para calcular uma aproximação desta integral através de uma quadratura, calculamos a função em pontos ( nós ) e efetuamos uma combinação linear das imagens da função nos pontos. Isto é, podemos escrever
                    = + Erro                      (2)
onde  Ak é o coeficiente da quadratura correspondente ao ponto xk. A soma  é chamada de quadratura. Isto é, quadratura é o nome geral dado às técnicas de integração numérica. A fórmula dos trapézios e a fórmula de Simpson são casos particulares de quadraturas. Estas são divididas em duas classes, conforme a localização dos pontos e os valores dos coeficientes : 
1. Numa classe os pontos onde a função é avaliada são uniformemente distribuídos no intervalo [a, b] ; em particular, nas fórmulas dos trapézios e de Simpson, a função é desenvolvida em séries de Taylor  e os coeficientes da quadratura são determinados de modo que a fórmula seja exata para polinômios de grau 1 e 3, respectivamente.
2. Em outra, os pontos são as raízes de polinômios especiais, os polinômios de Legendre e determinados de modo que a fórmula satisfaça um critério de exatidão pré-definido; os coeficientes são então calculados para que a fórmula satisfaça tal critério, como ocorre nas quadraturas de Gauss-Legendre, ou mais simplesmente, quadraturas gaussianas. Pode ser demonstrado que uma quadratura gaussiana com n pontos resulta em valores exatos para polinômios de grau menor ou igual a 2n – 1. 
As quadraturas gaussianas e fórmulas dos trapézios e de Simpson,  representantes das classes citadas acima,  são as mais  utilizadas atualmente nos mais variados problemas e aplicações que exigem o cálculo numérico de integrais. Como aplicações práticas, podemos citar o cálculo do comprimento de uma curva ou do volume de uma mina de carvão. Em problemas de Matemática Pura e Aplicada, podemos citar algumas funções matemáticas mais avançadas , como a função erro de Gauss a as funções de Fresnel. O estudo do erro destas fórmulas é bastante conhecido; em alguns casos, é possível estimar quantos pontos são necessários para que a quadratura satisfaça um critério de exatidão pré-fixado. Outras fórmulas de quadraturas são conhecidas e têm aplicações específicas em problemas nos quais uma propriedade particular pode tornar uma determinada fórmula mais atraente, como simetria, a localização e o  número de pontos. 
A regra do trapézio 
Consideremos a integral (1); dividimos o intervalo [a,b] em n intervalos iguais, cada um de comprimento igual a h = ( b – a ) / n, com pontos dados por xi= a + ih, i=0,...,n. Decorre da definição de h que x0=a  e b = xn . O comprimento h é chamado também de incremento da divisão, pois xi+1– xi = h. Faremos uma dedução intuitiva da fórmula. A área limitada pelo gráfico de y = f(x) e o eixo x , entre xi e xi+1, é igual a  . O  polinômio de primeiro grau que interpola f nos pontos (xi, yi ) e (xi+1, yi+1) é  p(x) = yi + (x-xi)(yi+1- yi)/ (xi+1-xi) . Se h é suficientemente pequeno, a integral da função é bem aproximada pela integral do polinômio interpolador. Temos, então,  = , aplicando diretamente o Teorema Fundamental do Cálculo, onde o símbolo  denota aproximadamente igual e yi=f(xi). Devido a linearidade da integral definida, isto é, , podemos escrever 
                                                                                                         (3)
Então, da linearidade e da igualdade (3), segue a fórmula dos trapézios: 
                         Ih =                               (4)
onde  Ih denota a fórmula dos trapézios, correspondente ao incremento h, assim chamada porque, se f(x)  0 , no intervalo [ a, b ] , aproxima a área sob uma curva pela soma da área de n trapézios.  A simplicidade desta fórmula é um dos fatores que a tornam atraente e utilizada em muitos problemas, apesar do erro de truncamento ser maior do que em outras fórmulas. O erro de truncamento é a soma das diferenças entre o valor exato da integral e a área do trapézio em cada intervalo da divisão. Usando o desenvolvimento da função f(x) em série de Taylor nos pontos extremos da divisão do intervalo [a,b], pode ser demonstrado que o erro de truncamento é proporcional a derivada segunda de f(x). Logo, o erro cometido pela fórmula dos trapézios é igual a zero, se a função integrando for um polinômio de grau menor ou igual a um, pois neste caso o trapézio coincide com o gráfico da função. Para estimar o erro de truncamento máximo, devemos estimar o valor máximo da segunda derivada de f(x) no intervalo de integração, o que nem sempre é trivial, podendo exigir algum esforço computacional. 
A extrapolação para o limite 
Para melhorar a aproximação obtida pela fórmula dos trapézios, podemos efetuar uma modificação baseada na idéia de aproximação sucessiva, com incrementos cada vez menores. De acordo com a fórmula do erro, se h = ( b–a ) / n, o erro de truncamento é proporcional a h2, ou seja, podemos escrever et = Ch2 ,onde esta constante C depende da segunda derivada de f(x), calculada em algum ponto do intervalo [a,b].     Suponhamos que a segunda derivada de f seja aproximadamente constante e tomemos um incremento diferente  k = (b – a) /m, onde m n. Então et = Ck2, Sendo I o valor exato da integral, Ih o resultado da fórmula dos trapézios com incremento h e Ik o resultado com o incremento k, podemos escrever:
                      I= Ih + Ch2            e              I = Ik + Ck2                                       (5)   Se subtrairmos estas duas equações, obtemos C =  . Substituindo esta fórmula, na primeira equação em (5), obtemos 
                                     I = Ih +                                                            (6) A fórmula (6) resulta numa melhor aproximação para I do que Ih ou Ik. Em geral, partindo de um valor de h e usando que acréscimos cada vez menores, por exemplo, pela divisão sucessiva de h por potências de 2, são geradas aproximações sucessivas de I, com erros de truncamento cada vez menores, desconsiderando-se o erro de arredondamento. O método acima é chamado de método de extrapolação para o limite de Richardson. 
A regra de Simpson
Considerando a integral (1), na fórmula dos trapézios definimos que o número de intervalos da divisão do intervalo de integração é igual a n = ( b – a ) / h. Suponhamos que n seja par; definimos então 
                                                  k = 2h                                                              (7)   e então aplicando a fórmula dos trapézios correspondente a estes dois incrementos:
             Ih=                          (8)
             Ik=                               (9)
De fato, na equação (9), o número de pontos é  a metade dos que temos em (8), pois o incremento k é o dobro de h. As equações (8) e (9) podem ser substituídas na equação (6) . Observando que k2 = 4h2 , obtemos após algumas simplificações a igualdade
               I =                     (10)
A equação (10) é chamada de fórmula de Simpson. Pode ser obtida a partir da interpolação por um polinômio de segundo grau que coincide com a função integrando em três pontos adjacentes, embora o esforço seja maior. Pode ser mostrado que o erro de truncamento é proporcional a quarta derivada de f(x), calculada em algum ponto do intervalo [ a, b]. Isto significa que para uma função bem comportada, o erro na fórmula de Simpson é menor do que na fórmula dos trapézios. Comparando as fórmulas para o erro, vemos que na fórmula de Simpson o erro é proporcional a h4, enquanto para a regra do trapézio o erro é proporcional a h2. Isto quer dizer que a regra de Simpson corresponde aos primeiros três termos da série de Taylor, logo esta fórmula é exata para polinômios de grau menor ou igual a três. Admitindo que a quarta derivada de f(x) é constante, podemos usar a extrapolação pelo limite para melhorar a fórmula de Simpson. De fato, supondo que 
                      I= Ih + Ch4            e              I = Ik + Ck4                                                 (11) decorre imediatamente que a fórmula de extrapolação de Richardson para a regra de Simpson é
                                                      I = Ih +                                                     (12)
Quadraturas de Gauss-Legendre
As quadraturas de Gauss-Legendre ou quadraturas gaussianas são fórmulas com a forma (2), porém são mais precisas que as anteriores, quando comparamos os resultados obtidos com o mesmo número de pontos. Nestas quadraturas, os pontos não são uniformemente distribuídos ao longo do intervalo de integração e são determinados, assim como os coeficientes, de modo que a quadratura tenha uma exatidão pré-determinada. A determinação da fórmula decorre de algumas propriedades dos polinômios de Legendre. Para o leitor interessado, citamos o livro de Stroud(1974), onde são deduzidas as propriedades das quadraturas de Legendre, que são ortogonais no intervalo [-1,1]. Uma propriedade fundamental é que estas fórmulas são exatas para polinômios de grau menor ou igual a  2n – 1 , onde n é o número de pontos da quadratura. Vamos apresentar aqui somente a dedução da fórmula para o caso n = 2. O problema continua sendo o de calcular  . Inicialmente, o intervalo de integração deve ser transformado no intervalo [ – 1, 1 ] , mediante a troca de variável : 
                         x  =                                                 (13)
donde
                                      =                                                       (14)
                                                      (15)
Para dois pontos, a fórmula de Gauss é :
                        I =  =                                  (16) onde A0 , A1, t0 e t1 são incógnitas a se determinar e independentes da função F(t). Como esta fórmula deve ser exata para polinômios de grau menor ou igual a 3, devemos considerar as funções  F(t) = tk , k =0, 1,2, 3.  Assim: 
                                                                             (17)
Teremos : 
k = 0 
k = 1 
k = 2 
k = 3 
Calculando as integrais, obtemos o sistema de equações não lineares: 
                                                                               (18)             Resolvendo o sistema (18), obtemos  A0 = A1  = 1,  t0 = – t1 = . Substituindo os valores encontrados na equação (16), temos a fórmula de Gauss para dois pontos: 
                              Ig = F( – )  +  F ( )                                           (19) Lembremos que esta fórmula é exata para polinômios de até terceiro grau. O erro de truncamento satisfaz  a fórmula et =  , onde  – 1  u  1. A fórmula geral para a quadratura gaussiana, baseada em propriedades dos polinômios de Legendre, é              
                                       I  =  =                                      (20) onde n é o número de pontos da quadratura, Ai são os coeficientes e os pontos ti  são as raízes dos polinômios de Legendre de grau n. Os primeiros polinômios de Legendre podem ser definidos recursivamente 
P0(t) = 1,  
P1(t) = t ,
 Pm(t) =  , se m  2. Por exemplo, fazendo m = 2, resulta P2(t) =  = 3t2/2 – 1/2  , cujas raízes são , como encontrado no cálculo acima. A fórmula para os pesos é 
                                                                          (21)  Como exemplo, considere n = 2, de modo que t0 = –  ,  (t) = 3t. Então (t0) = 3t0 = –  , de modo que  A0 =  = 1 , e A1 = 1, como acima.  Uma fórmula para o erro de truncamento é et = , onde u . Lembremos que todas estas propriedades decorrem da definição da quadratura gaussiana. No quadro 1, são dados os pontos e os coeficientes da quadratura gaussiana, para diferentes valores de n. 

	T	A	
n=2	0,57735 02692	1,00000 00000	
n=3	0,77459 66692
    0,00000 00000	0,55555 55556
0,88888 88889	
n=4	0,86113 63116
0,33998 10436	0,34785 48451
0,65214 51549	
n=6	0,93246 95142
0,66120 93865
0,23861 91861	0,17132 44924
0,36076 15730
0,46791 39346	
Quadro 1. Pontos e coeficientes da quadratura gaussiana

Quadros mais completos com os pontos e os coeficientes das quadraturas de Gauss para valores mais elevados de n são encontradas na literatura, para facilitar a aplicação das fórmulas, por exemplo em Abramowitz(1988). Além disto, os programas como Maple, Mathematica e outros calculam de modo bastante rápido estes valores, pois tem em suas bibliotecas as funções matemáticas mais avançadas, como os polinômios de Legendre. À parte o esforço para dedução destas fórmulas, devemos ter em mente que num primeiro momento, em particular num curso de graduação, devemos utilizar as fórmulas para mostrar a sua validade. Nos casos mais simples, podemos usar a interpretação geométrica para justificar alguns resultados para não exagerar no formalismo e usar a intuição, como no caso da fórmula dos trapézios. Como uma aplicação, mostremos o cálculo da eficiência luminosa. 
Cálculo da Eficiência Luminosa
Um corpo negro (radiador perfeito) emite energia em uma taxa proporcional à quarta potência de sua temperatura absoluta T ( medida em Kelvin) , de acordo com a equação de Stefan-Boltzmann. Uma parte desta energia forma o espectro visível, correspondendo ao comprimento de onde entre   a   cm. A eficiência luminosa é definida como a relação da energia no espectro visível para a energia total. Multiplicando por 100 para obter esta energia na forma percentual, o problema equivale a calcular 
                              EFF =                          (22)  
A notação EFF denota a eficiência luminosa, medida em W/cm2, x é o comprimento de onda, exp(x) denota a função exponencial ex. . A integral em (22) não pode ser calculada pelas técnicas usuais de integração ( substituição, integração por partes e outras ) , isto é, a função integrando não possui primitiva em termos das funções elementares do Cálculo. Assim, para calcular a eficiência luminosa , devemos calcular numericamente a integral. Para efeito de comparação, nós calculamos a eficiência luminosa, pelas três fórmulas apresentadas acima. Os resultados estão no quadro 2, supondo uma temperatura T=3500 K : 

N	Trapézios	Simpson	
10	14.48892779	14.51272461	
20	14.49496856	14.51266669 	
Quadro 2. Eficiência luminosa 

Usando a técnica da extrapolação pelo limite, obtemos os resultados do quadro 3.

Trapézios	Simpson	
14.49698215	14.51266283	
Quadro 3. Extrapolação pelo limite

Isto quer dizer que, usando a extrapolação pelo limite pela fórmula de Simpson, a eficiência luminosa é de 14,51266 %. 
Aplicando a quadratura gaussiana, obtemos os valores do quadro 4:

n	EFF	
2	14.46748120	
4	14.51267557	
Quadro 4. Quadraturas gaussianas 

A aplicação da quadratura de Gauss com 4 pontos nos diz que a eficiência é 14,51267%, um resultado que é equivalente ao obtido pela fórmula de Simpson com 20 pontos. Este exemplo mostra que a quadratura de Gauss é bem mais eficiente para o cálculo da integral proposta, uma vez que para obter um resultado com aproximadamente quatro casas decimais corretas, foram necessários apenas quatro pontos, contra vinte da fórmula de Simpson. Embora um pouco fora do escopo deste trabalho, observemos que no ponto de fusão do Tungstênio, cerca de 3600 K, somente 15% da energia irradiada é visível. Na aplicação de todas as fórmulas acima, os cálculos foram efetuados usando o Maple, que é um programa disponível aos alunos dos cursos de engenharia da PUCRS.
Conclusões
As diferentes técnicas de integração discutidas aqui podem ser utilizadas para o cálculo de integrais definidas. As fórmulas dos trapézios e de Simpson são mais práticas, pois os pontos são uniformemente distribuídos, obtidos a partir dos dados da integral. Nas aplicações onde os valores de uma função são obtidos experimentalmente (por exemplo, o volume de terra a ser extraído na construção de um túnel), estas fórmulas são bastante úteis. Por outro lado, em problemas onde se procura um resultado mais acurado, com menor esforço computacional, a quadratura gaussiana é bastante vantajosa, pois exige o cálculo da função em um menor número de pontos. De modo geral, o esforço computacional para implementar a quadratura gaussiana é a metade do esforço correspondente à fórmula de Simpson, para obter um resultado com erros equivalentes. Da mesma forma, a fórmula de Simpson requer a metade do esforço da regra do trapézio, para obter um resultado equivalente. A conclusão deste trabalho levou o aluno orientando a se perguntar, em outros problemas, qual método é o mais eficiente para resolver um dado problema.  O uso do Maple permitiu que os cálculos, bastante complexos pela natureza da função envolvida, fossem efetuados de modo simples, bastando para isto, saber definir a função na linguagem do programa e criar as variáveis necessárias para armazenar os resultados intermediários, o que exige um certo conhecimento de vetores. Algumas iniciativas do uso deste recurso tem sido reportadas, por exemplo, o estudo de séries de Fourier e o fenômeno de Gibbs, que são assuntos mais avançados, como em Alves (1997). Isto mostra que o Cálculo Numérico é uma disciplina que pode passar por modificações, em função do uso de recursos computacionais, pois problemas mais complexos podem ser resolvidos, comparados com os cursos “pré-computador”, como é citado no trabalho de Matta et al. (1997). Deste modo, conteúdos considerados mais avançados podem vir a ser abordados, pois o tempo usado em cálculos numéricos com máquinas lentas e com poucos recursos pode ser substituído pelo estudo de problemas mais complexos, deixando os cálculos a cargo do computador. O esforço computacional é minimizado com o uso de programas mais sofisticados, de modo que o tempo pode ser usado para o embasamento teórico, a  análise dos resultados e escolha de métodos mais eficientes. A integração numérica  é um caso típico de problema que hoje pode ser mais explorado pelos professores de Cálculo Numérico, pois com o uso de recursos computacionais mais sofisticados as quadraturas gaussianas podem ser usadas mais cedo na graduação, possibilitando aos estudantes um contato com uma técnica mais eficiente de integração. É claro que o uso de programas como o Maple e outros sistemas de computação algébrica exigem treinamento dos usuários, pois estes sistemas geralmente são concebidos de resolver desde os problemas mais elementares de Matemática até os mais avançados. Assim, os estudantes e professores que usam estes programas precisam estar constantemente a par das modificações que ocorrem de uma versão para outra, devem estar atentos para eventuais “bugs” e, principalmente, saber escrever os algoritmos de maneira a minimizar o tempo de processamento e o consumo de energia. 

Referências bibliográficas:
ALVES, G. O Maple na modernização do Cálculo. In: CONGRESSO BRASILEIRO DE ENSINO DE ENGENHARIA, 25., 1997, Salvador. Anais.... Salvador: ABENGE, 1997. v. 2. p. 919- 934.

ANTON , H. et alli. Cálculo: um novo horizonte. Porto Alegre: Bookman, 2007. 

ABRAMOWITZ, M.; STEGUN, I. Handbook of mathematical functions, with formulas, graphs and mathematical tables. New York:  Dover, 1998.

BARROSO, L. C. et alli. Cálculo Numérico com aplicações. São Paulo: Harbra, 1987. 

CLAUDIO, Dalcidio Moraes. Cálculo Numérico Computacional ( teoria e prática ): algoritmos em pseudo-linguagem. São Paulo: Atlas, 2000.

COURANT, R.  Cálculo Diferencial e Integral. Porto Alegre: Editora Globo, 1970. 

DORN, William S.; McCRACKEN, Daniel D. Cálculo numérico com estudos de casos em FORTRAN IV. São Paulo: Editora da USP, 1981.

MATTA, E. N. da.; SOARES, M. V.; LAURIA, D. O curso de Cálculo Numérico na Escola de Engenharia Mauá. In: CONGRESSO BRASILEIRO DE ENSINO DE ENGENHARIA, 25., 1997, Salvador. Anais.... Salvador: ABENGE, 1997. v. 2. p. 935- 950.

STEWART, J. Cálculo. São Paulo: Pioneira Thomson Learning, 2006.
STROUD, H. Numerical Quadrature and Solution of Ordinary Differential Equations: a textbook for a beginning course in numerical analysis. NewYork: Springer, 1974.