Integração Numérica

Neste capítulo vamos estudar métodos para aproximar a função primitiva F, resultante de integrar uma função conhecida f. Poderiamos encarar este problema numa perspectiva geral, em que o objectivo seria aproximar uma solução de uma equação diferencial, já que o este caso corresponde a encontrar F tal que F' = f.

Mais concretamente, basta-nos concentrar na integração de uma função f num certo intervalo [a, b]. A ideia base é aproveitar a aproximação por interpolação polinomial, que já estudámos, para obter uma aproximação razoável da função integranda através de polinómios, que são funções fáceis de integrar.
Um exemplo simples é aproximar a função por rectas interpoladoras nos pontos a e b,

claro que, neste caso, a aproximação pode estar desajustada, e podemos melhorar a aproximação, usando, por exemplo, um polinómio de grau superior ou um spline linear.

Dum modo geral, com vista a aproximar o integral
 

I(f ) =  ó
õ
b

a

f(x) dx 

vamos considerar fórmulas de integração (também designadas fórmulas de quadratura) do tipo:
 

Q(f ) =   n
S
i =1
Ai f( xi ), 

onde x0 , ... , xn são pontos conhecidos, pertencentes ao intervalo [a, b], designados por nós de integração,
e os coeficientes A0 , ... , An são coeficientes a determinar, independentes da função f, que designamos pesos.

Um critério habitual, para determinar os pesos na fórmula de quadratura, é exigir que se obtenha o valor correcto do integral de polinómios de um certo grau.

Assim, se a fórmula verificar

In( pm ) = I ( pm )
para qualquer polinómio pm de grau <m, dizemos que a fórmula de quadratura tem, pelo menos, grau m .

Se, para além disso, tivermos

In( pm+1 ) ¹ I ( pm+1 )
para um certo polinómio pm+1 de grau m+1, dizemos que a fórmula de quadratura tem grau m.



Método dos Coeficientes Indeterminados
Como, quer o integral, quer a fórmula de quadratura, são operadores lineares, exigir que a fórmula tenha pelo menos grau m:
In( a0 + ... + amxm ) = I ( a0 + ... + amxm )
equivale a
a0 In( 1 ) + ... + am In ( xm ) = a0 I ( 1 ) + ... + am I ( xm )
para quaisquer a0 , ... , am. Ou seja, equivale a exigir que
In ( 1 ) = I ( 1 )
In ( x ) = I ( x )
... 
In ( xm ) = I ( xm )
substituindo as expressões, obtemos um sistema
ì
ï
ï
ï
ï
ï
ï
í
ï
ï
ï
ï
ï
ï
î
 n
S
i =1
A ó
õ
b

a

1 dx = b-a
 n
S
i =1
Ai f( xi ) =  ó
õ
b

a

x dx = (b2-a2)/2
(...)
 n
S
i =1
Ai xim  ó
õ
b

a

xm dx = (bm+1-am+1)/(m+1)

Para que o sistema seja possível e determinado, impomos que m = n, obtendo-se assim

é
ê
ê
ê
ê
ë
1  1   ...   1
x0  x1  ...   xn
... ...
x0n  x1n  ... xnn
ù
ú
ú
ú
ú
û
é
ê
ê
ê
ê
ë
A0
A1
...
An
ù
ú
ú
ú
ú
û
 =  é
ê
ê
ê
ê
ë
b-a
(b2-a2)/2
...
(bn+1-an+1)/(n+1)
ù
ú
ú
ú
ú
û
onde a matriz do sistema é a transposta da matriz de Vandermonde que, como vimos, é invertível, se os nós forem distintos.
Temos, assim, dados os nós de integração, uma solução única para os pesos A0 , ... , An .
Reparando que o valor da fórmula de quadratura corresponde ao valor do integral do polinómio interpolador, obtemos, através da fórmula de interpolação de Lagrange :
 
Q(f ) = I(pn) =  ó
õ
b

a

 N
S
i =0
f(xi) li(x) dx =   N
S
i =0
 f(xi) ó
õ
b

a

li(x) dx
e assim, os pesos podem também ser dados através dos integrais dos polinómios base de Lagrange li( x ),
Ai ó
õ
b

a

li(x) dx


Regra dos Trapézios

Regra dos Trapézios (simples)

Trata-se do caso mais simples de Fórmula de Newton-Cotes fechada. Neste caso, consideramos n=1 e temos dois nós de integração:

x0 = a, x1 = b

e obtemos para os valores dos pesos:

