Coletânea de Exercícios de Cálculo Numérico com Soluções

Prévia do material em texto

DEPARTAMENTO
DE
ENGENHARIA MECÂNICA
MatemáticaComputacional
Colectânea de Exercícios
• Propostas de resolução
• Conceitos fundamentais
• Algoritmos básicos
• Gráficos ilustrativos
Leonel Fernandes
Miguel Matos Neves
Virgínia Infante
José Viriato
Licenciaturas em Engenharia Mecânica, Engenharia Aeroespacial e
Engenharia e Arquitectura Naval - Ano Lectivo 2007/08
Para uma melhor leitura imprimir este documento em frente e verso.
ÍNDICE
TERMINOLOGIA v
2 1
2.1 NÚMEROS INTEIROS:
Conversão da base b para a base 10 . . . . . . . . . . . . . . . . . . . . 1
2.2 NÚMEROS INTEIROS:
Conversão da base b1 para a base b2 . . . . . . . . . . . . . . . . . . . 2
2.3 NÚMEROS REAIS:
Conversão da base 10 para a base b < 10 . . . . . . . . . . . . . . . . . 3
2.4 NÚMEROS REAIS:
Conversão da base b1 para a base b2
(b1 6= 10, b2 6= 10 e ambas > 2) . . . . . . . . . . . . . . . . . . . . . . 4
2.5 NÚMEROS REAIS:
Conversão da base 10 para base b > 10 . . . . . . . . . . . . . . . . . . 6
2.6 SISTEMAS DE PONTO FLUTUANTE:
Erros de representação, Unidade de arredondamento, Overflow e Un-
derflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3 11
3.1 ARITMÉTICA EM SISTEMAS FP:
Cancelamento subtractivo . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2 FORMATO SIMPLES IEEE 754:
Erros de representação . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3 FORMATO SIMPLES IEEE 754:
Erros de representação e Operações 0/0 e 1/0 . . . . . . . . . . . . . . 17
4 19
4.1 CONDICIONAMENTO DE UMA FUNÇÃO . . . . . . . . . . . . . . 19
4.2 CONDICIONAMENTO DE UMA FUNÇÃO E
CANCELAMENTO SUBTRACTIVO . . . . . . . . . . . . . . . . . . 20
4.3 INTERPOLAÇÃO POLINOMIAL:
Formas de Lagrange e Newton . . . . . . . . . . . . . . . . . . . . . . 22
4.4 INTERPOLAÇÃO POLINOMIAL:
Algoritmo de Horner . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
5 27
5.1 INTERPOLAÇÃO POLINOMIAL:
Máximos e mínimos locais, Pontos de inflexão . . . . . . . . . . . . . . 27
ii ÍNDICE
5.2 INTERPOLAÇÃO POLINOMIAL:
Interpolação de Hermite, Nós duplos e diferença dividida confluente . 31
5.3 INTERPOLAÇÃO POLINOMIAL:
Interpolação de Hermite, Nós triplos e diferença dividida confluente . . 33
6 35
6.1 INTERPOLAÇÃO POLINOMIAL:
Spline quadrático . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
6.2 INTERPOLAÇÃO POLINOMIAL:
Spline cúbico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6.3 INTERPOLAÇÃO POLINOMIAL:
Erros em splines cúbicos . . . . . . . . . . . . . . . . . . . . . . . . . . 41
7 43
7.1 DIFERENCIAÇÃO NUMÉRICA:
Diferenças finitas de primeira e segunda ordens . . . . . . . . . . . . . 43
7.2 DIFERENCIAÇÃO NUMÉRICA:
Majorante do erro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
7.3 DIFERENCIAÇÃO NUMÉRICA:
Espaçamento óptimo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
7.4 DIFERENCIAÇÃO NUMÉRICA:
Ponto de inflexão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
7.5 DIFERENCIAÇÃO NUMÉRICA:
Espaçamento desigual . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
8 57
8.1 INTEGRAÇÃO NUMÉRICA:
Regras do trapézio simples e de Simpson simples e composta . . . . . 57
8.2 INTEGRAÇÃO NUMÉRICA:
Regra de Gauss-Legendre e de Gauss-Legendre-Lobatto . . . . . . . . 60
8.3 INTEGRAÇÃO NUMÉRICA:
Regras de Gauss-Legendre e do trapézio composta . . . . . . . . . . . 61
9 65
9.1 INTEGRAÇÃO NUMÉRICA:
Dedução e grau de regras de integração . . . . . . . . . . . . . . . . . 65
9.2 INTEGRAÇÃO NUMÉRICA:
Regra do trapézio corrigida . . . . . . . . . . . . . . . . . . . . . . . . 66
9.3 INTEGRAÇÃO NUMÉRICA:
Regra de Simpson adaptativa iterativa . . . . . . . . . . . . . . . . . . 68
10 75
10.1 DETERMINAÇÃO DE ZEROS DE EQUAÇÕES NÃO-LINEARES:
Métodos da Bissecção, Falsa Posição e Secante . . . . . . . . . . . . . 75
10.2 EQUAÇÕES NÃO-LINEARES:
Método da Bissecção e zeros múltiplos . . . . . . . . . . . . . . . . . . 79
ÍNDICE iii
10.3 EQUAÇÕES NÃO-LINEARES:
Método da Falsa Posição . . . . . . . . . . . . . . . . . . . . . . . . . . 80
11 83
11.1 EQUAÇÕES NÃO-LINEARES:
Natureza e localização de zeros de polinómios (metodologia) . . . . . . 83
11.2 EQUAÇÕES NÃO-LINEARES:
Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . 86
11.3 EQUAÇÕES NÃO-LINEARES:
Método do Ponto Fixo. Localização de zeros de polinómios e aceleração
de Aitken . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
12 95
12.1 SISTEMAS DE EQUAÇÕES LINEARES:
Método de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
12.2 SISTEMAS DE EQUAÇÕES LINEARES:
Método de Gauss com escolha de pivot total . . . . . . . . . . . . . . . 97
12.3 SISTEMAS DE EQUAÇÕES LINEARES:
Método de Gauss com escolha de pivot parcial . . . . . . . . . . . . . 101
13 107
13.1 SISTEMAS DE EQUAÇÕES LINEARES:
Métodos de Doolittle e Crout. Factorização LU . . . . . . . . . . . . . 107
13.2 SISTEMAS DE EQUAÇÕES LINEARES:
Factorização LU. Cálculo do determinante . . . . . . . . . . . . . . . 111
13.3 SISTEMAS DE EQUAÇÕES LINEARES:
Método de Choleski e factorização LDLT . . . . . . . . . . . . . . . . 113
13.4 SISTEMAS DE EQUAÇÕES LINEARES:
Método de Doolittle. Inversa, determinante e número de condição . . . 115
REFERÊNCIAS 119
TERMINOLOGIA
Listamos seguidamente as principais abreviaturas, símbolos e notação matemática
usados neste texto. No que respeita a matrizes e vectores, ambos são grafados a negro:
os primeiros em maiúsculas e os segundos em minúsculas. Também incluímos uma
tabela com as funções elementares, seus significados e respectivas expressões usadas
pela maior parte das linguagens de cálculo científico, nomeadamente em MATLAB.
Tabela de funções elementares
SÍMBOLO SIGNIFICADO EXPRESSÃO
expx ≡ ex Exponencial de x exp(x)
lnx ≡ loge x Logaritmo natural (base e) de x log(x) ou ln(x)
cosx Co-seno de x cos(x)
coshx ≡ (ex + e−x)/2 Co-seno hiperbólico de x cosh(x)
sinx ≡ sen x Seno de x sin(x)
sinhx ≡ (ex − e−x)/2 Seno hiperbólico de x sinh(x)
tanx ≡ tg x Tangente de x tan(x)
√
x ≡ x1/2 Raiz quadrada de x sqrt(x)
|x| Valor absoluto de x abs(x)
2
TÓPICOS: Números Inteiros: Conversão da base b para a base 10; Conversão
da base b1 para a base b2. Números Reais: Conversão da base 10 para a base b < 10;
Conversão da base b1 para a base b2 (b1 6= 10, b2 6= 10 e ambas > 2); Conversão da base
10 para base b > 10. Sistemas de Ponto Flutuante: Erros de representação,
Unidade de arredondamento, Overflow e Underflow.
LEITURAS RECOMENDADAS: Capítulo 1 (pp. 1–12) do livro [Pina(1995)]
2.1 NÚMEROS INTEIROS:
Conversão da base b para a base 10
Ache a representação decimal dos seguintes números inteiros:
a) (101101)2; b) (221)3; c) (427)8
CONCEITOS TEÓRICOS
Números inteiros positivos: Conversão da base b ≥ 2 para a base 10† (exp. (1.2.2),
p. 4, [Pina(1995)]):
(dndn−1 . . . d1d0)b = dnbn + dn−1bn−1 + · · ·+ d1b1 + d0b0
com 0 6 di < b− 1, i = 0, 1, . . . , n. Assim, existem (n+ 1) dígitos. ¤
RESOLUÇÃO
a)
(10110¸
d0
1)2 = 1× 25 + 0× 24 + 1× 23 + 1× 22 + 0× 21 + 1× 20
= 32 + 8 + 4 + 1 = (45)10 = 45
No MATLAB escrevemos bin2dec(’101101’), resultando ans = 45.
b)
(221)3 = 2× 32 + 2× 31 + 1× 30 = 25
†Para representar os inteiros de 0 a 10 necessitamos de quantos bits? Necessitamos de 4 bits ou
16 configurações, pois 23 < 10 < 24.
2 SEMANA 2
No MATLAB escrevemos base2dec(’221’,3), resultando ans = 25.
c)
(427)8 = 4× 82 + 2× 81 + 7× 80 = 4× 64 + 16 + 7 = 256 + 23 = 279
No MATLAB escrevemos base2dec(’427’,8), resultando ans = 279. ¥
2.2 NÚMEROS INTEIROS:
Conversão da base b1 para a base b2
Obtenha a representação octal (base 8) dos números inteiros:
a) (101001)2; b) (1000111)2
CONCEITOS TEÓRICOSNúmeros inteiros positivos: Conversão da base 10 para a base b > 2 (p. 5,
[Pina(1995)]):
dnb
n + dn−1bn−1 + · · ·+ d1b + d0
b
= dnbn−1 + dn−1bn−2 + · · ·+ d1 + d0
b
sendo d0 o resto da divisão. ¤
RESOLUÇÃO
a) Apresentamos duas resoluções alternativas:
Resolução 1:
Converter da base 2 para a base 10 e posteriormente para a base octal, ou seja,
(101001)2 → ( )10 → ( )8.
( )2 → ( )10:
Aplicamos a exp. (1.2.2), p. 4, [Pina(1995)], pelo que (101001)2 = (41)10 = 41.
( )10 → ( )8:
(101001)2 = 41 = (51)8, pois efectuamos as divisões (inteiras) e tomamos os restos
na ordem inversa:
41 |8
1 5 |8
5 0
Resolução 2:
Converter da base 2 para a base 2n, com n ≥ 2 e inteiro, através de agrupamentos
de n dígitos. Neste caso n = 3, pois na base 8 temos os dígitos 0, 1, 2, . . . , 6, 7, que
se representam na base 2 como (000)2, (001)2, (010)2, . . . , (110)2, (111)2, respectiva-
mente. Por outras palavras, necessitamos de 3 bits para representar os algarismos 0
a 7. Esquematicamente temos
(101001)2 → ( )8.
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 3
Assim, (101 | 001)2 → (22 + 20 | 20)8 → (5 | 1)8, ou seja,
(101001)2 = (51)8
b) À semelhança do que fizemos na alínea a), obtemos
Resolução 1:
(1000111)2 → (26 + 22 + 21 + 20)10 = 64 + 4 + 2 + 1 = 71→
71 |8
7 8 |8
0 1 |8
1 0
(1000111)2 = (107)8
Resolução 2:
(001 | 000 | 111)2 → (20 | 0 | 22 + 21 + 20)8 → (1 | 0 | 7)8,
ou seja (101001)2 = (107)8 ¥
2.3 NÚMEROS REAIS:
Conversão da base 10 para a base b < 10
Determine a representação binária do número real (0.5)10.
CONCEITOS TEÓRICOS
Números reais: Conversão da base 10 para a base b (p. 7, [Pina(1995)]):
Para converter um real da base 10 para a base b procede-se segundo 4 passos:
i) Separar a parte inteira da parte não-inteira;
ii) Converter a parte inteira;
iii) Converter parte não-inteira, multiplicando sucessivamente por b, usando a
componente não-inteira do resultado para continuar e guardando a parte inteira para
a representação;
iv) Reunir as duas partes. ¤
RESOLUÇÃO
Aplicando os conceitos teóricos acima, teremos
i) Parte inteira: 0; Parte não-inteira: .5
ii) Converter a parte inteira: (0)10 = (0)2
iii) Converter parte não-inteira:
0.5× 2 = 1.0 −→ 1
0.0× 2 = 0.
Portanto, obtemos (.1)2 como representação de (0.5)10 = 12
iv) Reunindo vem (0.5)10 = (0)2 + (.1)2 = (.1)2
¥
Facilmente de verifica que (0.1)2 = 0× 20 + 1× 2−1 = 12 = (0.5)10 ¡
4 SEMANA 2
2.4 NÚMEROS REAIS:
Conversão da base b1 para a base b2
(b1 6= 10, b2 6= 10 e ambas > 2)
Dado o número real 437.125 na base octal, determine a sua representação binária com
7 dígitos após a vírgula.
CONCEITOS TEÓRICOS
Ver pp. 7–8, [Pina(1995)] e resolução do Exercício 2.3.
Números reais: Conversão da base b para a base 10 (exp. (1.3.2), p. 7, [Pina(1995)]):
(dndn−1 . . . d1d0.d−1d−2 . . . d−k)b =dnbn + dn−1bn−1 + · · ·+ d1b1 + d0b0
+d−1b−1 + d−2b−2 + · · ·+ d−kb−k,
sendo (dndn−1 . . . d1d0)b e (d−1d−2 . . . d−k)b os dígitos das parte inteira e não-inteira
do número, respectivamente. ¤
RESOLUÇÃO
Este Exercício pode ser resolvido de duas formas, indicadas de seguida.
Resolução 1:
Esquematicamente, temos dois passos principais
(437.125)8 → ( )10 → ( )2
Conversão da base octal para a base 10:
Aplica-se a exp. (1.3.2), p. 7, [Pina(1995)], ou seja
4× 82 + 3× 81 + 7× 80 + 1× 8−1 + 2× 8−2 + 5× 8−3 = 287 + 0.166015625
= 287.166015625
Dado que devemos considerar somente 7 dígitos após a vírgula
(437.125)8 ≈ 287.1660156.
Conversão da base 10 para a base binária:
i) Parte inteira: 287; Parte não-inteira: .1660156
ii) Converter a parte inteira:
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 5
287 | 2
1 143 | 2
1 71 | 2
1 35 | 2
1 17 | 2
1 8 | 2
0 4 | 2
0 2 | 2
0 1 | 2
1 0
287 = (100011111)2
iii) Converter a parte não-inteira:
0.1660156× 2 = 0.3320312 −→ 0
0.3320312× 2 = 0.6640624 −→ 0
0.6640624× 2 = 1.3281248 −→ 1
0.3281248× 2 = 0.6562496 −→ 0
0.6562496× 2 = 1.3124992 −→ 1
0.3124992× 2 = 0.6249984 −→ 0
0.6249984× 2 = · · ·
(0.125)8 = (0.0010101 . . .)2
iv)
(437.125)8 = (100011111.0010101 . . .)2
Considerando 7 dígitos após a vírgula, o resultado é
(437.125)8 ≈ (100011111.0010101)2
Resolução 2:
Uma vez que necessitamos de 3 bits para representar os números 0 a 7 (ver Reso-
lução 2 do Exercício 2.2), cada dígito do número (437.125)8 corresponde a um grupo
de 3 dígitos (3 bits) na base 2. Assim,
(437.125)8 = (4 | 3 | 7 | . | 1 | 2 | 5)8
→ (100 | 011 | 111 | . | 001 | 010 | 101)2
Como são pedidos 7 dígitos após a vírgula, ficamos com
(437.125)8 ≈ (100011111.0010101)2 ¥
Verificação da parte não-inteira: (.0010101)2 = 1 × 2−3 + 1 × 2−5 + 1 × 2−7 =
0.1640625. Existe então um erro na representação da parte não-inteira, que deveria
ser 0.166015625. Isto deve-se a termos truncado após o dígito d−7, quando seriam
necessários 9 dígitos. ¡
6 SEMANA 2
2.5 NÚMEROS REAIS:
Conversão da base 10 para base b > 10
Dado o número real 539.125 na base decimal, achar a sua representação na base
hexadecimal.
CONCEITOS TEÓRICOS
Quando a base é tal que b > 10, há que recorrer a outros símbolos para represen-
tar os dígitos 10, 11, 12, . . . Uma possibilidade é recorrer às letras do alfabeto latino
A,B,C, . . .
No caso da base hexadecimal ou base 16, resulta, 0, 1, . . . , 9, A,B,C,D,E, F , cor-
respondendo as letras aos dígitos
A = 10, B = 11, C = 12, D = 13, E = 14, F = 15. ¤
Após esta ressalva, a conversão da base 10 para a base b > 10 processa-se segundo
os quatro passos referidos nos Conceitos Teóricos do Exercício 2.3.
RESOLUÇÃO
Esquematicamente
(539.125)10 → ( )16
i) Parte inteira: 539; Parte não-inteira: .125
ii) Converter a parte inteira:
539 |16
59 33 |16
11 1 2 |16
2 0
539 = (21B)16
iii) Converter a parte não-inteira:
0.125× 16 = 2.0 −→ 2
0.000× 16 = 0.0 −→ 0
(0.125)10 = (0.2)16
iv) Reunindo a parte inteira à parte não-inteira fica
(539.125)10 = (21B.2)16
obtendo-se uma representação exacta. ¥
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 7
2.6 SISTEMAS DE PONTO FLUTUANTE:
Erros de representação, Unidade de arredonda-
mento, Overflow e Underflow
Considere os sistemas de ponto flutuante FP (10, 5, 2, A) e FP (10, 5, 2, T ).
a) Represente o número −327.258 nesses sistemas.
b) Qual o maior número em módulo representável nesses sistemas?
c) Qual o menor número em módulo representável nesses sistemas?
c) Calcule os erros relativos de representação e compare-os com as respectivas unidades
de arredondamento.
CONCEITOS TEÓRICOS
Um sistema FP (b, p, q) é constituído por todos os números reais x da forma
([Pina(1995)], pp. 9-10)
x = ±mbt em que

