Aula Expositiva 13
4.4 Algoritmos Numéricos
4.4.1 Integração por Trapézios
4.4.2 Bisseção
4.4.3 Série de Taylor para exp(x) e
Cancelamento Catastrófico
DCC 001
Programação de Computadores
2º Semestre de 2011
Prof. Osvaldo Carvalho
DCC001 - 2011 - 2
1
Integração por Trapézios
DCC001 - 2011 - 2
2
Integração Numérica
 f x dx
 A integral indefinida de f  x 
 Queremos calcular
b
a
é dificilmente
obtida ou não existe
 O primeiro passo consiste em dividir o
intervalo a, b em n sub-intervalos iguais
a
x1
DCC001 - 2011 - 2
x
x2
x3
b
xn
xn+1
3
Aproximações por Retângulos
 Após a divisão em sub-intervalos, a área
embaixo da curva pode ser aproximada
por uma soma de retângulos
 Temos duas formas de definição de
somas:
Soma de Riemann pela esquerda
 Soma de Riemann pela direita

DCC001 - 2011 - 2
4
Área Coberta pela Soma de Riemann
pela Esquerda
O lado esquerdo de
cada retângulo
coincide com a curva
Área =
f(xi).∆x
DCC001 - 2011 - 2
xi
xi+1
5
Área Coberta pela Soma de Riemann
pela Direita
O lado direito de
cada retângulo
coincide com a curva
Área =
f(xi+1).∆x
DCC001 - 2011 - 2
xi
xi+1
6
Somas de Riemman
O lado esquerdo de
cada retângulo
coincide com a curva
DCC001 - 2011 - 2
O lado direito de
cada retângulo
coincide com a curva
7
Fórmulas para Somas de Riemann
Aesq

f  x1  x  f  x2  x    f  xn  x
n
  x  f  xi 
i 1
Adir

f x2  x  f x3  x    f xn1  x
n 1
  x  f xi 
i 2
DCC001 - 2011 - 2
8
Função LeftRiemannSum
function lrs = LeftRiemannSum(f,a,b,n)
// Soma de Riemann esquerda da função
// f entre os pontos a e b com n intervalos
x = linspace(a,b,n+1);
delta_x = (b-a)/n;
lrs = sum(f(x(1:n))) * delta_x;
endfunction
DCC001 - 2011 - 2
9
Função RightRiemannSum
function rrs = RightRiemannSum(f,a,b,n)
// Soma de Riemann direita da função
// f entre os pontos a e b com n intervalos
x = linspace(a,b,n+1);
delta_x = (b-a)/n;
rrs = sum(f(x(2:n+1)))* delta_x;
endfunction
DCC001 - 2011 - 2
10
Somas de Riemann com 16 subIntervalos
DCC001 - 2011 - 2
11
Área Coberta pela Soma de
Trapézios
Área do sub-intervalo =
∆x.(f(xi) + f(xi+1))/2
DCC001 - 2011 - 2
xi
xi+1
12
Fórmula para Aproximação por
Trapézios
A 

DCC001 - 2011 - 2
f  x1  f  x2 
2
x 
f  x2  f  x3 
2
x 
f  xn1  f  xn 
2
n
x 
 f  x1  f  xn1 

x 
  f xi 
2
i 2


f  xn  f  xn1 
2
x
13
Função TrapezoidalSum
function area = TrapezoidalSum(f,a,b,n)
// Calcula a área sob a curva f entre a e b,
// utilizando n pontos e a fórmula dos
// trapézios
x = linspace(a,b,n+1);
delta_x = (b-a)/n;
area = ( (f(x(1))+f(x(n+1)))/2 + ...
sum(f(x(2:n))) ...
)*delta_x;
endfunction
DCC001 - 2011 - 2
14
Comparação de Convergência
n
2
4
8
16
ERRO ESQUERDA
4.446396327e-001
2.092337399e-001
1.013895985e-001
4.989070473e-002
DCC001 - 2011 - 2
ERRO DIREITA
3.407585307e-001
1.834653418e-001
9.495994231e-002
4.828406570e-002
ERRO TRAPEZIO
5.194055103e-002
1.288419903e-002
3.214828114e-003
8.033195149e-004
15
Bisseção
DCC001 - 2011 - 2
16
Raízes de uma Função
 Problema: dada uma função f, contínua, encontrar
x tal que f(x) = 0
 Para algumas funções, como um polinômio de 2º
grau, este problema pode ser resolvido
analiticamente (por uma fórmula)
 Para outras não existe solução analítica, e
devemos usar métodos numéricos
 Ponto de partida: dois pontos a e b, tais que f(a)