A0 = A1 = ( b - a ) / 2

( isto pode ser obtido, quer através da resolução do sistema do método dos coeficientes indeterminados, quer através dos integrais dos polinómios de Lagrange, l0(x) e l1(x) )

Temos assim a Regra dos Trapézios (simples):

T( f ) =  ( f(a) + f(b) ) ( b - a ) / 2

que corresponde exactamente ao valor da área do trapézio definido pela recta interpoladora!


Teorema (do Valor Intermédio para Integrais):
Sejam f , g funções contínuas em [a,b].
Se g não muda de sinal no intervalo [a, b] temos :
 

ó
õ
b

a

f(x) g(x) dx = f(x) ó
õ
b

a

g(x)dx,   para certo x Î[a, b]

Erro da Regra dos Trapézios (simples)

Pretendemos agora determinar majorações para o valor absoluto do erro
E( f ) = I ( f ) - T ( f )
Sabemos que
E( f ) = I ( f ) - T ( f ) = I ( f ) - I ( p1 ) = I ( f - p1 )
Da fórmula do erro de interpolação temos
f (x) - p1(x) = f [ a, b, x ] ( x - a ) ( x - b )

e como ( x - a ) ( x - b ) não muda de sinal no intervalo [a, b] podemos aplicar o Teorema do Valor Intermédio para Integrais e obtemos
 

ó
õ
b

a

f[a, b, x] (x-a)(x-b) dx = f[a, b, x] ó
õ
b

a

(x-a) (x-b) dx,  para certo x Î[a, b]

e supondo que f é C2[a, b], obtemos a fórmula do erro:

E( f ) = - (b-a)3
  12
f '' (x),  para certo x Î[a, b]

Regra dos Trapézios (Composta)

Como é claro, se usassemos a regra dos trapézios simples para calcular o integral num intervalo [a, b], iamos obter uma aproximação grosseira do valor do integral, por isso, interessa-nos decompor esse intervalo em sub-intervalos cada vez mais pequenos, e nesses sub-intervalos aplicamos a regra dos trapézios simples.
Trata-se, neste caso, de fazer uma aproximação, da função integranda, usando um spline linear.

Para simplificar, consideramos que o tamanho desses sub-intervalos é constante = h.
Assim, definimos h = ( b - a ) / N, onde N é o número de sub-intervalos ( = número de nós - 1), e temos:

xi = a + i h
portanto, o valor do integral é igual à soma dos integrais nos sub-intervalos.
 
I(f ) =  ó
õ
b

a

f(x) dx =   N
S
i =1
ó
õ
xi

xi-1

f(x)dx, 

logo, aplicando a regra dos trapézios simples a cada um desses sub-intervalos, obtemos
 

TN(f ) =  h
 N
S
i =1
(f(xi-1)+f(xi))

Reparando que há termos que aparecem repetidos na soma, podemos simplificar a expressão, e obtemos a Regra dos Trapézios Composta:
 

TN(f ) = h ( f(a) + f(b)
       2 
 +   N-1
S
i =1
f(xi))

Erro da regra dos trapézios composta

Como já obtivemos a fórmula do erro para a regra simples, o erro da regra dos trapézios composta será a soma dos erros cometido em cada sub-intervalo, ou seja :
 
EN(f ) =  N
S
i =1
h3
12 
f '' (xi) = h3N
12 
( 1
N
 N
S
i =1
f '' (xi) )
.
como N h = b - a , e como podemos aplicar o teorema clássico do valor intermédio à média das 2as derivadas, obtemos a fórmula do erro da regra dos trapézios composta:
EN(f ) = (b-a) h2
   12 
f '' (x), onde [a,b] 

Regra de Simpson

Regra de Simpson (simples)

Tal como a Regra dos Trapézios, trata-se de outro exemplo de Fórmula de Newton-Cotes fechada, mas, ao invés de considerarmos a aproximação em cada sub-intervalo através de um polinómio interpolador do 1º grau (recta), podemos pensar numa aproximação um pouco melhor, considerando um polinómio interpolador do 2º grau (parábola). Para isso, ao considerarmos a regra de integração simples, precisamos de um ponto adicional, que será o ponto médio.
.
Neste caso,
n = 2 h = (b-a)/2 x0 = a x1= c = a+h x2 = b

a fórmula de integração será do tipo

Q( f ) = A0 f(a) + A1f(c) + A2 f(b)

e podemos obter os pesos A0, A1, A2, resolvendo o sistema linear