b−1 6 m 6 1− b−p
|t| 6 bq − 1
e ainda x = 0.
Portanto,
x = ±(.d−1d−2d−3 . . . d−p)b±(tq−1...t1t0),
sendo p um número finito de dígitos para a mantissa e q um número finito de dígitos
para o expoente, denominando-se b por base.
Salvo indicação em contrário, considera-se que a mantissa é normalizada, i.e.,
d−1 6= 0, exceptuando a representação do zero.
Truncatura (T ): Desprezam-se os dígitos do número real x que não cabem na
mantissa, i.e., os dígitos para além dos p primeiros não são incluídos na representação.
([Pina(1995)], p. 10)
Arredondamento (A): O número real x é representado pelo número do sistema
que lhe está mais próximo. ([Pina(1995)], pp. 10-11)
Limite de overflow: Maior número em módulo representável nesses sistemas, ou
seja, em termos computacionais, maior número que posso guardar em memória. Este
número não depende de se fazer arredondamento ou truncatura.
Limite de underflow: Menor número em módulo representável nesses sistemas.
Este número depende de se considerar a mantissa normalizada ou não.
Unidade de arredondamento, u: Majorante do erro relativo na representação de
um número num dado sistema FP (b, p, q),tal que ([Pina(1995)], p. 12)
u =

1
2b
1−p em FP (b, p, q, A)
b1−p em FP (b, p, q, T )
¤
8 SEMANA 2
RESOLUÇÃO
a) Concretizando para o caso em análise ficamos com
FP (10, 5, 2) : fl(x) = ±mb±t = ±(.d−1d−2d−3d−4d−5)10±(t1t0)
com 0 6 di 6 9, para i = −1,−2,−3,−4,−5.
Normalizando o número dado obtemos
x = −(.327258)10+3
verificando-se que tem 6 dígitos na mantissa.
Nos sistemas de ponto flutuante referidos teremos
fl(−327.258) = −(.32726)10+03 em FP (10, 5, 2, A)
fl(−327.258) = −(.32725)10+03 em FP (10, 5, 2, T )
b) O número limite de overflow é
(0.99999)10+99 em FP (10, 5, 2, A)
(0.99999)10+99 em FP (10, 5, 2, T )
sendo igual para T ou A.
c) Para o formato normalizado, o número limite de underflow é
(0.10000)10−99 em FP (10, 5, 2, A)
(0.10000)10−99 em FP (10, 5, 2, T )
sendo igual para T ou A.
Se tivermos um formato não normalizado este número é (0.00001)10−99. Neste
caso falamos em underflow gradual.
d) Atendendo à definição de unidade de arredondamento teremos neste caso
u =