< 0 e f(b) > 0, ou f(a) > 0 e f(b) < 0
DCC001 - 2011 - 2
17
F(a) < 0, F(b) > 0
1.5
1.0
0.5
a
0.0
3.0
3.5
4.0
4.5
5.0
5.5
b
6.0
6.5
-0.5
DCC001 - 2011 - 2
18
F(a) > 0, F(b)< 0
1.2
1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.5
a
1.0
b
1.5
2.0
2.5
3.0
3.5
-0.2
-0.4
-0.6
-0.8
DCC001 - 2011 - 2
19
Mais de uma raiz entre a e b
0.8
0.6
0.4
0.2
0.0
0
a
b
5
10
15
-0.2
-0.4
DCC001 - 2011 - 2
20
F(a)>0, F(b)>0 – Caso 1
0.8
0.6
0.4
0.2
0.0
0
5
10
15
-0.2
-0.4
DCC001 - 2011 - 2
21
F(a)>0, F(b)>0 – Caso 2
0.8
0.6
0.4
0.2
0.0
0
5
10
15
-0.2
-0.4
DCC001 - 2011 - 2
22
F(x) = 1/x
F(a) < 0, F(b) > 0, mas….
f(x) = 1/x
6
4
2
0
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
-2
-4
-6
DCC001 - 2011 - 2
23
F(a).F((a+b)/2) < 0
Antes, sabemos que
uma raiz está aqui
Teste simples para
sinais opostos
F(a)
0
m = (a+b)/2
b
a
F(m)
F(b)
Depois, sabemos que uma raiz está
aqui; podemos fazer b = m
DCC001 - 2011 - 2
24
F(b).F((a+b)/2) < 0
Antes, sabemos que
uma raiz está aqui
F(m)
F(a)
0
a
m = (a+b)/2
Depois, sabemos que
uma raiz está aqui;
podemos fazer a = m
DCC001 - 2011 - 2
b
F(b)
25
Parada do Algoritmo
 A cada passo o tamanho do intervalo que contém
a (uma) raiz é dividido por 2
 Quando devemos parar?
 Cálculos numéricos sempre são aproximados;
testar se f(x) == 0 pode levar a loops infinitos
 Solução: o algoritmo deve parar quando o
intervalo [a,b] for suficientemente pequeno, isto
é, quando b-a for menor que uma tolerância
fornecida pelo usuário.
DCC001 - 2011 - 2
26
A Função Bissecao
Cabeçalho
function r = Bissecao(f, a, b, tol)
// se f é contínua e se f(a).f(b) < 0,
esta
// função calcula a raiz r com
// precisão menor ou igual
// ao valor de tol
Tolerância
Função da qual se
endfunction
deseja conhecer
uma raiz
DCC001 - 2011 - 2
Intervalo que
contém a raiz
27
Teste da Função Bissecao
 Para testar a função Bissecao, nós
precisamos
 De
uma função contínua
 De um intervalo onde a função troca de
sinal
 E de conhecer o valor de uma raiz nesse
intervalo para poder verificar o resultado
DCC001 - 2011 - 2
28
A Função e^(x)sin(x)
π
a
DCC001 - 2011 - 2
b
29
A função exp_sin(x)
function y = exp_sin(x)
y = exp(-x) .* sin (x);
endfunction
 Reparem no uso de “.*”, ao invés de “*”
 Isso é essencial para que a função aceite um