é
ê
ê
ë
1  1 1
a  c  b
a2  c  b2
ù
ú
ú
û
é
ê
ê
ë
A0
A1
A2
ù
ú
ú
û
 =  é
ê
ê
ë
b - a
(b2 - a2)/2
(b3 - a3)/3
ù
ú
ú
û
.
ou através dos polinómios de Lagrange:
 
A0 ó
õ
b

a

l0(x) dx =  ó
õ
b

a

(x-a)(x-c)
(a-b)(a-c)
dx = h / 3 
A1 ó
õ
b

a

l1(x) dx =  ó
õ
b

a

(x-b)(x-c)
(b-a)(b-c)
dx = 4h / 3 
A2 ó
õ
b

a

l2(x) dx =  ó
õ
b

a

(x-a)(x-b)
(c-a)(c-b)
dx = h / 3 

Obtemos, assim, a Regra de Simpson (simples):

S( f ) = ( f (a) + 4 f (c) + f (b) ) h / 3


Erro da Regra de Simpson (simples)

Neste caso, usamos um polinómio interpolador do 2º grau, e sabemos da fórmula do erro de interpolação que:
f (x) - p2 (x) = f [ a, b, c, x ] (x - a) (x - b) (x - c)

portanto

ES(f ) = I(f ) - S(f ) = I(f) - I(p2) ó
õ
b

a

f[a, b, c, x] (x-a) (x-b) (x-c) dx
No entanto, aqui não podemos aplicar directamente o teorema do valor intermédio para integrais, porque (x-a)(x-b)(x-c) muda de sinal no intervalo [a, b].

 Introduzindo um ponto auxiliar d pertencente ao intervalo [a, b], podemos usar a definição de f [ a, b, c, d, x ] para obter
 
 

f [ a, b, c, x ] = f [ a, b, c, d ] + f [ a, b, c, d, x ] ( x - d )

Substituimos esta expressão no integral, e como

ó
õ
b

a