1
210
1−5 = (0.50000)10−4 em FP (10, 5, 2, A)
101−5 = (0.10000)10−3 em FP (10, 5, 2, T )
O erro relativo de representação é
e =
x˜− x
x
,
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 9
onde x˜ = fl(x) é a representação no sistema de ponto flutuante e x é o número real a
representar. Assim, os erros relativo neste caso são∣∣∣−(0.32726)103−(−327.258)−327.258 ∣∣∣ = (.61114)10−7 < (0.50000)10−4 em FP (10, 5, 2, A)∣∣∣−(0.32725)103−(−327.258)−327.258 ∣∣∣ = (.24446)10−6 < (0.10000)10−3 em FP (10, 5, 2, T )
¥
3
TÓPICOS: Aritmética em Sistemas FP: Cancelamento Subtractivo. For-
mato Simples IEEE 754: Erros de representação; Erros de representação e Ope-
rações 0/0 e 1/0.
LEITURAS RECOMENDADAS: Capítulo 1 (pp. 12–27) do livro [Pina(1995)]
3.1 ARITMÉTICA EM SISTEMAS FP:
Cancelamento subtractivo
Determine no sistema FP (10, 4, 2, T ) as raízes da equação
x2 + 0.7341x+ (0.600)10−4 = 0
considerando que não existem dígitos de guarda† no processamento das operações em
ponto flutuante.
a) Usando a fórmula resolvente. Indique os erros absolutos Ex1 e Ex2 .
b) Justifique a origem do erro relativo obtido na menor raiz (em módulo) e sugira
uma forma de melhoria numérica para a resolução deste problema.
CONCEITOS TEÓRICOS
As operações aritméticas no sistema de ponto flutuante FP desenvolvem-se de
acordo com os seguintes passos (DAONA):‡
1. Decomposição dos operandos nas respectivas mantissas e expoentes: (mbt);
2. Alinhamento das mantissas, no caso de soma ou subtracção. Por exemplo, para
t1 > t2
m1b
t1 +m2bt2 = (m1 +m2bt2−t1)bt1
(.1)10−1 + (.5)10−2 = [(.1) + (.5)10−2+1]10−1;
3. Operação com mantissas e expoentes;
4. Normalização da mantissa. Por exemplo, (1.1)10−1 = (0.11)100;
5. Arredondamento ou truncatura na mantissa.
†O processador pode ter mais dígitos que a memória, designando-se os dígitos adicionais por
dígitos de guarda.
‡Em geral, as operações em FP não respeitam as propriedades comutativa, distributiva ou asso-
ciativa.
12 SEMANA 3
Cancelamento subtractivo: Verifica-se quando se subtraem números muito próxi-
mos no sistema de ponto flutuante utilizado. ¤
RESOLUÇÃO
a) Sabemos que
ax2 + bx+ c = 0⇔ x1,2 = −b±
√
b2 − 4ac
2a
e para a = 1, é matematicamente equivalente utilizar
x2 + bx+ c = 0⇔ x1,2 = −b±
√
b2 − 4c
2
.
A sequência de realização deste cálculo em ponto flutuante é:
fl(b) = (.7341)100
fl(b2) = (.7341× .7341)× (100 × 100) = (.5389 *Truncar028)100 = (.5389)100
fl(c) = (.6000)10−4
fl(4) = (.4000)10+1
fl(2) = (.2000)10+1
fl(4c) = (.4000× .6000)× (10−4 × 10+1) = (.2400)10−3
fl(b2 − 4c) = (.5389)100 − (.2400)10−3
= (.5389− .0002 *10
0 em evidencia e truncar
400)100 = (.5387)100
fl
(√
b2 − 4c
)
= [(.5387)100]1/2 = (.7339 :
Truncar
618519)100 = (.7339)100
No cálculo de
√
b2 − 4c assumimos que √x = x1/2 é uma função implementada,
não sendo necessário decompô-la em funções elementares.
Para a primeira raiz:
fl
(
−b−
√
b2 − 4c
)
= −(.7341)100 − (.7339)100
= −(1.4680) º
Normalizar e truncar
|100 = −(.1468)101
fl
(
−b−√b2 − 4c
2
)
=
−(.1468)101
(.2000)101
= −(.7340)100 = fl(x1)
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 13
Para a segunda raiz:
fl
(
−b+
√
b2 − 4c
)
= −(.7341)100 + (.7339)100
= −(0.0002) º
Normalizar
|100 = −(.2000)10−3
fl
(
−b+√b2 − 4c
2
)
=
−(.2000)10−3
(.2000)101
= −(1.) º
Normalizar
|10−4 = −(.1000)10−3 = fl(x2)
Para calcular os erros cometidos em aritmética FP é necessário conhecermos os
valores exactos das raízes, com 5 dígitos na mantissa:
x1,2 =
−b±√b2 − 4c
2
=
−0.7341±√0.73412 − 4× 0.00006
2
x1 = −(0.73402)100 ∨ x2 = −(0.81742)10−4
pelo que os erros absolutos e relativos efectivamente cometidos são:
|Ex1 | =
∣∣−(0.7340)100 − [−(0.73402)100]∣∣ = (0.20000)10−4
|Ex2 | =
∣∣−(0.1000)10−3 − [−(0.81742)10−4]∣∣ = (0.18258)10−4,
|ex1 | =
∣∣∣∣Ex1x1
∣∣∣∣ = ∣∣∣∣ (0.2000)10−4−(0.73402)100
∣∣∣∣ = (0.27247)10−4
|ex2 | =
∣∣∣∣Ex2x2
∣∣∣∣ = ∣∣∣∣ (0.18258)10−4−(0.81742)10−4
∣∣∣∣ = (0.22336)100 ≡ 22.3%.
Apesar dos erros absolutos serem praticamente iguais, verificamos que o mesmo
não acontece com o erro relativo. A segunda raiz apresenta um erro relativo superior
em quatro ordens de grandeza ao correspondente erro da primeira raiz.
b) Onde ocorreu o erro? Não aplicamos nós o mesmo processo a x1 e a x2!? Então,
o que terá sucedido?
Sabemos que entre as operações +,−,×,÷, em FP existe um cálculo perigoso
que é a subtracção de valores muito próximos, já que os primeiros dígitos da
mantissa se anulam e após a normalização surgem zeros à direita que, possivelmente,
o não seriam caso se utilizasse um maior número de dígitos.
De facto, há uma subtracção de valores muito próximos quando calculamos x2 em
FP (10, 4, 2, T ),
fl
(
−b+
√
b2 − 4c
)
= −(.7341)100 + (.7339)100 = −(0.0002)100
= −(.2000)10−3
14 SEMANA 3
Mas em aritmética com um maior número de dígitos na mantissa teríamos
fl
(
−b+
√
b2 − 4c
)
= −.7341 +
√
.73412 − 4× 0.0006
= −(0.16348)10−3
Dividindo este número por 2 obteríamos um valor mais próximo do valor exacto.
Portanto, a origem do problema no cálculo de x2 foi o cancelamento subtrac-
tivo, que se verifica quando subtraímos números muito próximos em FP.
Como contornar o problema? Sabendo nós que o problema é a subtracção, o que
há a fazer é evitá-la. Existem duas alternativas:
Alternativa 1: Manipulando simbolicamente a equação de segundo grau genérica,
obtemos
ax2 + bx+ c = a(x− x1)(x− x2) = a(x2 − x1x− x2x+ x1x2)
= ax2 − a(x1 + x2)x+ ax1x2,
e concluímos que
c = ax1x2 ⇔ x2 = c
ax1
.
Substituindo pelos valores conhecidos, c, a e x1, no sistema FP dado, obtemos
fl
(
c
x1
)
=
(.6000)10−4
−(.7340)100 = −(.8174
>
Truncar
38)10−4 = −(.8174)10−4 = x2.
Alternativa 2: Neste caso manipulamos a fórmula resolvente, tal que
x2 =
−b+√b2 − 4c
2
=
−b+√b2 − 4c
2
× −b−
√
b2 − 4c
−b−√b2 − 4c =
b2 − (b2 − ¸
2
4c)
2(−b−√b2 − 4c)
=
2c
−b−√b2 − 4c =
c
x1
.
Notar que obtemos a mesma expressão da alternativa 1, pois a = 1. ¥
3.2 FORMATO SIMPLES IEEE 754:
Erros de representação
a) Represente o número 512.15 em formato simples IEEE 754 com truncatura.
b) Determine o erro de representação cometido na alínea anterior. Se preferir indique
um majorante adequado do erro cometido.
c) Será possível representar 10−32 e −10+42 no formato simples? Justifique.
CONCEITOS TEÓRICOS
Na norma IEEE 754 a mantissa é normalizada, ou seja, o primeiro bit é sempre 1
e diz-se implícito.
O expoente é enviesado, isto é, é dado por e− 127.
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 15
O maior expoente é 11111110 e o menor 00000001, estando os expoentes 11111111e 00000000 reservados para:
00000000 → Número desnormalizado, prevendo underflow gradual.
x = (−1)s(0.d−1 . . . d−23)2e−126 → Limite de underflow gradual :
2−23 × 2−126 = 2−149 = 1.4× 10−45
11111111 → Ocorrência de overflow. Neste caso podem surgir duas mensagens:
Se m = 0→ NaN
Se m 6= 0→ +INF
¤
RESOLUÇÃO
a) Esquematicamente
(512.15)10 → FP (2, 24, 8, T )
pois segundo a norma IEEE 754, no formato simples, b = 2, p = 24 e q = 8.
Representar a parte inteira da mantissa:
512 | 2
0 256 | 2
0 128 | 2
0 64 | 2
0 32 | 2
0 16 | 2
0 8 | 2
0 4 | 2
0 2 | 2
0 1 | 2
1 0
512 = (1000000000)2
Representar a parte fraccionária da mantissa:
0.15× 2 = 0.30 −→ 0
0.30× 2 = 0.60 −→ 0
0.60× 2 = 1.20 −→ 1|
0.20× 2 = 0.40 −→ 0|
0.40× 2 = 0.80 −→ 0|
0.80× 2 = 1.60 −→ 1|
0.60× 2 = 1.20 −→ 1
· · ·
16 SEMANA 3
0.15 = (0.001001 . . .)2
Representar toda a mantissa:
512.15 = (1 000 000 000.00 1001 1001 1001 . . .)2
A norma IEEE 754 requer um bit implícito de valor 1 e o expoente enviesado:
x = (−1)s2e−127(1.d−1d−2d−3 . . . d−23)2
Assim,
512.15 = (1.000 000 000 00 1001 1001 1001 :
Truncar
1001 . . .)2 × 29
= (1.000 000 000 00 1001 1001 1001)2 × 29
Representar o expoente (enviesado):
e− 127 = 9⇔ e = 136
136 | 2
0 68 | 2
0 34 | 2
0 17 | 2
1 8 | 2
0 4 | 2
0 2 | 2
0 1 | 2
1 0
136 = (10001000)2
Temos então a seguinte alocação dos 32 bits, no formato simples IEEE 754:
0 10001000 00000000000100110011001
Sinal Expoente Mantissa
1 bit 8 bits 23 bits
b) Erro de representação cometido:
O valor exacto é
(512.15) = (1.000 000 000 00 1001 1001 1001 |1001 1001 1001 1001 . . .)2 × 29
Na norma IEEE 754 teremos o valor
(512.15)IEEE 754 = (20 + 2−12 + 2−15 + 2−16 + 2−19 + 2−20 + 2−23)× 29
= 512.1499634
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 17
e o erro em valor absoluto é
|E| = |512.1499634− 512.15| = (3.7)10−5
ou, em alternativa, podemos realizar o cálculo estimando o valor truncado
|E| = (2−24 + 2−27 + 2−28 + 2−31 + 2−32 + . . .)× 29 = (3.7)10−5
O majorante do erro cometido em truncatura é†
|E|maj = bt−p = 210−24 = 2−14 = (6.10)10−5
c) Para responder à pergunta colocada temos que determinar os limites de overflow
e underflow:
Limite de overflow : (1.111 . . . 1)× 2(11111110)2−127 = (1.111 . . . 1)× 2(254−127)
= (2− 2−23)× 2127 = (3.4)10+38
Limite de underflow : (1.000 . . . 0)× 2(00000001)2−127 = (1.000 . . . 0)× 2(1−127)
= 20 × 2−126 = (1.2)10−38
Portanto,
10−32 tem representação aproximada no formato simples IEEE 754, dado que
1.2× 10−38 < 10−32 < (3.4)10+38.
10+42 não tem representação aproximada no formato simples IEEE 754, originando
overflow, pois
10+42 > (3.4)10+38.
¥
3.3 FORMATO SIMPLES IEEE 754:
Erros de representação e Operações 0/0 e 1/0
a) Estime o erro cometido na representação do número (0.1)10 no sistema de ponto
flutuante FP (2, 24, 8, T ).
b) Determine o número limite de overflow e a unidade de arredondamento do sistema
de ponto flutuante em formato simples na norma IEEE 754. Indique a mensagem
gerada em aritmética IEEE 754 pelas operações 0/0 e 1/0.
RESOLUÇÃO
†O expoente t refere-se à representação em FP. Não esquecer que o formato IEEE 754 considerado
tem um bit implícito, o que não acontece em FP. Por essa razão o expoente é 10 e não 9.
18 SEMANA 3
Representação de (0.1)10 em FP (2, 24, 8, T ):
0.1× 2 = 0.2 −→ 0
0.2× 2 = 0.4 −→ 0|
0.4× 2 = 0.8 −→ 0|
0.8× 2 = 1.6 −→ 1|
0.6× 2 = 1.2 −→ 1|
0.2× 2 = 0.4 −→ 0
· · ·
Detectamos a presença de um padrão repetitivo, ou seja,
0.1 = (0.0 0011 0011 0011 0011 0011 . . .)2
= (0.1100 1100 1100 1100 1100 1100)2 × 2−3
e pelo facto de truncar no dígito −24 obtemos um erro de representação de:
|E| = |fl(x)− x| ≈ bt−p = 2−3−24 = (7.4506)10−9
b)
Limite de overflow : (20 + 2−1 + . . .+ 2−23)× 2(27+26+...21+20−127)
= (2− u)× 2127 ≈ 2128 = (3.402823)10+38
Unidade de arredondamento para truncatura : u = 21−24 = (1.192092)10−7
A 0/0 corresponde a mensagem NaN (Not a Number) e a 1/0 corresponde a men-
sagem INF. ¥
A resposta na janela de comando do MATLAB para cada um dos casos é:
» 0/0
Warning: Divide by zero.
(Type "warning off MATLAB:divideByZero" to suppress this warning.)
ans =
NaN
» 1/0
Warning: Divide by zero.
(Type "warning off MATLAB:divideByZero" to suppress this warning.)
ans =
Inf ¡
4
TÓPICOS: Condicionamento de uma Função. Condicionamento de
uma Função e Cancelamento Subtractivo. Interpolação Polinomial:
Formas de Lagrange e Newton; Algoritmo de Horner.
LEITURAS RECOMENDADAS: Capítulo 1 (pp. 23–27) e Capítulo 2 (pp. 35–
52) do livro [Pina(1995)]
4.1 CONDICIONAMENTO DE UMA FUNÇÃO
Estime o número de condição da função f (x) = tg(x2) na origem. A função é bem
ou mal condicionada nesse ponto? Justifique.
CONCEITOS TEÓRICOS
Número de condição da função f (x): Para analisar em aritmética de Ponto Flu-
tuante a influência do erro de representação do argumento x no cálculo de uma dada
função f (x), a Análise de Erros Directa pode não ser a mais adequada. A alternativa
é uma Análise de Erros Indirecta, a qual fornece uma indicação do factor de ampli-
ficação dos erros relativos baseada no número de condição de uma função num dado
ponto.
Considerando que x sofre um erro de representação, o argumento utilizado no
cálculo será x˜ e o erro de representação (erro absoluto) será E = x˜− x. Do teorema
do valor médio, e dado que x e x˜ podem considerar-se suficientemente próximos e f
suficientemente regular, obtém-se a seguinte relação:
f (x˜)− f (x) = f ′ (ξ)× (x˜− x)⇒ f (x˜)− f (x)
f (x)
=
f ′ (ξ)× (x˜− x)
f (x)
× x
x
ξ ∈ {[xmin − xmax] : xmin = min (x˜, x) e xmax = max (x˜, x)}
∣∣∣∣f (x˜)− f (x)f (x)
∣∣∣∣ ≈ ∣∣∣∣f ′ (x)× xf (x)
∣∣∣∣× ∣∣∣∣ x˜− xx
∣∣∣∣
ou, por outras palavras,
ef ≈ condf (x)× ex,
20 SEMANA 4
considerando-se geralmente no lugar de ex o valor da unidade de arredondamento u.
Diz-se que uma função é bem condicionada se condf (x) for pequeno e diz-se mal
condicionada nos restantes casos.
Para exemplificar melhor a distinção, recorre-se à seguinte comparação baseada
num cálculo em FP (2, 24, 8, A) onde a unidade de arredondamento é:
u =
1
2
b1−p =
1
2
21−24 ≈ (0.6) 10−7.
Com este valor de erro relativo de representação, obtém-se no caso de
condf (x) = 104 um ef = 10−3 ou seja 0.1% (f(x) é bem condicionada)
e no caso de
condf (x) = 106 um ef = 10−1 ou seja 10% (f(x) é mal condicionada).
Portanto, uma função ser bem condicionada ou mal condicionada é uma caracte-
rística relativa. ¤
RESOLUÇÃO
condf (x) =
∣∣∣∣f ′ (x)× xf (x)
∣∣∣∣ =
∣∣∣∣∣2× x×
(
1 + tg2(x2)
)× x
tg(x2)
∣∣∣∣∣
Lembrando-nos que na vizinhança de x = 0 se tem tg(x2) ≈ (x)2 facilmente se
conclui que condf(x = 0) = 2. Uma via mais trabalhosa é fazer o que se segue!
Em x = 0, teremos
condf (x = 0) =
∣∣∣∣∣2× 02 ×
(
1 + tg2(02)
)
tg(02)
∣∣∣∣∣ = 00
A indeterminação pode ser levantada aplicando sucessivamente a Regra de L’Hôpital-
Cauchy:
condf (x = 0) = lim
x→0
∣∣∣∣2x2 + 2x2tg2(x2)tg(x2)
∣∣∣∣
= lim
x→0
∣∣∣∣∣4x + 4xtg2(x2) + 8x3tg(x2)
(
1 + tg2(x2)
)
2x
(
1 + tg2(x2)
) ∣∣∣∣∣
=
4
2
= 2
A função é bem condicionada em torno da origem, pois os erros são apenas ampli-
ados para o dobro.
Por exemplo, em precisão simples — onde o majorante do erro relativo de repre-
sentação de x é de aproximadamente u = (0.6)10−7 — o erro relativo no cálculo da
função em FP será cerca de (1.2)10−7, ou seja, aproximadamente 0.000012%. ¥
4.2 CONDICIONAMENTO DE UMA FUNÇÃO E
CANCELAMENTO SUBTRACTIVO
a) Mostre que f(x) = sin(x)− x é bem condicionada na vizinhança da origem.
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 21
b) Calcule o valor de f(10−2) em FP (10, 3, 2, T ).
c) Determine os erros absolutos e relativos cometidos.
d) Em face do resultadoda alínea c) proponha uma forma alternativa para o cálculo
de f(10−2).
RESOLUÇÃO
a) A aplicação directa da expressão origina uma indeterminação que se pode levan-
tar aplicando sucessivamente a regra de L’Hôpital-Cauchy, em três aplicações. Uma
alternativa é utilizar a expansão em série
cond (f) =
∣∣∣∣f ′ (x)xf (x)
∣∣∣∣ = ∣∣∣∣ (cos (x)− 1)xsin (x)− x
∣∣∣∣ =
∣∣∣∣∣∣
([
1− x22! +O
(
x4
)]− 1)x[
x− x33! +O (x5)
]− x
∣∣∣∣∣∣
=
∣∣∣∣∣ x
3
2! +O
(
x5
)
x3
3! +O (x5)
∣∣∣∣∣ =
∣∣∣∣∣ 12 +O
(
x2
)
1
6 +O (x2)
∣∣∣∣∣ ≈ 121
6
= 3
A função é bem condicionada em torno da origem pois os erros são apenas ampli-
ados para o triplo. Por exemplo, em precisão simples — onde o majorante do erro
relativo de representação de x é de u = (0.6)10−7 — o erro relativo no cálculo da
função em FP será cerca de (1.8)10−7, ou seja, aproximadamente 0.000018%.
b) Para calcular f(x) = sin(x)− x no ponto x = 10−2 em FP (10, 3, 2, T ) tem-se:
fl(x) = fl(10−2) = +(0.100)10−01
Assumindo que sin(x) esteja implementada computacionalmente como função:
fl(sin(x)) = sin(0.100× 10−01) = 0.999 :Truncar9833× 10−02 = +(0.999)10−02
fl(sin(x)− x) = 0.999× 10−02 − 0.100× 10−01 = (0.099− 0.100)º
Alinhar× 10−01
= −0.001× 10−01
= −(0.100)10−03
Como iremos ver de seguida, a truncatura para um número reduzido de dígitos,
seguidamente agravada por um alinhamento das mantissas, é potencialmente perigosa
para a precisão dos resultados, ocorrendo cancelamento subtractivo.
c) O valor exacto, na precisão da calculadora,† é
f(10−2) = sin(10−02)− 10−02 = −(0.16666)10−6
pelo que
E = −0.100× 10−03 − (−0.166× 10−06) = −0.000099834
†Não esquecer de efectuar os cálculos em radianos e não em graus ou grados!
22 SEMANA 4
|e| =
∣∣∣∣ Ef(10−2)
∣∣∣∣ = 0.000099834/(0.16666× 10−06) = 599.03 ≡ 59903%
Obtemos, assim, um valor muito elevado.
Observa-se que f ser bem condicionada na vizinhança deste ponto não significa
imune ao cancelamento subtractivo, realçando-se neste exemplo que se tratam de dois
conceitos diferentes.
d) Para evitar o cancelamento subtractivo quando x é pequeno, pode utilizar-
se neste caso a expansão em série da função seno:
sin (x)− x =
[
x− x
3
3!
+
x5
5!
−O (x7)]− x ≈ −x3
6
fl(sin(x) − x) = fl(−x3/6) = −(0.100 × 10−01)/(0.600 × 10+01) = 0.166 × 10−06
obtendo-se a máxima precisão em FP (10, 3, 2, T ). ¥
No caso de f(x) = ln(1+x) com x < |u|, para evitar erros de FP usa-se o primeiro
termo da expansão ln(1 + x) = x− x2/2 + x3/3− x4/4 + . . .
Em MATLAB temos os logaritmos natural, log, de base 10, log10 e de base 2, log2,
enquanto u é obtido por eps. ¡
4.3 INTERPOLAÇÃO POLINOMIAL:
Formas de Lagrange e Newton
Determine o polinómio interpolador dos valores (0, 1); (1, 3) e (2, 2).
a) Na forma de Lagrange.
b) Na forma de Newton com centros nos nós.
RESOLUÇÃO
a) Na forma de Lagrange este polinómio obtém-se como se segue:
n é o número de pontos, ou nós, menos 1, ou seja o grau do polinómio (desde que
os pontos (xk, yk) dados não sejam colineares).
pn (x) =
n∑
k=0
Lk (x) yk com Lk =
n∏
i=0 i 6=k
x− xi
xk − xi
L0 =
x− x1
x0 − x1
x− x2
x0 − x2
=
(
x− 1
−1
)(
x− 2
−2
)
L1 =
x− x0
x1 − x0
x− x2
x1 − x2
= −x (x− 2)
L2 =
x− x0
x2 − x0
x− x1
x2 − x1
=
x
2
(x− 1)
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 23
Obtemos:
p2(x) = L0y0 + L1y1 + L2y2
p2(x) =
(x− 1) (x− 2)
2
− 3x (x− 2) + x (x− 1)
Podemos verificar rapidamente que p2(0) = 1 sendo L1 = L2 = 0; p2(1) = 3 com
L0 = L2 = 0 e p2(2) = 2 com L0 = L1 = 0. Tal é ilustrado pela representação gráfica
dos polinómios na Figura 4.3.1.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
−0.5
0
0.5
1
1.5
2
2.5
3
3.5
p(x)
L0(x)
L1(x)
L2(x)
Figura 4.3.1: Polinómio interpolador e polinómios de Lagrange
b) Na forma de Newton com centros nos nós o polinómio é obtido por:
pn (x) = y [x0] + y [x0, x1]W0 + y [x0, x1, x2]W1 + ...+ y [x0, x1, ..., xn]Wn−1
com
y [xk] = yk e y [xi, xi+1, ..., xk−1, xk] =
y [xi+1, ..., xk]− y [xi, ..., xk−1]
xk − xi
onde as diferenças divididas se obtêm como indicado na tabela seguinte:
x y[·] y[·, ·] y[·, ·, ·]
0 1
2
1 3 −3/2
−1
2 2
e Wi = (x− x0) . . . (x− xi−1) são os polinómios nodais.
24 SEMANA 4
O polinómio que se obtém na forma de Newton com centros nos nós é
p2(x) = 1 + 2x− 32x(x− 1)
Podemos verificar rapidamente que p2(0) = 1 sendo W0 =W1 = 0, p2(1) = 3 com
W1 = 0 e p2(2) = 1 + 4− 3 = 2. Na Figura 4.3.2 observa-se a representação de p(x)
e dos polinómios nodais Wi.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
−3
−2
−1
0
1
2
3
4
p(x)
2x
1
−
3
2
x(x− 1)
Figura 4.3.2: Polinómio interpolador e respectivas parcelas ¥
As diferentes representações polinomiais obtidas nas alíneas a) e b) são exacta-
mente o mesmo polinómio interpolador, pois este é único (ver Teorema 2.2.2. da
Unicidade, na p. 44 em [Pina(1995)]). ¡
4.4 INTERPOLAÇÃO POLINOMIAL:
Algoritmo de Horner
Determine o polinómio interpolador dos valores (0, -1); (1, 1); (2, 4) e (3, 2).
a) Na forma de Newton com centros nos nós.
b) Na forma de Lagrange.
c) Calcule p(4) pelo algoritmo de Horner.
ALGORITMOS
Cálculo do valor do polinómio num dado ponto x através de 3 algoritmos distintos
Algoritmo 1 - Cálculo do valor de polinómios num dado ponto x:
Ler/Introduzir valor de n, a0, a1, . . ., an
Ler/Introduzir valor de x
Inicializar p = a0
Ciclo de i = 1 até n
p = p + ai× xi (1 flop soma + 1 flop produto + (i− 1) flops potência)
Repetir ciclo
Escrever p
Logo o número de flops é (2+ 3+ 4+ . . .+ (n+ 2)) = n[(n+ 2) + 2]/2 = n2/2+ 2n.
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 25
Algoritmo 2 - Técnica de Horner para cálculo do valor de um polinómio num
dado ponto x:
Ler/Introduzir valor de n, a0, a1, . . ., an
Ler/Introduzir valor de x
Inicializar p = an
Ciclo de i = n - 1 até 0
p = ai + p × x (1 flop soma + 1 flop produto)
Repetir ciclo
Escrever p
Logo o número de flops é 2[(n− 1) + 1] = 2n, o que torna o algoritmo 1 «impró-
prio»
em termos de métodos computacionais, porque sendo o número de operações aritmé-
ticas maior, considera-se também menos preciso.
Algoritmo 3 - Algoritmo de Horner para polinómios na forma de Newton com
centros nos nós:
Ler/Introduzir valor de n, a0, a1, . . ., an, c1, . . ., cn
Ler/Introduzir valor de x
Inicializar p = an
Ciclo de i = n - 1 até 0
p = ai + p × (x - ci+1) (2 flops soma + 1 flop produto)
Repetir ciclo
Escrever p
Logo o número de flops é 3[(n− 1) + 1] = 3n. ¤
RESOLUÇÃO
a) O polinómio interpolador na forma de Newton determina-se recorrendo à tabela
de diferenças divididas
x y[·] y[·, ·] y[·, ·, ·] y[·, ·, ·, ·]
0 −1
2
1 1 1/2
3 −1
2 4 −5/2
−2
3 2
p3(x) = −1 + 2x+ 12x(x− 1)− x(x− 1)(x− 2)
b) A expressão para determinar o polinómio interpolador na forma de Lagrange é
pn (x) =
n∑
k=0
Lk (x) yk com Lk =
n∏
i=0 i 6=k
x− xi
xk − xi
26 SEMANA 4
pelo que
L0 =
(
x− 1
−1
)(
x− 2
−2
)(
x− 3
−3
)
L1 =
x
2
(x− 2)(x− 3)
L2 = −
x
2
(x− 1) (x− 3)
L3 =
x
6
(x− 1) (x− 2)
obtendo-se
p3(x) =
(x− 1)(x− 2)(x− 3)
6
+
x
2
(x− 2)(x− 3)
− 2x(x− 1)(x− 3) + x
3
(x− 1)(x− 2)
c) Aplica-se aqui o Algoritmo 3 ao polinómio da alínea a)
p3 (x) = −1 + 2x+
1
2
x (x− 1)− x (x− 1) (x− 2)
Logo n = 3 e x = 4. Precisamos também de C1 = 0, C2 = 1 e C3 = 2. Os centros
contam-se de 1 a n, não existindo C0.
Seguindo o algoritmo, obtém-se sucessivamente:
p = a3 = −1
p = a2 + p(x− C3) = (1/2) + (−1)(4− 2) = −(3/2)
p = a1 + p(x− C2) = 2− (3/2)(4− 1) = −(5/2)
p = a0 + p(x− C1) = −1− (5/2)(4− 0) = −11
Pelo que, o valor de p em x = 4 é p(4) = −11. ¥
5
TÓPICOS: Interpolação Polinomial: Máximos e mínimos locais, Pontos de
inflexão;Interpolação de Hermite, Nós duplos e triplos e diferença dividida confluente.
LEITURAS RECOMENDADAS: Capítulo 2 (pp. 52–73 excepto secção 2.5
e pp. 77–80) do livro [Pina(1995)]
5.1 INTERPOLAÇÃO POLINOMIAL:
Máximos e mínimos locais, Pontos de inflexão
Considere a tabela de valores de uma função y = f(x):
x 0 1 3 4
y −21 5 −15 35
a) Determine as aproximações no intervalo [0, 4] para a localização de um máximo
local, um mínimo local e um ponto de inflexão da função.
b) Obtenha uma interpolante de f(x) na forma de Newton com centros em x0 =
1/2, x1 = 1 e x2 = 3/2.
c) Se a tabela tivesse um maior número de pontos, por exemplo superior a 5, o
método da alínea a) seria recomendável? Justifique.
RESOLUÇÃO
a) As aproximações são aqui obtidas por interpolação polinomial na forma de
Newton, com centros nos nós.
x y[·] y[·, ·] y[·, ·, ·] y[·, ·, ·, ·]
0 −21
26
1 5 −12
−10 8
3 −15 20
50
4 35
O polinómio obtido é
p3(x) = −21 + 26(x− 0)− 12(x− 0)(x− 1) + 8(x− 0)(x− 1)(x− 3)
= −21 + 62x− 44x2 + 8x3
28 SEMANA 5
Há que ter cuidado em colocar como coeficientes do polinómio os elementos da
primeira diagonal (a negro), na ordem indicada.
Por exemplo, se escolhêssemos os elementos da diagonal inferior 35, 50, 20, 8 obte-
ríamos o polinómio 35 + 50(x− 4) + 20(x− 3)(x− 4) + 8(x− 3)(x− 4)(x− 1) que é
uma interpolante correcta (Porquê? Que centros foram agora utilizados?).
Deste modo, se pretendermos que os centros sejam 1, 3 e 4, deixando de fora o 0,
teremos de construir uma nova tabela de diferenças divididas, tal que
x y[·] y[·, ·] y[·, ·, ·] y[·, ·, ·, ·]
1 5
−10
3 −15 20
50 8
4 35 12
14
0 −21
obtendo-se o polinómio:
p3(x) = 5−10(x−1)+20(x−1)(x−3)+8(x−1)(x−3)(x−4), que após expansão na
forma de potências simples fica p3(x) = −21+62x− 44x2+8x3. Como esperávamos,
o polinómio é o mesmo.
Podemos verificar rapidamente que as condições de interpolação são verificadas
para o polinómio obtido, pois, por exemplo, p3(0) = −21 e p3(1) = 70− 65 = 5.
0 0.5 1 1.5 2 2.5 3 3.5 4
−30
−20
−10
0
10
20
30
40
Figura 5.1.1: Polinómio interpolador
As primeira e segunda derivadas do polinómio interpolador são:
p′3(x) = 62− 88x+ 24x2, p′′3(x) = −88 + 48x
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 29
Donde se conclui que os pontos de estacionaridade são:
p′3(x) = 62− 88x+ 24x2 = 0⇔ x =
88±√882 − 4× 24× 62
2× 24
⇔ x1 = 0.95142 ∨ x2 = 2.7152
O ponto de inflexão é:
p′′3(x) = −88 + 48x = 0⇔ x = 1.8333
Como p3 é de grau 3 e p′′3 é negativa para x < 1.8333 e positiva para x > 1.8333
então 0.95142 corresponde a um máximo local e 2.7152 a um mínimo local.
b) Na alínea a) obteve-se o polinómio na forma de Newton com centros nos nós
da tabela. Para considerar a sua representação nos nós 1/2, 1 e 3/2 há que proceder
a uma mudança de nós.
Alternativa 1 - Aplicar um algoritmo de mudança de nós.
Algoritmo de Horner para mudança de centros de um polinómio (Algoritmo 2.2.2,
p. 42 de [Pina(1995)]):†
Introduzir o centro c e retirar o centro cn
Ler/Introduzir valor de n, a0, a1, . . ., an, c1, . . ., cn
Inicializar a
′
n = an
Ciclo de i = n - 1 até 0
a
′
i = ai + a
′
i+1 × (c - ci+1) (2 flops soma + 1 flop produto)
Repetir ciclo
Escrever a
′
Utilize-se a representação com centros na origem:
p3(x) = −21 + 62x− 44x2 + 8x3
e introduzam-se sucessivamente os centros pretendidos. Tal é feito introduzindo estes
centros por ordem inversa à definida.
Introdução do centro c = 3/2 (n = 3)
a3 = 8 a′3 = 8
a2 = −44 c3 = 0 c = 3/2 −44 + 8× (3/2− 0) a′2 = −32
a1 = +62 c2 = 0 c = 3/2 +62− 32× (3/2− 0) a′1 = 14
a0 = −21 c1 = 0 c = 3/2 −21 + 14× (3/2− 0) a′0 = 0
p3(x) = 0 + 14 (x− 3/2)− 32 (x− 3/2) (x− 0) + 8 (x− 3/2) (x− 0)2
†O polinómio inicial deve estar na forma de potências simples.
30 SEMANA 5
Introdução do centro c = 1 (n = 3)
a3 = 8 a′3 = 8
a2 = −32 c3 = 0 c = 1 −32 + 8× (1− 0) a′2 = −24
a1 = +14 c2 = 0 c = 1 +14− 24× (1− 0) a′1 = −10
a0 = 0 c1 = 3/2 c = 1 0− 10× (1− 3/2) a′0 = 5
p3(x) = 5− 10(x− 1)− 24(x− 1)(x− 3/2) + 8(x− 1)(x− 3/2)(x− 0)
Introdução do centro c = 1/2 (n = 3)
a3 = 8 a′3 = 8
a2 = −24 c3 = 0 c = 1/2 −24 + 8× (1/2− 0) a′2 = −20
a1 = −10 c2 = 3/2 c = 1/2 −10− 20× (1/2− 3/2) a′1 = 10
a0 = 5 c1 = 1 c = 1/2 5 + 10× (1/2− 1) a′0 = 0
p3(x) = 0 + 10(x− 1/2)− 20(x− 1/2)(x− 1) + 8(x− 1/2)(x− 1)(x− 3/2)
Alternativa 2 - Como a representação polinomial é única, podemos fazer:
p3(x) = −21 + 62x− 44x2 + 8x3
p3(0.5) = −21 + 62× 0.5− 44× 0.25 + 8× 0.125 = −21 + 31− 11 + 1 = 0
p3(1.0) = −21 + 62− 44 + 8 = 5
p3(1.5) = −21 + 62× 1.5− 44× 1.52 + 8× 1.53 = 0
Utilizamos, então, estes pontos para construir a nova representação, usando como
ponto auxiliar x = 3:
x y[·] y[·, ·] y[·, ·, ·] y[·, ·, ·, ·]
0.5 0
10
1 5 −20
−10 8
1.5 0 0
−10
3 −15
Obtemos exactamente a representação esperada devido à unicidade acima referida:
p3(x) = 0 + 10(x− 1/2)− 20(x− 1/2)(x− 1) + 8(x− 1/2)(x− 1)(x− 3/2)
Omesmo sucede com o ponto auxiliar x = 4, pois, neste caso, a tabela de diferenças
divididas é
x y[·] y[·, ·] y[·, ·, ·] y[·, ·, ·, ·]
0.5 0
10
1 5 −20
−10 8
1.5 0 8
14
4 35
ou o ponto auxiliar x = 0:
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 31
x y[·] y[·, ·] y[·, ·, ·] y[·, ·, ·, ·]
0.5 0
10
1 5 −20
−10 8
1.5 0 −24
14
0 −21
d) Se a tabela tivesse um maior número de pontos, por exemplo superior a 5, o
método da alínea a) não seria recomendável porque:
- obteríamos uma equação polinomial de elevado grau cuja derivada seria morosa
de resolver.
- as oscilações são maiores em polinómios interpoladores de elevado grau.
Nesse caso seria preferível utilizar um método de diferenças finitas apropriado
seguido de interpolação. ¥
5.2 INTERPOLAÇÃO POLINOMIAL:
Interpolação de Hermite, Nós duplos e diferença
dividida confluente
a) Construa o polinómio de menor grau que interpola a função f(x) = ln(1 + x) e a
sua primeira derivada f ′(x) nos extremos do intervalo [0, 1].
b) Majore o erro cometido.
RESOLUÇÃO
a) Interpolar com valores de f e f ′ em cada nó corresponde à interpolação de
Hermite, a qual se obtém aplicando as expressões do Exemplo 2.6.1, pp. 71–72 de
[Pina(1995)], [U0, U1, V0 e V1] com f(0) = 0, f(1) = ln(2), f ′(0) = 1 e f ′(1) = 1/2.
Contudo, podemos também utilizar a tabela de diferenças divididas com a técnica
de nós múltiplos. De facto, da definição de diferença dividida podemos constatar que
f [x0, x1, ..., xk] =
1
k!
f (k) (ξ)→ f [x0, x1] = f (x0)− f (x1)(x0 − x1) ,
obtendo-se a diferença dividida confluente
lim
x1→x0
f [x0, x1] = f [x0, x0] = f ′ (x0) .
De igual modo, obtemos
f [x0, x0, x0] =
f ′′ (x0)
2!
.
Podemos então construir a tabela seguinte:
32 SEMANA 5
x y[·] y[·, ·] y[·, ·, ·] y[·, ·, ·, ·]
0 0
1
0 0 ln(2)− 1
ln(2) 3/2− 2ln(2)
1 ln(2) 1/2− ln(2)
1/2
1 ln(2)
p3(x) = x+ [ln(2)− 1]x2 +
[
3
2
− 2 ln(2)
]
x2(x− 1)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
(a) Polinómio
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.5
0.55
0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
1
(b) Primeira derivada
Figura 5.2.1: Polinómio interpolador e respectiva primeira derivada
b) Análise de Erros de Interpolação - Quanto é que no máximo um polinómio
interpolador se afasta da respectiva função.
Pelo Teorema 2.4.1, p. 55, [Pina(1995)], sabemos que
en(x) ≡ f(x)− pn(x) = 1(n+ 1)!f
(n+1) (ξ)Wn (x)
Com base neste Teorema obtemos o majorante do erro dado por (2.4.3):
‖en‖∞ ≤ 1(n+ 1)!‖f
(n+1)‖∞‖Wn‖∞
A estimativa do erro (2.4.4)
‖en‖∞ ≤ 14 (n+ 1)‖f
(n+1)‖∞hn+1 = 14 (3 + 1)6× 1 =
6
16
=
3
8
= 0.3750,
não deve ser utilizada quando se pede um majorante, porque, tratando-se de uma
interpolação da função e da derivada, i.e., uma interpolação de Hermite,considera-se
o número de pontos apenas 2, logo n = 1, e o grau do polinómio é 2n + 1, podendo
mostrar-se que o erro é dado, neste caso, por (p. 73, [Pina(1995)]):
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 33
en (x) =
(Wn(x))
2
(2n+ 2)!
f (2n+2) (ξ) ≤
max
[xmin,xmax]
W 2n (x)
(2n+ 2)!
‖f (2n+2)‖∞.
Concretizando para o caso presente, obtemos
e (x) =
W 21 (x)
(2 + 2)!
‖f (2+2) (ξ)‖∞ ≤
max
x∈[0,1]
[
(x− 0)2 (x− 1)2
]
24
∥∥∥[−6 (1 + x)−4]∥∥∥
∞
=
1
16
24
6 =
1
64
= 0.15625× 10−1,
pois
f(x) = ln(1 + x)
f ′(x) = (1 + x)−1
...
f (iv) = −6(1 + x)−4
sendo
∥∥f (iv)∥∥ monótona decrescente em [−1,+∞] e, logo, no intervalo [0, 1] tem má-
ximo em x = 0. ¥
5.3 INTERPOLAÇÃO POLINOMIAL:
Interpolação de Hermite, Nós triplos e diferença
dividida confluente
Construa o polinómio de menor grau que satisfaz as seguintes condições de interpo-
lação: p(1) = 2, p′(1) = 0, p′′(1) = 2, p(2) = 2 e p(3) = 3.
RESOLUÇÃO
a) Tendo por base a definição de diferença dividida:
f [x0, x1, ..., xk] =
1
k!
f (k) (ξ) ,
verificamos que
f [x0, x0, x0] =
f ′′ (x0)
2!
.
Utilizando a tabela de diferenças finitas com a técnica de nós múltiplos (neste
caso, nó triplo):
34 SEMANA 5
x y[·] y[·, ·] y[·, ·, ·] y[·, ·, ·, ·] y[·, ·, ·, ·, ·]
1 2
0
1 2 1
0 −1
1 2 0 5/8
0 1/4
2 2 1/2
1
3 3
obtemos a representação
p4(x) = 2 + (x− 1)2 − (x− 1)3 + 58(x− 1)
3(x− 2) ¥
Como suplemento à resolução do problema, mostramos o gráfico do polinómio e
as suas primeira e segunda derivadas, obtidas analiticamente (Figuras 5.3.1 e 5.3.2 ).
Observamos que as condições de interpolação são respeitadas.
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
1.8
2
2.2
2.4
2.6
2.8
3
Figura 5.3.1: Polinómio interpolador
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
−0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
(a) Primeira derivada
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
−2
0
2
4
6
8
10
12
14
(b) Segunda derivada
Figura 5.3.2: Derivadas analíticas do polinómio interpolador ¡
6
TÓPICOS: Interpolação Polinomial: Spline quadrático; Spline cúbico; Erros
em splines cúbicos.
LEITURAS RECOMENDADAS: Capítulo 2 (pp. 77–87) do livro [Pina(1995)]
6.1 INTERPOLAÇÃO POLINOMIAL:
Spline quadrático
Construa o spline quadrático para f(x) = sin(pi/2× x) em [0, 1] com malha uniforme
de 4 nós. Determine uma aproximação para f ′(0.5) e f ′′(0.5). Obtenha os respectivos
erros absoluto e relativo.
CONCEITOS TEÓRICOS
Os splines correspondem a aproximações de funções por troços.† No caso de
um spline quadrático esta aproximação é obtida entre cada dois pontos consecutivos
através da expressão para cada subintervalo ou, mais prosaicamente, troço (Exp.
(2.7.3), p. 80, [Pina(1995)]):
x ∈ [xi−1, xi], i = 1, 2, . . . , n