vetor como argumento de entrada e produza
um vetor como argumento de saída
DCC001 - 2011 - 2
30
Testando a Função Bissecao
clear
exec("exp_sin.sci");
exec("bissecao.sci");
tolerancia = input("\nTolerância = ");
while tolerancia > 0
raiz = bissecao(exp_sin,2,4, tolerancia);
printf(" Raiz = %12.10f;\n
Pi =
%12.10f\n",raiz,%pi);
tolerancia = input("\nTolerância = ");
end
DCC001 - 2011 - 2
31
A função Bissecao
Loop
function r = bissecao(f,a,b,tol)
while b-a > tol
// Redução do intervalo
end
r = (a+b)/2;
endfunction
DCC001 - 2011 - 2
32
A função Bissecao
Redução do Intervalo
// Redução do intervalo
m = (a+b)/2; //Ponto médio
if f(a)*f(m) <= 0 then
// [a,m] contém uma raiz
b = m;
end
if f(m)*f(b) <= 0 then
// [m,b] contém uma raiz
a = m;
end
DCC001 - 2011 - 2
33
Teste
Tolerância = 1.0e-3
Raiz = 3.1411132813;
Pi = 3.1415926536
Tolerância = 1.0e-6
Raiz = 3.1415925026;
Pi = 3.1415926536
Tolerância = 1.0e-10
Raiz = 3.1415926536;
Pi = 3.1415926536
DCC001 - 2011 - 2
34
Série de Taylor para exp(x) e
Cancelamento Catastrófico
DCC001 - 2011 - 2
35
Série de Taylor para exp(x)
 Do cálculo sabe-se que, para qualquer x,
exp(x) pode ser calculado pela soma de
infinitos termos
DCC001 - 2011 - 2
36
Série de Taylor para exp(x)
 A partir do ponto onde n >= x, n!
cresce mais rapidamente que xn
 Na prática a soma é feita até que o
valor absoluto do termo calculado
seja menor que uma tolerância
desejada
DCC001 - 2011 - 2
37
Função expTaylor
Primeira Versão
function y = expTaylor(x,tol)
// Calcula a soma dos termos
// da série de Taylor até o primeiro
// termo com valor absoluto menor
// que a tolerância tol
endfunction
DCC001 - 2011 - 2
38
Testando a função expTaylor
exec("expTaylor.sci");
tol = input("\ntol = ");
x = input("\nx = ");
while x ~= 999
expCalc = expTaylor(x,tol);
printf ("\n%15g %15.8e %15.8e %15.8e\n", ...
x,exp(x),expCalc,exp(x)-expCalc
)
x = input("\nx = ");
end
Função Scilab, muito
DCC001 - 2011 - 2
confiável
39
Cálculo dos Termos da Série de Taylor
t0
DCC001 - 2011 - 2
t1
t2
t3
40
A Função expTaylor
function y = expTaylor(x,tol)
Termo = 1;
y = Termo;
i = 1;
while abs(Termo) >= tol
Termo = Termo * x / i;
y = y + Termo;
i = i+1;
end
endfunction
DCC001 - 2011 - 2
41
Teste da função expTaylor
x positivo
 Bons resultados para x positivo
 Erro 16 ordens de grandeza menor que os
valores calculados
tol = 1.0e-40
x = 1
1 2.71828183e+000 2.71828183e+000 -4.44089210e-016
x = 10
10 2.20264658e+004 2.20264658e+004 7.27595761e-012
x = 30
30 1.06864746e+013 1.06864746e+013 -3.90625000e-003
DCC001 - 2011 - 2
42
Teste da função expTaylor
x negativo
 Péssimos resultados para x negativo
 x=-20: erro da mesma ordem de grandeza
dos valores
 x=-30: valor calculado negativo!
tol = 1.0e-40
x = -1
-1 3.67879441e-001 3.67879441e-001 -1.11022302e-016
x = -10
-10 4.53999298e-005 4.53999296e-005 1.39453573e-013
x = -20
-20 2.06115362e-009 5.62188447e-009 -3.56073085e-009
x = -30
-30 9.35762297e-014 -3.06681236e-005 3.06681237e-005
DCC001 - 2011 - 2
43
E agora?
 A fórmula para a série de Taylor é
provada matematicamente há
(literalmente) séculos
 A função exp_Taylor é uma
implantação simples e direta da
fórmula
 O que aconteceu?
DCC001 - 2011 - 2
44
Origem: Aritmética de Ponto Flutuante
 Uso de um número fixo de bits para
representação da mantissa
 Aritmética de números com grandes diferenças de
ordem de grandeza não funciona como esperado:
-->eps = 1.0e-23;
-->y = 1.0e23;
-->x = y + eps;
-->x == y
ans =
x
é
igual
a
y
bit
por
bit!
T
DCC001 - 2011 - 2
45
Soma de Números de Ordens de
Grandeza muito Diferentes
-->eps = 1.0e-23;
-->y = 1.0e23;
-->x = y + eps;
-->x == y
ans =
T
DCC001 - 2011 - 2
46
Dízimas Periódicas
 0,1 = 0.00011001100110011001100….
 Ou seja, 0,1 não tem representação exata em
binário
-->0.1+0.2 == 0.3
ans =
F
DCC001 - 2011 - 2
47
Cancelamento Catastrófico
 Ocorre em subtrações de números com valores
absolutos próximos, sempre que se usa um
número fixo de algarismos significativos, binários,
decimais ou em qualquer outra base
a
b
a+b
a-b
DCC001 - 2011 - 2
Precisão 6 Precisão 3
3,14159
3,14
3,13000
3,13
6,27159
6,27
0,01159
0,01
Diferença Dif/Precisão 6
0,00159
0,0506%
0,00000
0,0000%
0,00159
0,0254%
0,00159
13,7187%
48
Valores dos Termos da Série para x = 20
DCC001 - 2011 - 2
49
Valor Absoluto dos Termos da Série para x = -20
~4.3e007
15 ordens de
grandeza maior que o
resultado correto
DCC001 - 2011 - 2
50
Conclusões
 Cuidado ao somar números de tamanhos
muito diferentes e ao subtrair números de
tamanhos próximos
 A aritmética de ponto flutuante é melindrosa;
use funções de bibliotecas sempre que
possível
 Por outro lado, não se deixe levar pelo
pessimismo: programas numéricos funcionam
como esperado na maior parte dos casos
DCC001 - 2011 - 2
51