Métodos Computacionais Numéricos e Algébricos
com Maple
Integração Numérica
Fernando Deeke Sasse
CCT - UDESC
[email protected]
Integração Numérica
Quadratura de Gauss
Métodos de Newton-Cotes
Introdução
As fórmulas de Newton-Cotes formam os esquemas mais comuns de integração numérica.
Eles consistem na estratégia de substituir um função complicada ou que não possui primitiva
ou dada numa tabela, por uma função aproximada que é fácil de integrar, ou seja,
,
sendo
um polinômio da forma
.
O método mais simples consiste em escolher um polinômio constante, que corresponde à soma
de Riemann, um conceito usado para introduzir o conceito de integral definida.
Normalmente a aproximação é feita num pequeno intervalo. Uma integral definida num
intervalo grande é então realizada somando-se as integrais aproximadas em um dado número
de subintervalos. Tal soma de Riemann pode ser realizada de diferentes modos,
correspondentes à posição escolhida que define a altura de cada retângulo. Por exemplo,
consideremos a integral:
.
O valor exato desta integral é
> with(Student[Calculus1]):
> ApproximateInt(x^2, 0..1, 'partition' = 4, 'method' =
lower, 'partitiontype' = subintervals, 'output' = 'plot',
'boxoptions' = ['filled' = ['transparency'=.5]]);
1
0
1
x
1
f x dx, where f x = x2 and
A lower Riemann sum approximation of
0
the partition is uniform. The approximate value of the integral is
0.2187500000. Number of partitions used: 4.
A soma de Riemann superior é dada por
> value(%);
1
0
1
x
1
f x dx, where f x = x2 and
A lower Riemann sum approximation of
0
the partition is uniform. The approximate value of the integral is
0.2187500000. Number of partitions used: 4.
> ApproximateInt(x^2, 0..1, 'partition' = 4, 'method' =
upper, 'partitiontype' = subintervals, 'output' = 'plot',
'boxoptions' = ['filled' = ['transparency'=.5]]);
1
0
1
x
1
f x dx, where f x = x2 and
An upper Riemann sum approximation of
0
the partition is uniform. The approximate value of the integral is
0.4687500000. Number of partitions used: 4.
Uma melhor aproximação consiste na soma de Riemann de ponto médio:
> ApproximateInt(x^2, 0..1, 'partition' = 4, 'method' =
midpoint, 'partitiontype' = subintervals, 'output' =
'plot', 'boxoptions' = ['filled' = ['transparency'=.5]]);
1
0
1
x
1
f x dx, where f x = x2
A midpoint Riemann sum approximation of
0
and the partition is uniform. The approximate value of the integral is
0.3281250000. Number of partitions used: 4.
(1.2.1.1)
Método Trapezoidal
Um método que dá resultados ainda melhores que o da soma de Riemann é o chamado método
do trapézio, que consiste em aproximar a integral da função em cada intervalo por um
polinômio de primeiro grau. A integral da função num intervalo qualquer é dada então pela
soma das áreas de todos os trapézios. Por exemplo consideremos novamente a integral:
Tomando novamente 4 intervalos temos
> ApproximateInt(x^2, 0..1, 'partition' = 4, 'method' =
trapezoid, 'partitiontype' = subintervals, 'output' =
'plot', 'boxoptions' = ['filled' = ['transparency'=0.9]]);
1
0
1
x
1
f x dx using trapezoid rule, where f x = x2 and
An approximation of
0
the partition is uniform. The approximate value of the integral is
0.3437500000. Number of partitions used: 4.
O valor exato da integral é
.
Deduziremos agora uma fórmula para calcular a integral de uma função qualquer
, não
divergente em um dado intervalo
. Num dado subintervalo
a área do
correspondente trapézio é dada por
.
Suponhamos que temos
pontos igualmente espaçados
.
Consequentemente temos subintervalos de tamanho dado por
.
A integral de
no intervalo
é dada por
.
Usando a aproximação do trapézio em cada subintervalo temos
=
.
Exemplo 1 Calcule a integral
utilizando o método do trapézio com
Definamos o integrando:
> f:=x->x^2;
intervalos.
(1.2.2.1)
O espaçamento dos intervalos, para
> n:=10:a:=0:b:=1.:
> h:=(b-a)/n;
, é dado por
(1.2.2.2)
Definimos a malha de
pontos:
> X:=[seq(a+i*h,i=0..n)];
(1.2.2.3)
Implementemos a fórmula do método do trapézio:
> S:=h/2*(f(X[1])+f(X[n+1])+2*sum(f(X[k]),k=2..n));
(1.2.2.4)
Definamos um procedimento que tem como entrada uma função, um intervalo, e o número de
nós.
> restart:
> trapezio:=proc(f,a,b,n)
local h,k,X,S,g;
g:=unapply(f,x):
h:=(b-a)/n;
X:=[seq(a+i*h,i=0..n)];
S:=evalf(h/2*(g(X[1])+g(X[n+1])+2*sum(g(X[k]),k=2..n)))
;
end:
> G:=x^2;
(1.2.2.5)
Portanto, com
,
> trapezio(G,0,1,20);
0.3337500000
Com
0,
> trapezio(G,0,1,10000);
0.3333333350
Observamos que a convergência para o valor correto 1/3 é muito lenta.
(1.2.2.6)
(1.2.2.7)
Erro no método trapezoidal simples
Esta dedução segue este vídeo (http://www.youtube.com/watch?v=NAIGuh5I3Io) do Prof. L.
A. F. Coelho. Consideremos o erro cometido na aplicação do método trapezoidal aplicado a
um único intervalo. Se pudéssemos saber que a primitiva de
para a integral no intervalo
:
é
, poderíamos escrever
.
A integral aproximada é dada pela fórmula da área do trapézio:
.
O erro cometido na integração numérica é então definido por
.
Usando a expansão em série de Taylor temos
Como
,
Portanto, a integral exata torna-se
Por outro lado,
de modo que
e
.
Truncando a série em
por
temos que, de acordo com o teorema de Taylor, o erro é dado
,
onde
.
Erro no método trapezoidal composto
Somando os erros individuais para os diversos trapézios temos, usando o fato de que para
nós igualmente espaçados,
/ ,
,
sendo um número no intervalo
no intervalo
e escrever:
. Na prática podemos estimar o erro maximizando
Na fórmula acima o somatório começa em
, pois temos
é definido dentro do intervalo
.
Exemplo 1 Utilizando a regra do trapézio, calcule
pontos e intervalos. Cada
com erro de no máximo
.
Solução
Para calcular o erro máximo cometido devemos maximizar a segunda derivada de
, ou seja,
> restart:
> f:=exp(-x^3+x^2);
(1.2.2.8)
A segunda derivada desta função é dada por
> d2f:=diff(f,x$2);
(1.2.2.9)
Façamos o gráfico desta função:
> plot(abs(d2f),x=1..2);
3
2
1
1
2
x
Vemos imediatamente que o valor máximo procurado
,
ou seja,
> M:=abs(evalf(subs(x=1,d2f)));
(1.2.2.10)
Como
,
devemos agora calcular :
> eq:=10^(-5)=(2-1)^3*M/(12*n^2);
(1.2.2.11)
> sol:=solve(eq,n);
(1.2.2.12)
Portanto, se
usando a fórmula
o erro será garantidamente menor que
. Calculemos agora a integral,
.
Temos aqui :
> n:=159: a:=1:b:=2: h:=(b-a)/n:
> F:=x->exp(-x^3+x^2);
(1.2.2.13)
> X:=[evalf(seq(a+h*i,i=0..159))]:
Portanto,
> Iap:=h/2*(F(X[1])+F(X[160])+2*sum(F(X[i]),i=2..159));
(1.2.2.14)
O Maple fornece uma aproximação ótima para esta integral:
> Iap2:=evalf(int(F(x),x=1..2));
(1.2.2.15)
Comparando este resultado com o do método do trapézio temos:
> Erro:=Iap-Iap2;
(1.2.2.16)
Ou seja, aproximadamente
<
.
Exemplo 2 Utilizando a regra do trapézio, calcule
com erro de no máximo
/2.
> restart:
> f:=sin(x)/x;
(1.2.2.17)
A segunda derivada desta função é dada por
> d2f:=diff(f,x$2);
(1.2.2.18)
Façamos o gráfico desta função:
> plot(abs(d2f),x=0..1);
0
1
x
Vemos imediatamente que o valor máximo procurado
,
ou seja,
> M:=abs(evalf(limit(d2f,x=0)));
(1.2.2.19)
Como
,
devemos agora calcular :
> eq:=10^(-5)/2=(1-0)*M/(12*n^2);
(1.2.2.20)
> sol:=solve(eq,n);
(1.2.2.21)
Portanto, se
o erro será garantidamente menor que
usando a fórmula
/2 Calculemos agora a integral,
.
Temos aqui :
> n:=75: a:=0:b:=1: h:=(b-a)/n:
> F:=x->sin(x)/x;
(1.2.2.22)
> X:=[evalf(seq(a+h*i,i=0..75))]:
Portanto,
> Iap:=h/2*(limit(F(t),t=0)+F(X[76])+2*sum(F(X[i]),i=2..75))
;
(1.2.2.23)
(1.2.2.23)
O Maple fornece uma aproximação ótima para esta integral:
> Iap2:=evalf(int(F(x),x=0..1));
(1.2.2.24)
Comparando este resultado com o do método do trapézio temos:
> Erro:=Iap-Iap2;
(1.2.2.25)
Ou seja, aproximadamente
<
.
Exercício
1. Calcule a integral
utilizando o método trapezoidal composto com erro menor que
.
Métodos de Simpson-Kepler
O método do trapézio consistia em aproximar a área sob a curva, entre
por um
polinômio interpolante de primeiro grau (reta). A partir de agora vamos generalizar este
procedimento utilizando polinômios de graus superiores.
Método de Simpson 1/3 simples
Aqui aproximamos a função
na Fig. 1 Ou seja,
por um polinômio do segundo grau
,
com
=
.
Fig. 1
, como mostrado
Necessitamos de três pontos que serão interpolados por
ponto médio
. Além de e
escolheremos o
.
Os coeficientes deste polinômio são definidos pelas condições de interpolação:
,
,
Mais explicitamente,
> restart:
> p:=x->a0+a1*x+a2*x^2;
(1.2.3.1)
> EQ:=[p(a)=f(a),p((a+b)/2)=f(m),p(b)=f(b)];
(1.2.3.2)
> sol:=solve(EQ,[a0,a1,a2]);
(1.2.3.3)
O polinômio interpolante é então dado por
> assign(sol):
> p(x);
(1.2.3.4)
A integral correspondente é dada por
> int(p(x),x=a..b);
(1.2.3.5)
> simplify(%);
(1.2.3.6)
(1.2.3.6)
> factor(%);
(1.2.3.7)
Portanto,
=
Como
.
,
=
Método de Simpson 1/3 composto
O método desenvolvido acima pode agora ser aplicado de maneira composta. Na Fig. 2
ilustramos tal aplicação para
.
Fig. 2 (Extraído de Pauls Online Notes)
Temos então três parábolas que interpolam cada tríade
.
e a integral total é a soma das áreas sob cada parábola:
,
com
,
,
.
Portanto,
.
Notemos que os pontos com índices ímpares têm coeficiente 2 e os índices têm coeficiente 4.
De modo geral temos
.
Para implementar esta fórmula no Maple, devemos notar que os pontos da malha têm índices
que começam com índice 1 e não 0. Vamos então reescrever a fórmula acima como
.
Exemplo 3
Resolvamos novamente a integral do Exemplo 1
,
agora utilizando o método de Simpson 1/3. Tomemos 158 intervalos, para fazer uma
comparação deste resultado com o do Exemplo 1.
> restart:
> f:=x->exp(-x^3+x^2);
(1.2.3.8)
> a:=1:n:=158:h:=evalf(1/n):
> x:=[evalf(seq(a+h*i,i=0..n))]:
A integral aproximada é então dada por
> Iap1 :=(1/3)*h*(f(x[1])+4*(sum(f(x[2*j]),j=1..(1/2)*n))+2*
(sum(f(x[2*k+1]),k=1..(1/2)*n-1))+f(x[n+1]));
(1.2.3.9)
O valor dado pelo Maple é:
> Iap2:=evalf(int(f(t),t=1..2));
(1.2.3.10)
O erro é dado por
> Iap1-Iap2;
(1.2.3.11)
Ou seja, o erro diminuiu quadraticamente, relativamente àquele cometido pelo método
trapezoidal, com um número similar de intervalos.
De fato, o erro no método de Simpson 1/3 composto é dado por
.
Exemplo 4
Resolvamos a integral
,
pelo método de Simpson 1/3 composto, com erro de até
> restart:
> f:=exp(-sin(x));
.
(1.2.3.12)
> d4:=diff(f,x$4);
(1.2.3.13)
> plot(d4,x=0..1);
0
1
x
Certamente o maior valor em módulo desta função é 3, em
> m:=evalf(subs(x=1,d4));
.
(1.2.3.14)
Então
> eq:=10^(-5)=1/(180*n^4)*m;
(1.2.3.15)
e o número de intervalos deve ser
> nn:=fsolve(eq,n,0..10);
Portanto, devemos ter ao menos
Definamos
> F:=unapply(f,x);
(1.2.3.16)
( deve ser sempre um número par neste método).
(1.2.3.17)
Com
> a:=0:n:=8:h:=evalf(1/n);
(1.2.3.18)
(1.2.3.18)
> x:=[evalf(seq(a+h*i,i=0..n))];
(1.2.3.19)
a integral aproximada é então dada por
> Iap1 :=(1/3)*h*(F(x[1])+4*(sum(F(x[2*j]),j=1..(1/2)*n))+2*
(sum(F(x[2*k+1]),k=1..(1/2)*n-1))+F(x[n+1]));
(1.2.3.20)
A integral dada pelo Maple é
> Iap2:=evalf(int(F(t),t=0..1));
(1.2.3.21)
O erro é
> Iap2-Iap1;
(1.2.3.22)
que é menor do que o limite superior estabelecido (como esperado).
Método de Simpson-Kepler 3/8 simples
Aqui aproximamos a função
por um polinômio do terceiro grau
pontos entre
e
, como mostrado na Fig. 3.
, que intepola 4
(x1 ,y1 )
(x 2 ,y 2 )
(x ,y )
y = p(x)
x0
(x 3 ,y 3 )
y = f(x)
x1
x2
x3
x
Fig. 3 Regra de Simpson 3/8
Ou seja,
,
com
=
.
Suponhamos que os quatro pontos que serão interpolados por
são igualmente espaçados:
,
Ou seja,
Fazendo
e
interpolação:
,
,
.
, os coeficientes deste polinômio são definidos pelas condições de
,
,
Mais explicitamente,
> restart:
> X:=[a,a+h,a+2*h,a+3*h];
(1.2.3.23)
> Y:=map(f,X);
(1.2.3.24)
Geramos agora o polinômio que interpola estes pontos:
> p:=interp(X,Y,x);
(1.2.3.25)
Integrando no intervalo
temos
> I1:=int(p,x=a..a+3*h);
(1.2.3.26)
Este resultado pode ser consideravelmente simplificado:
> I2:=simplify(I1);
(1.2.3.27)
Portanto,
=
,
com
.
Método de Simpson-Kepler 3/8 composto
Fazendo a composição do método de Simpson-Kepler 3/8 no intervalo
intervalos de comprimento
, com
,
temos
.
O erro tem a seguinte expressão:
,
com
. Este erro é da mesma ordem daquele cometido no método de Simpson 1/3.
Notemos que o número de intervalos agora pode ser ímpar, mas deve ser divisível por 3.
Exemplo 5
Resolvamos a integral
,
pelo método de Simpson 3/8 composto, com erro de até
> restart:
> f:=exp(-sin(x));
.
(1.2.3.28)
> d4:=diff(f,x$4);
(1.2.3.29)
> plot(d4,x=0..1);
0
1
x
Certamente o maior valor em módulo desta função é 3, em
> eq:=10^(-5)=1/(80*n^4)*3;
. Então
(1.2.3.30)
e o número de intervalos deve ser
> nn:=fsolve(eq,n,0..10);
Portanto, devemos ter ao menos
método). Definamos
> F:=unapply(f,x);
(1.2.3.31)
( deve ser sempre um número divisível por 3 neste
(1.2.3.32)
Com
> a:=0:n:=9:m:=n/3;h:=evalf(1/n);
(1.2.3.33)
> x:=[evalf(seq(a+h*i,i=0..n))];
(1.2.3.34)
a integral aproximada é então dada por (lembrando que devemos deslocar o índice de uma
unidade para fazer a implementação):
> Iap1 :=(3/8)*h*(sum(F(x[3*k-2])+3*F(x[3*k-1])+3*F(x[3*k])+
F(x[3*k+1]),k=1..m));
(1.2.3.35)
A integral dada pelo Maple é
> Iap2:=evalf(int(F(t),t=0..1));
(1.2.3.36)
O erro é
> Iap2-Iap1;
(1.2.3.37)
que é menor do que o limite superior estabelecido (como esperado).
Exercícios
1. Use a regra de Simpson 1/3, com
de
a
, para calcular a integral de
. Calcule o erro verdadeiro e o erro estimado.
2. Resolva a integral
(a) Analiticamente
(b) Usando uma única aplicação da regra do trapézio,
(c) Usando múltiplas aplicações da regra do trapézio, com
e
.
(d) Uma única aplicação da regra de Simpson 1/3.
(e) Múltiplas aplicações da regra de Simpson 1/3, com
.
Em cada item estime o erro e compare com o erro exato, comparando com (a).
3. Calcule a integral dos dados tabelados abaixo, usando (a) a regra trapezoidal, (b) a regra de
Simpson 1/3 e (c) o polinômio interpolador dos dados.
4. Calcule as integrais do exercício 1, com
, usando Simpson 1/3 e estime o
erro cometido. Compare com os valores fornecidos pelo sistema.
5. Determine
(i)
,
(ii)
(iii)
com um erro de no máximo 0.001,utilizando:
(a) O método dos trapézios
(b) O método de Simpson 1/3
(c) Utilize o comando de integração numérica do sistema e verifique se o erro real é realmente
menor ou igual ao previsto.