Si(x) = yi−1 +mi−1(x− xi−1) + Mi2 (x− xi−1)2
Mi =
mi−mi−1
hi
hi = xi − xi−1
mi = 2
(
yi−1−yi
hi
)
−mi−1
¤
RESOLUÇÃO
Construção do spline:
Neste caso, f é contínua e continuamente diferenciável.
Uma malha uniforme, isto é, nós igualmente espaçados, com quatro pontos no
intervalo [0, 1] terá os nós 0, 1/3, 2/3 e 1.
†Para uma definição mais rigorosa de splines ver p. 77 de [Pina(1995)].
36 SEMANA 6
Podemos calcular os momentos mi se conhecermos o valor de um deles. Em geral
m0 define-se comom0 = f ′(x0). Neste problema temosm0 = f ′(x0) = pi/2×cos(pi/2×
x0) = pi/2 = 1.5708.
Obtemos os restantes mi aplicando a fórmula recorrente mi = 2
(
yi−1−yi
hi
)
−mi−1.
i xi−1 yi−1 xi yi mi−1 mi Mi
1 0 0 13 0.5 1.5708 1.4292 −0.42478
2 13 0.5
2
3 0.86602 1.4292 0.76695 −1.9868
3 23 0.86602 1 1 0.76695 0.036899 −2.1901
Então, ficamos com o spline:
x ∈ [0, 1/3] : S1 = 0 + 1.5708(x− 0) + −0.424782 (x− 0)
2
= 1.5708x− 0.2124x2
x ∈ [1/3, 2/3] : S2 = 0.5 + 1.4292
(
x− 1
3
)
+
−1.9868
2
(
x− 1
3
)2
= 0.5 + 1.4292
(
x− 1
3
)
− 0.99340
(
x− 1
3
)2
x ∈ [2/3, 1] : S3 = 0.86602 + 0.76695
(
x− 2
3
)
+
−2.1901
2
(
x− 2
3
)2
= 0.86602 + 0.76695
(
x− 2
3
)
− 1.0950
(
x− 2
3
)2
A Figura mostra graficamente o spline.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figura 6.1.1: Spline constituído por três troços
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 37
Determinação de uma aproximação para f ′(x):
E′(0.5) = f ′(0.5)− S′(0.5) = f ′(0.5)− S′2(0.5)
= pi/2× cos(pi/2× 1/2)− 1.0981 = 0.012621
E′′(0.5) = f ′′(0.5)− S′′(0.5) = f ′′(0.5)− S′′2 (0.5)
= −(pi/2)2 sin(pi/2× 1/2) + 1.9868 = 0.24208
e′(0.5) = E′(0.5)/f ′(0.5) = 0.012621/(pi/2× cos(pi/2× 1/2)) = 1.1363× 10−2,
ou seja, 1.14%.
e′′(0.5) = E′′(0.5)/f ′′(0.5) = 0.24208/(−(pi/2)2 sin(pi/2× 1/2)) = −0.13875,
ou seja, 13.9%.
f(x) e S(x) são praticamente coincidentes.
Justifica-se um erro maior em S′′(0.5) do que em S′(0.5), porque nos splines
quadráticos S′(x) é linear e S′′(x) é constante por troços (Ver Figura 6.1.2).
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
(a) Primeira derivada
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
−2.5
−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
(b) Segunda derivada
Figura 6.1.2: Derivadas dos três troços do spline ¥
6.2 INTERPOLAÇÃO POLINOMIAL:
Spline cúbico
Construir o spline cúbico que interpole os pontos: (0, 0); (0.5, 0.0625); (1, 1); (1.5, 5.0625)
e (2, 16).
CONCEITOS TEÓRICOS
Os splines correspondem a aproximações de funções por troços.† No caso de um
spline cúbico esta aproximação é obtida entre cada dois pontos consecutivos através
da expressão para cada subintervalo ou, mais prosaicamente, troço (Exp. (2.7.10),
†Para uma definição mais rigorosa de splines ver p. 77 de [Pina(1995)].
38 SEMANA 6
p. 82, [Pina(1995)]):
x ∈ [xi−1, xi], i = 1, 2, . . . , n