(x-a) (x-b) (x-c) dx = 0
obtemos, fazendo d tender para c, e aplicando o teorema do valor intermédio para integrais, porque (x-a)(x-b)(x-c)2 já não muda de sinal no intervalo [a, b]:
admitindo que a função f está, pelo menos, em C4[ a, b ].
(Observamos que a repetição dos nós nas diferenças divididas leva a expressões mal definidas, que são indeterminações, pelo que esta notação só faz sentido como um limite, sabendo que a função f é diferenciável.
Por exemplo:
f [a, b, c, c, x] =    lim
c'®c
f [a, b, c, c', x]
e recursivamente, através da definição de diferenças divididas, podemos obter uma expressão em função de f [ c, c ] = f ' (c) ).


Portanto, se f pertence a C4[ a, b ], o Erro da Regra de Simpson (simples) fica:
 

ES(f ) = h5
90
 f (4)(x), onde [a,b]

Observação:
Por construção, sabiamos que a Regra de Simpson era exacta para polinómios do 2º grau, mas esta fórmula do erro revela-nos que ela é exacta mesmo para polinómios do 3º grau! Portanto, enquanto que a Regra dos Trapézios tem apenas grau 1, a Regra de Simpson tem grau 3.


Regra de Simpson (composta)


Regra de Simpson aplicada a dois sub-intervalos.

Convém referir que, enquanto a Regra dos Trapézios composta corresponde a fazer a aproximação da função integranda através de um spline linear, no caso da Regra de Simpson composta, a aproximação feita não corresponde a um spline de grau 2, pois não exigimos regularidade da derivada nos nós. Essa regularidade não é necessária quando integramos. Aliás, geometricamente depreende-se que, exigir a regularidade da função aproximadora, nos nós, não traz aparentes vantagens para a aproximação da área delimitada pelo gráfico da função...

Para aplicar a regra de Simpson usando sub-intervalos, devemos considerar um número impar de nós N+1, de forma que ao dividirmos o intervalo [ a, b ] em N/2 sub-intervalos, obtemos os nós

xi = a + i h
para i = 0, ..., N, com h = (b - a)/N.

Assim, podemos considerar três nós em cada sub-intervalo [ x2k-2, x2k ] :

x2k-2, x2k-1, x2k
para k = 0, ..., N/2, e aplicamos a regra de Simpson simples a cada um destes sub-intervalos.

Como
 

I(f ) =  ó
õ
b

a

f(x) dx =   N/2 
S
i =1
ó
õ
x2i

x2i-2

f(x)dx, 

obtemos

SN(f ) =  h
 N/2
S
i =1
(f(x2i-2) + 4 f(x2i-1) + f(x2i))

e simplificando os termos repetidos, temos a Regra de Simpson Composta:
 
 


Erro da Regra de Simpson Composta

Tal como no caso da Regra dos Trapézios composta, o erro da Regra de Simpson composta, resulta da soma dos erros em cada sub-intervalo, ou seja:
 
EN(f ) =  N/2 
S
i =1
h5
90 
f (4)(xi) = h5 N/2
 90 
( 2
N
 N/2
S
i =1
f (4)(xi) )
.a partir desta expressão, e de forma análoga, obtemos, a Fórmula do Erro para Regra de Simpson Composta:
 
EN(f ) = (b-a) h4
   180 
f (4)(x), onde [a,b] 


Fórmulas de Gauss

Até aqui vimos apenas as regras de Newton-Cotes, nas quais os nós de quadratura estão igualmente espaçados no intervalo de integração. Essa suposição condiciona o grau de quadratura da fórmula, já que sendo dados os nós z1, ... , zn na fórmula de quadratura

Q(f) = A1 f(z1 )+ ... + An f(zn )
os únicos graus de liberdade resultam de n incógnitas, que são os pesos A1, ... , An.
Para aumentar o grau de quadratura de uma fórmula podemos admitir que a própria localização dos nós é uma incógnita, duplicando assim o número de graus liberdade (já que z1, ... , zn passam a ser também incógnitas), mas a priori ficamos com um sistema não linear para resolver.
Normalmente consideramos a integração num intervalo de referência [-1,1], notando que podemos reduzir qualquer integração num intervalo [a, b] a este intervalo de referência com uma simples mudança de variável:
 
 I( f ) =  ó
õ
b

a

f(x) dx = ½ (b-a) ó
õ
1

-1

f (a + ½(t+1)(b-a) ) dt



Exemplos:


(I) Começamos por reparar que, no intervalo [-1, 1], se pretendermos utilizar apenas
um nó de integração, z1 ficamos com Q(f) = A1 f(z1 ), e tendo duas incógnitas A1 , z1
podemos exigir que se trate de uma fórmula de grau 1, resolvendo as duas equações:
I(1) = Q(1), I(x) = Q(x),
Retiramos da primeira equação 2 = b-a = I(1) = Q(1) = A e da segunda equação
0 = I(x) = Q(x) = A1 z1 = 2z1 , ou seja z1 = 0.
Portanto, no intervalo [-1, 1], obtemos uma fórmula de grau 1,
Q( f ) = 2 f(0),
que será válida num intervalo qualquer [a, b], efectuando a mudança de variável já citada, ficando com
Q( f ) = ( ½ (b-a) ) 2 f(a + ½ (b-a)) = (b - a) f ( ½(a + b) ),
que é designada por regra do ponto médio.
Tal como no caso da regra dos trapézios ou da regra de Simpson podemos definir uma regra do ponto médio composta, utilizando esta regra em cada subintervalo.

Observação: Poderá então colocar-se a questão de preferir a regra do ponto médio à regra dos trapézios, já que ambas as regras são exactas para polinómios de grau 1, e tendo em conta que a regra do ponto médio simples utiliza apenas um ponto, enquanto que a regra dos trapézios utiliza 2. No entanto, reparamos que, aplicadas a N subintervalos, a regra composta do ponto médio exige o cálculo da função em N pontos enquanto que a regra dos trapézios composta exige o cálculo em N+1 pontos, ou seja apenas mais um!



(II) Considerando agora 2 nós de integração, temos quatro variáveis, pois
Q(f) = A1 f(z1 ) + A2f(z2 ) ,
verifica as quatro equações
I(1) = Q(1), I(x) = Q(x), I(x2) = Q(x2), I(x3) = Q(x3),
desde que
ì
ï
ï
í
ï
ï
î
2  = A1 + A2
0 = A1 z1 + A2 z2
2/3 = A1 z12 + A2z22
0 = A1 z13 + A2z23
Este sistema pode ser resolvido facilmente começando por usar a segunda e quarta equação para obter
A1 z13 + A2 z2z12 = A1 z13 + A2 z23<=> A2 z2 (z12 - z22) = 0.
- caso A2 z2 = 0 teríamos uma contradição (pela 2ª e 3ª equação, 2/3 = A1 z12 = 0 z1= 0).
- caso z1 = z2 teríamos uma contradição (pela 2ª e 3ª equação, 2/3 = (A1 + A2) z12 = 0 z1= 0).
Resta a hipótese z1 = - z2 que implica (A1 - A2) z1= 0, logo A1 = A2 .
Assim, pela 1ª equação, A1 = A2 = 1, e pela 3ª equação, 2 z12 = 2/3 <=> z1= + 1 / Ö3 .
Concluímos que Q(f ) = f (-1/ Ö3 ) + f ( 1 / Ö3 ).
Por construção, trata-se de uma fórmula de quadratura que tem pelo menos grau 3,
concluindo-se que, como 2/5 = I(x4) ¹ Q(x4) = 2/9, tem exactamente grau 3 (tal como a regra de Simpson).

Observação: Tal como na comparação entre a regra do ponto médio e a regra dos trapézios concluímos que as regras compostas apenas diferiam no cálculo de um nó suplementar, para a regra dos trapézios composta, também se considerarmos uma regra composta baseada na fórmula de quadratura Q que acabamos de deduzir, verificamos que a regra de Simpson composta apenas exigiria o cálculo suplementar de um nó, e tem a vantagem de se considerar os nós igualmente espaçados.


Regras de Gauss-Legendre

Iremos agora verificar que é possível encontrar, de forma muito mais simples, os valores dos nós  z1, ..., zn e dos pesos A1, ..., An que permitem obter uma fórmula de quadratura de grau 2n -1, aproveitando os 2n graus de liberdade existentes.

Teorema:
Seja Pn o polinómio de Legendre de grau n, e z1, ..., zn as suas raízes.
A fórmula de quadratura Q(f) = A1 f(z1 )+ ... + An f(zn ), em que os pesos são

Ak ó
õ
1

-1

lk(x) dx ó
õ
1

-1

    n
 P
i=1, i ¹k
x - zi
zk
- zi
dx
tem grau 2n -1.
Demonstração:
Pretendemos ver que, sendo p2n-1 um polinómio de grau < 2n-1,
temos I(p2n-1) = Q(p2n-1).
Para isso consideramos a divisão de  p2n-1 por Pn, o que nos dá
p2n-1 (x) = qn-1(x) Pn(x)+rn-1(x),
em que qn-1 é o quóciente e rn-1 é o resto (ambos polinómios de grau < n -1).
Como Pn(zk) = 0,  obtemos imediatamente Q(p2n-1) = Q(rn-1).
Por outro lado, como rn-1 tem grau n-1, o polinómio interpolador em quaisquer
n pontos coincide com rn-1 (em particular, considerando os pontos z1, ..., zn).
Assim,  rn-1 pode ser escrito de acordo com a fórmula de Lagrange,
rn-1(x) = rn-1(z1) l1(x)+...+rn-1(zn) ln(x).
Portanto, I(rn-1) = rn-1(z1) I(l1) +...+rn-1(zn) I(ln) = Q(rn-1), já que Ak = I(lk).
Concluímos assim que Q(p2n-1) = Q(rn-1) = I(rn-1).
Resta ver que I(rn-1) = I(p2n-1) para concluirmos o resultado.
Ora, I(p2n-1) = I(qn-1 Pn)+ I(rn-1), e como I(qn-1 Pn) representa o produto
interno entre qn-1 e Pn , que são ortogonais, concluímos que I(qn-1 Pn) = 0.
Falta apenas verificar que não é exacta para um polinómio de grau 2n.
Para esse efeito, basta observar que Q(Pn2) = 0 ¹  I(Pn2) = ||Pn||2 .

Observação:
O mesmo tipo de raciocínio poderia ser facilmente aplicado a outros polinómios ortogonais, com outros produtos internos.
Por exemplo, definindo

J( f ) =  ó
õ
1

-1

w(x) f(x) dx
em que w(x) = 1 / Ö(1-x2) é o peso usado no produto interno que define os polinómios ortogonais de Chebyshev, então uma fórmula de quadratura que seria exacta para J( p3 ), em que  p3é um polinómio de grau < 3, seria
Q(f) = A1 f(z1 ) + A2 f(z2 ),
em que z1= -1/Ö2, z2= 1/Ö2, são as raízes do polinómio de Chebyshev T2(x) = 2x2 -1,
e em que
Ak ó
õ
1

-1

w(x) lk(x) dx ó
õ
1

-1

Ö(1-x2)     n
 P
i=1, i ¹k
x - zi
zk
- zi
dx
para k=1, 2 .