Si(x) =Mi−1
(xi−x)3
6hi
+Mi
(x−xi−1)3
6hi
+
(
yi−1 −Mi−1 h
2
i
6
)
xi−x
hi
+
(
yi −Mi h
2
i
6
)
x−xi−1
hi
hi = xi − xi−1
Mi = S′′i (xi) = S
′′(xi)
Os Mi, designados momentos, são dados pela exp. (2.7.13), p. 82 de [Pina(1995)]
obtida impondo S′(xi−) = S′(xi+), isto é, impondo continuidade da primeira deri-
vada nos nós:
hi
6
Mi−1 +
hi + hi+1
3
Mi +
hi+1
6
Mi+1 =
yi+1 − yi
hi+1
− yi − yi−1
hi
,
com i = 1, 2, . . . , n− 1
Obtemos então n − 1 equações com n + 1 incógnitas Mi, i. e., M0, . . . ,Mn, pelo
que são necessárias duas condições suplementares. Há várias possibilidades, embora
as três a seguir descritas sejam as mais usadas:
1) Spline completo: S′1(x0) = y′0, S′n(xn) = y′n.
2) Spline natural: S′′1 (x0) =M0 = 0, S′′n(xn) =Mn = 0.
Este tipo de spline é indicado quando não conhecemos as derivadas y′0 e y′n. No en-
tanto, estas condições podem reduzir a precisão do spline.
3) Spline periódico: y0 = yn, S′(x0) = S′(xn), M0 =Mn. ¤
RESOLUÇÃO
Construção do spline:
Neste caso vamos recorrer ao spline natural, n = 4 e hi = h = 0.5. Como o
espaçamento h é constante, a expressão simplifica-se e com i = 1, 2, 3, ficamos com:
M0 = 0
1
2
M0 + 2M1 +
1
2
M2 = 3
(
y2 − 2y1 + y0
h2
)
1
2
M1 + 2M2 +
1
2
M3 = 3
(
y3 − 2y2 + y1
h2
)
1
2
M2 + 2M3 +
1
2
M4 = 3
(
y4 − 2y3 + y2
h2
)
M4 = 0
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 39
ou, na forma matricial,
1 0 0 0 0
1/2 2 1/2 0 0
0 1/2 2 1/2 0
0 0 1/2 2 1/2
0 0 0 0 1


M0
M1
M2
M3
M4

=

0
3
[
(y2 − 2y1 + y0)/h2
]
3
[
(y3 − 2y2 + y1)/h2
]
3
[
(y4 − 2y3 + y2)/h2
]
0

Substituindoos valores de yi, i = 0, 1, 2, 3, 4 e eliminando as primeiras e últimas
linhas e colunas obtemos sucessivamente,†
2 1/2 0 10.5
1/2 2 1/2 37.5
0 1/2 2 82.5

→

2 0.5 0 10.5
0 1.875 0.5 34.875
0 0.5 2 82.5

→
→

2 0.5 0 10.5
0 1.875 0.5 34.875
0 0 1.8667 73.2

Portanto, os momentos são
M3 = 73.2/1.8667 = 39.214
M2 = (34.875− 0.5M3)/1.875 = 8.1429
M1 = (10.5− 0.5M2)/2 = 3.2143
†O leitor poderá consultar o Capítulo 12 referente a Sistemas de Equações Lineares para compre-
ender a obtenção dos valores M1, M2 e M3.
40 SEMANA 6
O primeiro troço do spline corresponde a x ∈ [x0, x1], sendo
S1(x) =M0
(x1 − x)3
6h
+M1
(x− x0)3
6h
+
(
y0 −M0h
2
6
)
x1 − x
h
+
(
y1 −M1h
2
6
)
x− x0
h
= 0 + 3.2143
(x− 0)3
6× 12
+
(0− 0) 0.5− x1
2
+(
0.0625− 3.2143(
1
2 )
2
6
)
x− 0
1
2
,
ou seja,
S1(x) = 1.0714x3 − 0.14286x para 0 6 x 6 0.5.
Os restantes troços determinam-se da mesma forma, obtendo-se
S2 = 1.6428x3 − 0.85714x2 + 0.2857x− 0.071428 para 0.5 6 x 6 1
S3 = 10.357x3 − 26.999x2 + 26.428x− 8.7857 para 1 6 x 6 1.5
S4 = −13.071x3 − 78.428x2 − 131.71x+ 70.286 para 1.5 6 x 6 2
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
−2
0
2
4
6
8
10
12
14
16
Figura 6.2.1: Spline cúbico natural resultante
¥
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 41
6.3 INTERPOLAÇÃO POLINOMIAL:
Erros em splines cúbicos
Pretende-se interpolar por um spline cúbico numa malha uniforme no intervalo [0, 1]
a função f(x) = exp(1 + x). Determine o número de nós necessário para obter um
erro relativo em valor absoluto inferior a 10−3.
CONCEITOS TEÓRICOS
O cálculo de erros de interpolação com splines cúbicos tem por base o Teorema
2.7.4, p. 86, [Pina(1995)]:
Teorema 2.7.4: Seja f ∈ C4(Ω¯) e S o spline cúbico satisfazendo qualquer das
condições suplementares referidas nesta subsecção [pp. 82–83, [Pina(1995)]]. Então,
‖f − S‖∞ 6
5
384
∥∥D4f∥∥∞ h4
‖D(f − S)‖∞ 6
(√
3
216
+
1
24
)∥∥D4f∥∥∞ h3∥∥D2(f − S)∥∥∞ 6 ( 112 + 13(h/h)
)∥∥D4f∥∥∞ h2∥∥D3(f − S)∥∥∞ 6 12 (1 + (h/h)2) ∥∥D4f∥∥∞ h
em que
h = min
16i6n
hi ¤
RESOLUÇÃO
O erro relativo em valor absoluto é tal que
|e(x)| =
∣∣∣∣f(x)− S(x)f(x)
∣∣∣∣ 6 maj |f(x)− S(x)|min |f(x)|
sendo maj |f(x)− S(x)| ≡ ‖f − S‖∞.
Pelo Teorema 2.7.4, fica
5
384
∥∥D4f∥∥∞ h4
min |f(x)| 6 ² = 10
−3 →
5
384 ‖exp(1 + x)‖∞ h4
min | exp(1 + x)| 6 ² = 10
−3
⇔
5
384 maxx∈[0,1]
|exp(1 + x)|h4
min
x∈[0,1]
|exp(1 + x)| 6 ² = 10
−3
⇔
5
384 exp(2)h
4
exp(1)
6 ² = 10−3,
pois,
∣∣D4f ∣∣ = exp(1 + x) é monótona crescente no intervalo [0,+∞] e, logo, tem um
mínimo em x = 0 e um máximo em x = 1 no intervalo [0, 1].
42 SEMANA 6
Assim, resulta
h <
(
10−3 exp(1)
5
384 exp(2)
)1/4
h < 0.40998→ h = b− a
N
=
1− 0
N
< 0.40998→ N > 2.43914,
ou seja, 3 intervalos. Assim, são necessários pelo menos 4 nós. ¥
7
TÓPICOS: Diferenciação Numérica: Diferenças finitas de primeira e segunda
ordens; Majorante do erro; Espaçamento óptimo; Ponto de inflexão; Espaçamento
desigual.
LEITURAS RECOMENDADAS: Capítulo 3 (pp. 97–111) do livro [Pina(1995)]
7.1 DIFERENCIAÇÃO NUMÉRICA:
Diferenças finitas de primeira e segunda ordens
Considere a seguinte tabela de valores para a função y = coshx:
x 0.1 0.2 0.3 0.4
y 1.0050 1.0201 1.0453 1.0811
Determine os valores de f ′(0.2) e f ′′(0.2) pelas várias fórmulas de diferenças finitas
e compare os resultados com os valores exactos.
CONCEITOS TEÓRICOS
Consultar as Secções 3.2 e 3.3 de [Pina(1995)]. ¤
−4 −3 −2 −1 0 1 2 3 4
−15
−10
−5
0
5
10
15
Figura 7.1.1: Funções cosh(x) (—) e sinh(x) (- - -)
44 SEMANA 7
RESOLUÇÃO
Valores exactos:
y = coshx =
ex + e−x
2
, y′ = sinhx =
ex − e−x
2
, y′′ = coshx =
ex + e−x
2
= y
y′(0.2) = sinh(0.2) = 0.20134, y′′(0.2) = cosh(0.2) = 1.0201
Valores por diferenças finitas:
Para o cálculo da primeira derivada, vamos usar diferenças finitas de primeira e
segunda ordens.
Diferenças finitas de 1.a ordem Dhf(x) com h = 0.1:
Progressiva (exp. (3.2.8)) f(x+h)−f(x)h =
f(0.2+0.1)−f(0.2)
0.1 = 0.2520
Regressiva (exp. (3.2.10)) f(x)−f(x−h)h =
f(0.2)−f(0.2−0.1)
0.1 = 0.1510
Central (exp. (3.2.12)) f(x+h)−f(x−h)2h =
f(0.2+0.1)−f(0.2−0.1)
2×0.1 = 0.2015
Diferenças finitas de 2.a ordem Dhf(x) com h = 0.1:
Progressiva (Exp. (3.2.16)) −3f(x)+4f(x+h)−f(x+2h)2h = ... = 0.1990
Regressiva (Exp. (3.2.20)) 3f(x)−4f(x−h)+f(x−2h)2h = ... = 0.20150
Central (Exp. (3.2.12)) Expressão e valor iguais a 1.a ordem
Para calcular a diferença finita regressiva de 2.a ordem precisamos de saber o valor
da função em x = 0, que neste caso é f(0) = cosh(0) = 1.
A Figura abaixo mostra as tangentes, isto é, polinómios de ordem um, dados pelas
diferenças finitas de primeira ordem.
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 45
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
0.96
0.98
1
1.02
1.04
1.06
1.08
1.1
R
P
C
Figura 7.1.2: Diferenças Finitas: (P) Progressiva, (R) regressiva e (C) central
A melhor aproximação foi obtida com a Diferença Finita Central, e com a Dife-
rença Regressiva de 2.a ordem, mas esta última requer um maior número de pontos e
inclusive um valor de f(x) que não está na tabela dada.
Para a segunda derivada, vamos considerar apenas a diferença finita central (ob-
viamente de segunda ordem).
Exp. (3.3.1) com h = 0.1:
D2hf(x) =
f(x+ h)− 2f(x) + f(x− h)
h2
=
f(0.3)− 2f(0.2) + f(0.1)
0.12
=
1.0453− 2× 1.0201 + 1.005
0.12
= 1.0100,
ou seja, sendo o valor exacto y′′(0.2) = 1.0201, esta aproximação dispõe de dois dígitos
correctos. ¥
7.2 DIFERENCIAÇÃO NUMÉRICA:
Majorante do erro
Obtenha o valor de f ′(1.0) para a função f(x) = exp(−x) usando diferenças finitas
progressivas e passo h = 0.001. Determine o erro efectivamente cometido e compare-o
com o majorante teórico.
CONCEITOS TEÓRICOS
Consultar as Secções 3.2 e 3.3 de [Pina(1995)]. ¤
46 SEMANA 7
RESOLUÇÃO
Valores aproximado e exacto:
Dhf(x) =
f(x+ h)− f(x)
h
=
f(1.001)− f(1.0)
0.001
=
e−1.001 − e−1
0.001
= −0.3676956
f ′(x) = −ex ⇒ f ′(1.0) = −e−1.0 = −0.3678794
Erro efectivamente cometido:
E(1.0) = f ′(1.0)−Dhf(1.0) = −0.3678794− (−0.3676956)
= −0.0001838 = −(1.838)10−4
O majorante teórico do erro (erro devido apenas à aproximação da derivada) é
dado pela exp. (3.2.9), p. 100 de [Pina(1995)]: †
e′1(x) = −
1
2
hf ′′(ξ) em que ξ ∈ [x0, x1] = Ω¯
|e′1(x)| 6
1
2
h ‖f ′′(x)‖∞
Neste problema, f ′′(x) = e−x em Ω¯ = [1.0, 1.0 + h]. Concretizando, e atendendo
a que a função é monótona, o máximo encontra-se num extremo do intervalo (ver
Figura 7.2.1), fica
‖f ′′(x)‖Ω¯∞ = f ′′(1.0) = f(1.0) = 0.367879
|e′1(1.0)| 6
1
2
× 0.001× 0.367879 = (1.839)10−4
Constatamos que, de facto, |E(1.0)| 6 |e′1(1.0)|, sendo praticamente igual.
†Utiliza-se para o erro de interpolação a letra minúscula e de forma a manter a coerência com a
nomenclatura utilizada no livro de [Pina(1995)]. Note-se no entanto que se trata de um erro absoluto
e por isso deveríamos utilizar uma letra maiúscula.
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 47
−1.5 −1 −0.5 0 0.5 1 1.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
Figura 7.2.1: Funções exp(−x) (—) e exp(x) (- - -) ¥
7.3 DIFERENCIAÇÃO NUMÉRICA:
Espaçamento óptimo
Dada a seguinte tabela de tempos e posições:
t 0 1/3 2/3 1
y(t) 0 sen(pi/6) sen(pi/3) 1
a) Determine uma aproximação da velocidade e da aceleração nos 2 nós interiores
usando diferenças finitas centrais.
b) Deduza uma estimativa para o espaçamento óptimo no cálculo da aceleração
por diferenças finitas centrais.
RESOLUÇÃO
a) Antes de mais, devemos observarque h = 1/3 e que os 2 pontos interiores são
t = 1/3 e t = 2/3.
Expressão para cálculo das velocidades por diferenças finitas centrais (p. 101, exp.
(3.2.12)):
Dhf(x) =
f(x+ h)− f(x− h)
2h
Expressão para cálculo das acelerações por diferenças finitas centrais (p. 103, exp.
(3.3.1)):
D2hf(x) =
f(x− h)− 2f(x) + f(x+ h)
h2
48 SEMANA 7
Concretizando para os valores dados na tabela e recorrendo à calculadora, obte-
mos, sucessivamente:
Dhf(1/3) =
f(1/3 + 1/3)− f(1/3− 1/3)
2× 1/3 =
f(2/3)− f(0)
2× 1/3
=
sen(pi/3)− 0
2× 1/3 = 1.299
Dhf(2/3) =
f(2/3 + 1/3)− f(2/3− 1/3)
2× 1/3 =
f(1)− f(1/3)
2× 1/3
=
1− sen(pi/6)
2× 1/3 = 0.75
D2hf(1/3) =
f(1/3− 1/3)− 2f(1/3) + f(1/3 + 1/3)
(1/3)2
=
f(0)− 2f(1/3) + f(2/3)
(1/3)2
=
0− 2sen(pi/6) + sen(pi/3)
(1/3)2
= −1.2057
D2hf(2/3) =
f(2/3− 1/3)− 2f(2/3) + f(2/3 + 1/3)
(1/3)2
=
f(1/3)− 2f(2/3) + f(1)
(1/3)2
=
sen(pi/6)− 2sen(pi/3) + 1
(1/3)2
= −2.088457
b) Neste caso, o espaçamento h óptimo é desconhecido.
O valor exacto da segunda derivada f ′′(x) relaciona-se com o valor por diferenças
finitas centrais e o erro de aproximação, e′′2(x), inerente ao cálculo por este processo,
do seguinte modo:
f ′′(x) =
f(x− h)− 2f(h) + f(x+ h)
h2
− 1
12
h2f (4)(η)
=
f(x− h)− 2f(h) + f(x+ h)
h2
− e′′2(x)
Em representação em ponto flutuante, temos de contar ainda com os erros de
arredondamento, e˜(x), no cálculo do valor da própria função, ou seja,
f(x) = f˜(x) + e˜(x) e, portanto
f ′′(x) =
f˜(x− h)− 2f˜(h) + f˜(x+ h)
h2
+
e˜(x− h)− 2e˜(h) + e˜(x+ h)
h2
− 1
12
h2f (4)(η)
=
f˜(x− h)− 2f˜(h) + f˜(x+ h)
h2
+ E
onde
E =
e˜(x− h)− 2e˜(h) + e˜(x+ h)
h2
− 1
12
h2f (4)(η)
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 49
Considerando∥∥∥f (4)∥∥∥
∞
6M4 e |e˜(x− h)− 2e˜(h) + e˜(x+ h)| 6 4²
onde M4 é o majorante das quartas derivadas e ² é um parâmetro que depende da
forma como a função f é calculada, ou seja, é da ordem de grandeza da unidade de
arredondamento do computador utilizado, podemos fazer
|E| 6 4²
h2
+
M4h
2
12
O segundo membro é mínimo quando a sua primeira derivada, em ordem a h, for
nula, isto é, quando
− 8²
h3
+
M4h
6
= 0⇒ h =
(
48²
M4
)1/4
Portanto, o valor de h deduzido será a estimativa para o espaçamento óptimo no
cálculo da aceleração por diferenças finitas centrais. ¥
Como complemento teórico à resolução deste problema, sugerimos a observação
atenta da Figura 3.5.1, na p. 108. do livro [Pina(1995)]. ¡
7.4 DIFERENCIAÇÃO NUMÉRICA:
Ponto de inflexão
Obtenha o ponto de inflexão de uma função f no intervalo [0, 1] sendo conhecidos os
seguintes valores:
x 0.00 0.25 0.50 0.75 1.00
f(x) 0.00 0.56 0.72 1.34 4.00
RESOLUÇÃO
Podemos considerar os seguintes processos para resolver o problema:
i) Determinar o polinómio interpolador e depois derivar analiticamente;
ii) Determinar as derivadas aproximadasD2hf(x) pela tabela e interpolar os valores
obtidos;
iii) Determinar os splines, seguido de derivação analítica.
Em qualquer dos casos obtemos uma equação polinomial a solucionar. De notar
que existem outros processos, igualmente válidos para resolver este problema.
Processo i:
Esquematicamente temos
pn(x)→ p′′n(x)→ p′′n(x) = 0
Assim, primeiro determinamos o polinómio interpolador p4(x) por diferenças di-
vididas:
50 SEMANA 7
x y[·] y[·, ·] y[·, ·, ·] y[·, ·, ·, ·] y[·, ·, ·, ·, ·]
0.00 0.00
2.24
0.25 0.56 −3.2
0.64 9.17(3)
0.50 0.72 3.68 7.68
2.48 16.85(3)
0.75 1.34 16.32
10.64
1.00 4.00
O polinómio é então
p(x) = 2.24x− 3.2x(x− 0.25) + 9.17(3)x(x− 0.25)(x− 0.50)
+ 7.68x(x− 0.25)(x− 0.50)(x− 0.75)
que, desenvolvido e agrupado na forma de potências simples de x corresponde a
p(x) = 3.4(6)x− 4.8x2 − 2.34(6)x3 + 7.68x4.
Derivando sucessivamente, obtemos
p′(x) = 3.4(6)− 9.6x− 7.04x2 + 30.72x3
p′′(x) = −9.6− 14.08x+ 92.16x2
O ponto de inflexão é tal que
p′′(x) = 0⇔ x = −b±
√
b2 − 4ac
2a
ou seja, x = 0.408054 ∨ x = −0.255276
Portanto, o ponto de inflexão no intervalo [0, 1], obtido por este processo, é
x ≈ 0.40805.
Processo ii:
Esquematicamente temos
D2hf(x)→ p′′n(x)→ p′′n(x) = 0
A segunda derivada é obtida por diferenças finitas centrais:
D2hf(x) =
f(x− h)− 2f(x) + f(x+ h)
h2
tal que h = 0.25. Assim, obtemos
x 0.00 0.25 0.50 0.75 1.00
y 0.00 0.56 0.72 1.34 4.00
y′′ −6.4 7.36 32.64
e a correspondente tabela de diferenças divididas para determinar a interpolação de
D2hf(x):
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 51
x y′′[·] y′′[·, ·] y′′[·, ·, ·]
0.25 −6.4
55.04
0.50 7.36 92.16
101.12
0.75 32.64
O polinómio que interpola a segunda derivada é, então
p′′(x) = −6.4 + 55.04(x− 0.25) + 92.16(x− 0.25)(x− 0.50)
que, desenvolvido e agrupado em potências de x corresponde a
p′′(x) = −8.64− 14.08x+ 92.16x2
O ponto de inflexão é tal que
p′′(x) = 0⇔ x = −b±
√
b2 − 4ac
2a
ou seja, x = 0.391960 ∨ x = −0.239182
Portanto, o ponto de inflexão no intervalo [0, 1], obtido por este processo, é
x ≈ 0.39196.
Processo iii:
Esquematicamente temos
Mi → S′′i (x)→ S′′i (x) = 0
Reparar que não necessitamos de determinar o próprio spline cúbico Si(x), bas-
tando, para o que pretendemos, calcular S′′i (x).
Dado que desconhecemos as derivadas y′ nos extremos x0 e x4, temos M0 =M4 =
0, isto é, impomos S′′1 (x0) = 0 =M0 e S′′4 (x4) = 0 =M4.†
Para determinar M1, M3 e M3 faz-se
h
6
Mi−1 +
2h
3
Mi +
h
6
Mi+1 =
yi+1 − 2yi + yi−1
h
com i = 1, 2, 3
1
2
Mi−1 + 2Mi +
1
2
Mi+1 =
3
h2
(yi+1 − 2yi + yi−1) com i = 1, 2, 3
o que origina, na forma matricial, 2 1/2 01/2 2 1/2
0 1/2 2

 M1M2
M3
 = 3
h2
 y2 − 2y1 + y0y3 − 2y2 + y1
y4 − 2y3 + y2
 =
 −19.222.08
97.92

†Como veremos, ao impor estas igualdades, os resultados deterioram-se.
52 SEMANA 7
Resolvendo o sistema de equações, por exemplo, através do método de Gauss e
substituição descendente, obtemos M1M2
M3
 =
 48.61711.3714
−9.9428

Dado que S′′i =Mi, é de esperar que o ponto de inflexão esteja no segundo troço,
pois M2 e M3 têm sinais diferentes.
Então, como a segunda derivada do i-ésimo troço de spline é tal que (Exp. (2.7.7),
p. 81, [Pina(1995)])
S′′i (x) =Mi−1
xi − x
h
+Mi
x− xi−1
h
obtemos para o segundo troço, isto é, para o subintervalo [0.25, 0.5]
S′′2 (x) =M1
x2 − x
h
+M2
x− x1
h
= −9.94280.50− x
0.25
+ 1.3714
x− 0.25
0.25
O ponto de inflexão é tal que
S′′2 (x) = 0⇔ −9.9428
0.50− x
0.25
+ 1.3714
x− 0.25
0.25
= 0⇔ x = 0.469697.
Portanto, o ponto de inflexão no intervalo [0, 1], obtido por este processo, é
x ≈ 0.46970.
Comparação entre os três processos:
Observando as Figuras 7.4.1(a) e 7.4.1(b), que mostram o polinómio obtido pelo
processo i), é possível concluirmos que o ponto de inflexão positivo se encontra na
vizinhança de 0.4. No entanto, existe um outro ponto de inflexão negativo, perto de
0.25.
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−3
−2
−1
0
1
2
3
4
5
(a)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
(b)
Figura 7.4.1: Polinómio de quarta ordem obtido pelo processo i
Os valores do ponto de inflexão obtidos pelos processos i e ii são bastante seme-
lhantes, pois os polinómios de grau 2 obtidos são quase coincidentes (Figura 7.4.2).
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 53
O resultado com o processo iii não será o melhor, dado que impusemos um va-
lor nulo às segundas derivadas nas extremidades do intervalo (S′′i (0) = 0 = M0 e
S′′i (1) = 0 =M4).
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
−20
−10
0
10
20
30
40
50
60
70
Processo i
Processo ii
Processo iii
Figura 7.4.2: Polinómios obtidos por derivaçãoanalítica (Processos i e ii) e por
spline cúbico natural (Processo iii) )
A título comparativo e para reflexão , mostramos o polinómio obtido pelo processo
i e o spline cúbico natural na Figura 7.4.3.† Porque se verificam grandes discrepâncias
nos intervalos [0, 0.25] e [0.75, 1]?
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
3
3.5
4
Processo iii
Processo i
Figura 7.4.3: Polinómios de quarta ordem (Processo i) e spline cúbico (Processo iii)
¥
†Relembramos que não é necessário construir o spline para calcular o ponto de inflexão.
54 SEMANA 7
7.5 DIFERENCIAÇÃO NUMÉRICA:
Espaçamento desigual
Dados 3 pontos, (x1, y1), (x2, y2), (x3, y3), com espaçamento desigual entre as suas
abcissas, deduza uma fórmula para calcular uma aproximação da segunda derivada
num ponto qualquer do intervalo [x1, x3]. Considerando os pontos (1, 1), (2, 4), (5, 5),
utilize a fórmula para calcular y′′(2) e y′′(3).
Em qual dos pontos, x2 ou (x1 + x3)/2, é de esperar um erro absoluto menor?
Justifique (admita que as derivadas de y = f(x) se mantêm na mesma ordem de
grandeza em [x1, x3]).
RESOLUÇÃO
Pretendemos deduzir uma fórmula para calcular D2hf ≈ f ′′(x) = y′′.
De entre as várias possibilidades, a dedução da diferença finita com espaçamento
desigual é a mais simples. De notar que o recurso a splines não é o mais adequado, pois
não se conhece a função e estes requerem o conhecimento de m0 (splines quadráticos)
ou M0 e Mn (splines cúbicos). Assim,
D2h = p
′′
2(x) = 2f [x1, x2, x3] = 2
f [x1, x2]− f [x2, x3]
x1 − x3
= 2
f(x1)−f(x2)
x1−x2 −
f(x2)−f(x3)
x2−x3
x1 − x3
Como D2h não depende de x, temos
f ′′(2) ≡ f ′′(3) ≈ 2
1−4
1−2 − 4−52−5
1− 5 = −
4
3
A expressão do erro de D2h é
e′′2 = f [x1, x2, x3, x, x, x]W2(x)
+ 2f [x1, x2, x3, x, x]W ′2(x)
+ f [x1, x2, x3, x]W ′′2 (x)
onde, usando
f [x1, x2, x3, . . . , xk] =
f (k)(ξ)
k!
vem
e′′2 =
f (5)(r)
120
W2(x)
+
f (4)(ξ)
12
W ′2(x)
+
f (3)(η)
6
W ′′2 (x)
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 55
Como as derivadas mantêm a mesma ordem de grandeza, ficamos com
e′′2 6
M
12
(
1
10
|W2(x)|+ |W ′2(x)|+ 2|W ′′2 (x)|
)
onde
W2 = (x− x1)(x− x2)(x− x3) = (x− 1)(x− 2)(x− 5)
W ′2 = (x− x2)(x− x3) + (x− x1)(x− x3) + (x− x1)(x− x2)
W ′′2 = 2(x− x3) + 2(x− x2) + 2(x− x1) = 6x− 16
Compilando os valores numa tabela, obtemos
x |W2(x)| |W ′2(x)| |W ′′2 (x)| e′′2
x = x2 = 2 0 3 4 1112M
x = (x1 + x3)/2 = 3 4 4 2 8.412M
Concluímos então que é de esperar um erro absoluto menor em x, embora em
termos de majorante a diferença encontrada não seja significativa. ¥
8
TÓPICOS: Integração Numérica: Regras de Newton-Cotes simples e com-
postas; Regra de Gauss-Legendre-Lobatto; Regras de Gauss-Legendre e do trapézio
composta. Erros absoluto e relativo. Estimativa e majorante do erro.
LEITURAS RECOMENDADAS: Capítulo 4 (pp. 115–139) do livro [Pina(1995)]
8.1 INTEGRAÇÃO NUMÉRICA:
Regras do trapézio simples e de Simpson simples
e composta
(a) Calcule, pela regra do trapézio simples e pela regra de Simpson simples, o valor
aproximado do integral
I =
∫ 3
2
exp (1 + cosx)dx.
(b) Idem pela regra de Simpson composta com 2 subintervalos.
(c) Estime o erro absoluto cometido na alínea anterior.
RESOLUÇÃO
(a) As regras referidas correspondem, respectivamente, às expressões (4.2.6) e
(4.2.7), p. 120, de [Pina(1995)]:
TIh(f) =
b− a
2
[f(a) + f(b)]
SIh(f) =
b− a
6
[
f(a) + 4f(
a+ b
2
) + f(b)
]
onde a e b são os limites de integração.
Concretizado para a função f(x) = exp (1 + cosx), com a = 2 e b = 3 obtemos
TIh(f) =
3− 2
2
[f(2) + f(3)] = 1.401495675
SIh(f) =
3− 2
6
[
f(2) + 4f(
2 + 3
2
) + f(3)
]
= 1.280503052
(b)
58 SEMANA 8
Aplicação directa da exp. (4.5.10): A exp. para o cálculo com a regra de
Simpson composta com N subintervalos é (Exp. (4.5.10), p. 138, de [Pina(1995)]):
ScIh(f) =
h
6
[
f(a) + f(b) + 2
N−1∑
i=1
f(ai) + 4
N∑
i=1
f(ai−1 + h/2)
]
onde h = (b− a)/N .
No caso presente, como N = 2 teremos
ScIh(f) =
h
6
[f(a) + f(b) + 2f(a1) + 4(f(a0 + h/2) + f(a1 + h/2))]
com a = 2, b = 3, h = (3−2)/2 = 0.5, a0 = 2 = a, a1 = a0+h = 2.5, a2 = a1+h = b.
Assim, vem, sucessivamente
ScIh(f) =
0.5
6
{f(2) + f(3) + 2f(2.5) + 4[f(2 + 0.5/2) + f(2.5 + 0.5/2)]}
=
0.5
6
{f(2) + f(3) + 2f(2.5) + 4[f(2.25) + f(2.75)]}
= 1.279922906
Resolução alternativa e justificação da exp. (4.5.10): Neste caso conside-
ramos 2 intervalos, usamos a regra de Simpson simples em cada um deles e somamos
os valores obtidos, pois sabemos que
I =
∫ b
a
f(x)dx =
∫ a1
a
f(x)dx+
∫ b
a1
f(x)dx
=
∫ a1
a0
f(x)dx+
∫ a2
a1
f(x)dx =
∫ 2.5
2
f(x)dx+
∫ 3
2.5
f(x)dx.
Portanto, aplicando a regra de Simpson simples a cada um dos integrais acima,
obtemos
ScIh(f) =
2.5− 2
6
[f(2) + 4f(2.25) + f(2.5)]
+
3− 2.5
6
[f(2.5) + 4f(2.75) + f(3)] = 1.279922906
(c) Antes de mais, há que notar que se pretende uma estimativa, pois desconhe-
cemos o valor exacto do integral.
Os erros de integração das regras de Simpson simples e composta são dados, res-
pectivamente pelas expressões (4.2.18), p. 124 e (4.5.11), p. 138, de [Pina(1995)]:
SEh(f) = − 12880f
(4)(ξ)(b− a)5 = I − SIh(f) (8.1.1)
ScEh(f) = −b− a2880 f
(4)(ξ)h4 = I − ScIh(f) onde h = b− a
N
(8.1.2)
A exp. (8.1.2) pode escrever-se da seguinte forma:
ScEh(f) = −b− a2880 f
(4)(ξ)
(
b− a
N
)4
= − 1
N4
1
2880
f (4)(ξ) (b− a)5 (8.1.3)
EXERCÍCIOS DE MATEMÁTICA COMPUTACIONAL 59
pelo que, e como N = 2, fica
I − SIh(f) = − 12880f
(4)(ξ)(b− a)5 = −c (8.1.4)
I − ScIh(f) = − 124
1
2880
f (4)(ξ) (b− a)5 = − 1
16
c (8.1.5)
Subtraindo (4) a (5), membro a membro, fica
I − SIh(f)− I + ScIh(f) = −c+ 124 c ⇔
ScIh(f)− SIh(f) =
(
1
16
− 1
)
c
(8.1.6)
ou seja,
c =
(
1
16
− 1
)−1 [
ScIh(f)− SIh(f)
]
=
(
−15
16
)−1 [
ScIh(f)− SIh(f)
]
= −16
15
[
ScIh(f)− SIh(f)
]
(8.1.7)
Substituindo c em (5), obtemos
I − ScIh(f) = − 116
(
−16
15
)[
ScIh(f)− SIh(f)
]
=
ScIh(f)− SIh(f)
15
(8.1.8)
Então, substituindo os valores calculados nas alíneas anteriores, obtemos uma
estimativa de erro absoluto cometido na alínea b) de 3.9× 10−5.† ¥
A Figura 8.1.1 mostra a negro a diferença entre as áreas abaixo da curva da função
f(x) = exp(1+cos(x)) e dos polinómios de grau 1 (Figura 8.1.1(a): Regra do trapézio
simples) e grau 2 (Figura 8.1.1(b): Regra de Simpson simples). É notório o melhor
desempenho da regra de Simpson que, recordamos, é de ordem 3.
2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3
0
0.5
1
1.5
2
2.5
(a) Regra do trapézio simples
2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3
0
0.5
1
1.5
2
2.5
(b) Regra de Simpson simples
Figura 8.1.1: Diferenças entre integração exacta e com regras do trapézio e Simpson
simples
¡
†Admitimos que a estimativa de f (4) é válida em ambos os casos, sendo o máximo no intervalo
[2, 3].
60 SEMANA 8
8.2 INTEGRAÇÃO NUMÉRICA:
Regra de Gauss-Legendre e de Gauss-Legendre-
Lobatto
Pretende-se calcular o seguinte integral
I =
∫ 2
0
sinh(x)dx.
(a) Obtenha uma aproximação com 2 pontos de Gauss-Legendre.
(b) Idem, com 3 pontos de Gauss-Legendre-Lobatto.
(c) Calcule o erro absoluto e relativo cometido em cada caso.
RESOLUÇÃO
(a) Mudança de coordenadas (exp. (4.1.2), p. 116):
T (ξ) =
a(1− ξ)
2
+
b(1 + ξ)
2
=
0(1− ξ)
2
+
2(1 + ξ)
2
= 1 + ξ
J(ξ) =
b− a
2
=
2
2
= 1
Então,
I =
∫ 2
0
sinh(x)dx =
∫ 1
−1
sinh(1 + ξ)× 1dξ =
∫ 1
−1
sinh(1 + ξ)dξ.
Com 2 pontos, temos as abcissas e os pesos (Tabela 4.4.1, p. 132)
ξ1 = −0.5773502692, A1 = 1,
ξ1 = +0.5773502692, A2 = 1.
e, portanto
I =
∫ 2
0
sinh(x)dx