Calculo Numerico Um Livro Colaborativo V
Calculo Numerico Um Livro Colaborativo V
Calculo Numerico Um Livro Colaborativo V
Um Livro Colaborativo
Versão Python
21 de junho de 2017
Organizadores
ii
Licença
iii
Nota dos organizadores
http://www.ufrgs.br/numerico
iv
Prefácio
v
Sumário
Capa i
Organizadores ii
Licença iii
Prefácio v
Sumário x
1 Introdução 1
vi
SUMÁRIO vii
6 Interpolação 156
6.1 Interpolação polinomial . . . . . . . . . . . . . . . . . . . . . . . . . 157
6.2 Diferenças divididas de Newton . . . . . . . . . . . . . . . . . . . . 161
6.3 Polinômios de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . 163
6.4 Aproximação de funções reais por polinômios interpoladores . . . . 164
6.5 Interpolação linear segmentada . . . . . . . . . . . . . . . . . . . . 167
6.6 Interpolação cúbica segmentada - spline . . . . . . . . . . . . . . . . 168
6.6.1 Spline natural . . . . . . . . . . . . . . . . . . . . . . . . . . 171
6.6.2 Spline fixado . . . . . . . . . . . . . . . . . . . . . . . . . . 173
6.6.3 Spline not-a-knot . . . . . . . . . . . . . . . . . . . . . . . . 174
6.6.4 Spline periódico . . . . . . . . . . . . . . . . . . . . . . . . . 175
Colaboradores 320
Introdução
1
2 Cálculo Numérico
Equações não polinomiais podem ser ainda mais complicadas de resolver exata-
mente, por exemplo:
cos(x) = x e xex = 10
Para resolver o problema de valor inicial
y ′ + xy = x,
y(0) = 2,
2 /2
podemos usar o método de fator integrante e obtemos y = 1 + e−x . Já o cálculo
da solução exata para o problema
y ′ + xy = e−y ,
y(0) = 2,
não é possível.
Da mesma forma, resolvemos a integral
Z 2
2
xe−x dx
1
Representação de números e
aritmética de máquina
3
4 Cálculo Numérico
O sistema de numeração posicional também pode ser usado com outras bases.
Vejamos a seguinte definição.
x = (1001,101)2
= 1 · 23 + 0 · 22 + 0 · 21 + 1 · 20 + 1 · 2−1 + 0 · 2−2 + 1 · 2−3
= 8 + 0 + 0 + 1 + 0,5 + 0 + 0,125 = 9,625.
Verifique no computador!
Exemplo 2.1.4 (Sistema octal). No sistema octal a base é b = 8. Por exemplo:
Verifique no computador!
Exemplo 2.1.5 (Sistema hexadecimal). O sistema de numeração cuja a base é
b = 16 é chamado de sistema hexadecimal. Neste, temos o conjunto de algarismos
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E, F }. Convertendo o número (E2AC)16 para
a base 10 temos
(E2AC)16 = 14 · 163 + 2 · 162 + 10 · 161 + 12 · 160
= 57344 + 512 + 160 + 12 = 58028.
Verifique no computador!
Observação 2.1.2. Python tem algumas sintaxes para representar números em
algumas bases. Por exemplo, temos:
>>> print(0b1001) #bin -> dec
9
>>> print(0157) #oct -> dec
111
>>> print(0xbeba) #hex -> dec
48826
Nos exemplos acima vimos como converter números representados em um sis-
tema de numeração de base b para o sistema decimal. Agora, vamos estudar
como fazer o processo inverso. Isto é, dado um número decimal (X)10 queremos
escrevê-lo em uma outra base b, isto é, queremos obter a seguinte representação:
(X)10 = (dn dn−1 · · · d0 ,d−1 · · · )b
= dn · bn + dn−1 · bn−1 + · · · + d0 · b0 + d−1 · b−1 + d−2 · b−2 + · · ·
X i = dn · bn + · · · + dn−1 bn−1 · · · + d1 · b1 + d0 · b0
e
d−1 d−2
Xf =+ 2 + ···
b1 b
Nosso objetivo é determinar os algarismos {dn , dn−1 , ...}.
Primeiramente, vejamos como tratar a parte inteira X i . Calculando sua divisão
de X i por b, temos:
Xi d0
= + d1 + d2 · b1 + · · · + dn−1 · bn−2 + dn · bn−1 .
b b
9,625 = 9 + 0,625.
9=4·2+1
= (2 · 2 + 0) · 2 + 1
= 23 + 1.
Ou seja, temos que 0,625 = (0,101)2 . Em Python, podemos computar esta con-
versão da parte fracionária da seguinte forma
>>> x=0.625
>>> d = int(2*x); x = 2*x - d; print("d = %d, x = %f" % (d,x))
d = 1, x = 0.250000
>>> d = int(2*x); x = 2*x - d; print("d = %d, x = %f" % (d,x))
d = 0, x = 0.500000
>>> d = int(2*x); x = 2*x - d; print("d = %d, x = %f" % (d,x))
d = 1, x = 0.000000
>>> print(bin(9))
0b1001
>>> print(oct(111))
0157
>>> print(hex(48826))
0xbeba
Exercícios resolvidos
ER 2.1.1. Obtenha a representação do número 125,583̄ na base 6.
Solução. Decompomos 125,583̄ nas suas partes inteira 125 e fracionária 0,583̄.
Então, convertemos cada parte.
Conversão da parte inteira. Vamos escrever o número 125 na base 6. Para
tanto, fazemos sucessivas divisões por 6 como segue:
125 = 20 · 6 + 5 (125 dividido por 6 é igual a 20 e resta 5)
= (3 · 6 + 2) · 6 + 5 = 3 · 62 + 2 · 6 + 5,
logo 125 = (325)6 .
Estes cálculos podem ser feitos em Python com o auxílio do comando % e da
função int. Com o primeiro calculamos o resto da divisão entre dois números,
enquanto que a segunda retorna a parte inteira de um número dado. No nosso
exemplo, temos:
>>> q = 125; d0 = (q % 6); print(q,d0)
>>> q = int(q/6); d1 = (q % 6); print(q,d1)
>>> q = int(q/6); d2 = (q % 6); print(q,d2)
Verifique!
Conversão da parte fracionária. Para converter 0,583̄ para a base 6, faze-
mos sucessivas multiplicações por 6 como segue:
0,583 = 3,5 · 6−1 (0,583 multiplicado por 6 é igual a 3,5)
= 3 · 6−1 + 0,5 · 6−1
= 3 · 6−1 + (3 · 6−1 ) · 6−1
= 3 · 6−1 + 3 · 6−2 ,
logo 0,583 = (0,33)6 . As contas feitas aqui, também podem ser computadas em
Python. Você sabe como? ♦
ER 2.1.2. Obtenha a representação na base 4 do número (101,01)2 .
Solução. Começamos convertendo (101,01)2 para a base decimal:
(1001,101)2 = 1 · 22 + 1 · 20 + 1 · 2−2 = 5,25.
Então, convermetos 5,25 para a base 4. Para sua parte inteira, temos
5 = 1 · 4 + 1 = (11)4 .
Para sua parte fracionária, temos
0,25 = 1 · 4−1 = (0,1)4 .
Logo, (101,01)2 = (11,1)4 . Verifique estas contas no computador! ♦
Exercícios
a) (100)2
b) (100)3
c) (100)b
d) (12)5
e) (AA)16
f) (7,1)8
g) (3,12)5
a) (25,13)8
b) (101,1)2
c) (12F,4)16
d) (11,2)3
a) 7,6 na base b = 5
b) 29,16 na base b = 6
x = (−1)s (M )b × bE ,
onde (M )b = (d1 ,d−1 d−2 d−3 · · · )b , com d1 6= 034 , s é 0 para positivo e 1 para
negativo, E é o expoente.
3
Em algumas referências é usado Mb = (0,d−1 d−2 d−3 ·)b .
4
No caso de x = 0, Mb = (0,00 · · · )b .
2.2.1 Exercícios
>>> int(-0.675*1e2)/1e2
• Por arredondamento:
e, em notação normalizada:
>>> 2**2 == 4
True
>>> np.sqrt(3)**2 == 3
False
s dn−2 ··· d1 d0
Observação 2.4.1. Note que todo registro começando com 1 será um número
negativo.
[11111111] = −27 + 26 + · · · + 2 + 1 = −1
..
.
[10000001] = −27 + 1 = −127
[10000000] = −27 = −128
[01111111] = 26 + · · · + 2 + 1 = 127
..
.
[00000010] = 2
[00000001] = 1
[00000000] = 0
• (−1)d31 (d30 d29 · · · d17 d16 , d15 d14 · · · d1 d0 )2 se o sinal for representado por um
dígito. Observe que nesse caso o zero possui duas representações possíveis:
[10000000000000000000000000000000]
e
[00000000000000000000000000000000]
• (d30 d29 · · · d17 d16 )2 − d31 (215 − 2−16 ) + (0,d15 d14 · · · d1 d0 )2 se o sinal do nú-
mero estiver representado por uma implementação em complemento de um.
Observe que o zero também possui duas representações possíveis:
[11111111111111111111111111111111]
e
[00000000000000000000000000000000]
• (d30 d29 · · · d17 d16 )2 − d31 215 + (0,d15 d14 · · · d1 d0 )2 se o sinal do número estiver
representado por uma implementação em complemento de dois. Nesse caso
o zero é unicamente representado por
[00000000000000000000000000000000]
Observe que 16 dígitos são usados para representar a parte fracionária, 15 são para
representar a parte inteira e um dígito, o d31 , está relacionado ao sinal do número.
representa o número
O expoente deslocado
Uma maneira de representar os expoentes inteiros é deslocar todos eles uma
mesma quantidade. Desta forma permitimos a representação de números negativos
e a ordem deles continua crescente. O expoente é representado por um inteiro sem
sinal do qual é deslocado o BIAS.
Tendo |E| dígitos para representar o expoente, geralmente o BIAS é predefi-
nido de tal forma a dividir a tabela ao meio de tal forma que o expoente um seja
representado pelo sequência [100 . . . 000].
Exemplo 2.4.6. Com 64 bits, pelo padrão IEEE754, temos que |E| := 11. As-
sim, (100 0000 0000)2 = 210 = 1024. Como queremos que esta sequência represente
o 1, definimos BIAS := 1023, pois
1024 − BIAS = 1.
Com 32 bits, temos |E| := 8 e BIAS := 127. E com 128 bits, temos |E| := 15
e BIAS := 16383.
Com |E| = 11 temos
[111 1111 1111] = reservado
[111 1111 1110] = 2046 − BIAS = 102310 = EM AX
..
.=
[100 0000 0001] = 210 + 1 − BIAS = 210
[100 0000 0000] = 210 − BIAS = 110
[011 1111 1111] = 1023 − BIAS = 010
[011 1111 1110] = 1022 − BIAS = −110
..
.=
[000 0000 0001] = 1 − BIAS = −1022 = EM IN
[000 0000 0000] = reservado
Casos especiais
O zero é um caso especial representado pelo registro
[0|000 0000 0000|0000 0000 0000...0000 0000]
Os expoentes reservados são usados para casos especiais:
• c = [0000...0000] é usado para representar o zero (se m = 0) e os números
subnormais (se m 6= 0).
• c = [1111...1111] é usado para representar o infinito (se m = 0) e NaN (se
m 6= 0).
Os números subnormais6 tem a forma
x = (−1)s (0.m1 m2 · · · m51 m52 )2 × 21−BIAS .
6
Note que poderíamos definir números um pouco menores que o M IN R.
Exemplo 2.4.7. Com 64 bits, temos que o épsilon será dado por
1 → (1.0000 0000....0000)2 × 20
ǫ → +(0.0000 0000....0001)2 × 20 = 2−52
(1.0000 0000....0001)2 × 20 6= 1
Assim, ǫ = 2−52 .
Exercícios
|x − x|.
|x − x| |10 − 9,999|
= ≈ 0,0000999... < 5 × 10−4 .
|x| 10
|1 − 1,5|
= 5 × 10−1 < 5 × 100 .
|1|
Exercícios
a) x = π = 3,14159265358979 . . . e x̄ = 3,141
b) x = 1,00001 e x̄ = 1
c) x = 100001 e x̄ = 100000
a) 1,7888544
b) 1788,8544
c) 0,0017888544
d) 0,004596632
e) 2,1754999 × 10−10
f) 2,1754999 × 1010
a) 3276.
b) 42,55.
c) 0,00003331.
∆ = b2 − 4 · a · c
= 0,300000 × 103 × 0,300000 × 103
+ 0,400000 × 101 × 0,100000 × 101 × 0,140000 × 10−1
= 0,900000 × 105 + 0,560000 × 10−1
= 0,900001 × 105
e as raízes:
√
−0,300000 × 103 ± ∆
x1 ,x2 =
0,200000 × 101
√
−0,300000 × 103 ± 0,900001 × 105
=
0,200000 × 101
−0,300000 × 103 ± 0,300000 × 103
=
0,200000 × 101
e
−0,300000 × 103 + 0,300000 × 103
x̃2 = = 0,000000 × 100
0,200000 × 101
No entanto, os valores das raízes com seis dígitos significativos livres de erros de
arredondamento, são:
Observe que a primeira raiz apresenta seis dígitos significativos corretos, mas a
segunda não possui nenhum dígito significativo correto.
√ Observe que isto acontece porque b é muito maior que 4ac, ou seja, b ≈
2
−b − b + 4ac b c
x̃1 = 2b
=− +
2a a b
0,300000 × 103 0,140000 × 10−1
= − −
0,100000 × 101 0,300000 × 103
= −0,300000 × 10 − 0,466667 × 10−4
3
= −0,300000 × 103
−b + b − 4ac
x̃2 = 2b
2a
4ac
= −
4ab
c −0,140000 × 10−1
= − =− = 0,466667 × 10−4
b 0,300000 × 103
Ou seja,
√ os erros na entrada√serão diminuídos pela metade. De fato, usando
y = 2 = 1,4142136... e y ∗ = 1,999 = 1,41386..., obtemos
√ √
∆y 2 − 1,999
= √ ≈ 0,000250031...
y 2
Exemplo 2.8.2. Considere a função f (x) = 10
1−x2
e x∗ = 0,9995 com um erro
absoluto na entrada de 0,0001.
Calculando y ∗ = f (x∗ ) temos
10
y∗ = ≈ 10002,500625157739705173
1 − (0,9995)2
Mas qual é a estimativa de erro nessa resposta? Quantos dígitos significativos
temos nessa resposta?
Sabendo que f ′ (x) = −20x/(1 − x2 )2 , o número de condicionamento é
xf ′ (x) 2x2
κf (x) :=
=
f (x) 1 − x2
Exemplo 2.8.3. Seja f (x) = x exp(x). Calcule o erro absoluto ao calcular f (x)
sabendo que x = 2 ± 0,05.
Solução. Temos que x ≈ 2 com erro absoluto de δx = 0,05. Neste caso, calculamos
δf , isto é, o erro absoluto ao calcular f (x), por:
δf = |f ′ (x)|δx .
δf = |(1 + x)ex | · δx
= |3e2 | · 0,05 = 1,1084.
δy
= 0,03 ⇒ δy = 2 · 0,03 = 0,06
|y|
Aplicando a expressão para estimar o erro em f temos
δf = ∂f
∂x x
∂f
δ + ∂y δy
2e4
= 27
· 0,3 + 2 9+1
9
e4 · 0,06 = 8,493045557
∂A ab
= √ 2 ,
∂a 2 a − b2
√
∂A a2 − b2 b2
= − √ 2 ,
∂b 2 2 a − b2
e substituindo na estimativa para o erro δA em termos de δa = 0,01 e δb = 0,01:
∂A ∂A
+
δA ≈ δ δ
∂a a ∂b b
√ √
3 5 5
≈ · 0,01 + · 0,01 = 0,01565247584
5 10
0,01
Em termos do erro relativo temos erro na hipotenusa de 3
≈ 0,333%, erro no
cateto de 0,01
2
= 0,5% e erro na área de
0,01565247584
√
2 32 −22
= 0,7%
2
Exercícios
lim |x̃n | = ∞
x→∞
n
Figura 2.2: Gráfico de 1 + n1 em função de n em escala linear-logarítmica vari-
ando de 100 até 1018 . Veja o exemplo 2.9.3.
n n
n 1+ 1
n
n 1+ 1
n
Para compreendermos melhor por que existe um limiar N que, quando atin-
gido torna a expressão do exemplo acima identicamente igual a 1, observamos a
sequência de operações realizadas pelo computador:
n → 1 + n → (1 + n) − 1. (2.3)
n
Finalmente, notamos que quando tentamos calcular 1 + 1
n
para n grande, existe
perda de significância no cálculo de 1 + 1/n.
Exemplo 2.9.4 (Analogia da balança). Observe a seguinte comparação interes-
sante que pode ser feita para ilustrar os sistemas de numeração com ponto fixo e
flutuante: o sistema de ponto fixo é como uma balança cujas marcas estão igual-
mente espaçadas; o sistema de ponto flutuante é como uma balança cuja distância
entre as marcas é proporcional à massa medida. Assim, podemos ter uma ba-
lança de ponto fixo cujas marcas estão sempre distanciadas de 100g (100g, 200g,
300g, ..., 1Kg, 1,1Kg,...) e outra balança de ponto flutuante cujas marcas estão
distanciadas sempre de aproximadamente um décimo do valor lido (100g, 110g,
121g, 133g, ..., 1Kg, 1,1Kg, 1,21Kg, ...) A balança de ponto fixo apresenta uma
resolução baixa para pequenas medidas, porém uma resolução alta para grandes
medidas. A balança de ponto flutuante distribui a resolução de forma proporcional
ao longo da escala.
Seguindo nesta analogia, o fenômeno de perda de significância pode ser inter-
pretado como a seguir: imagine que você deseje obter o peso de um gato (apro-
ximadamente 4Kg). Dois processos estão disponíveis: colocar o gato diretamente
na balança ou medir seu peso com o gato e, depois, sem o gato. Na balança
de ponto flutuante, a incerteza associada na medida do peso do gato (sozinho) é
aproximadamente 10% de 4Kg, isto é, 400g. Já a incerteza associada à medida da
uma pessoa (aproximadamente 70Kg) com o gato é de 10% do peso total, isto é,
aproximadamente 7Kg. Esta incerteza é da mesma ordem de grandeza da medida
a ser realizada, tornado o processo impossível de ser realizado, já que teríamos
uma incerteza da ordem de 14Kg (devido à dupla medição) sobre uma grandeza
de 4Kg.
Exercícios resolvidos
ER 2.9.1. Deseja-se medir a concentração de dois diferentes oxidantes no ar. Três
sensores eletroquímicos estão disponíveis para a medida e apresentam a seguintes
respostas:
De forma que
−1
[A] S11 S12 v1 1 S22 −S12 v1
= =
[B] S21 S22 v2 S11 S22 − S12 S21 −S21 S11 v2
Portanto
S22 v1 − S12 v2
[A] =
S11 S22 − S12 S21
−S21 v1 + S11 v2
[B] =
S11 S22 − S12 S21
1 ∂[A] S22
= −
[A] ∂S11 S11 S22 − S12 S21
1 ∂[A] v2 S21 [A] S22
= − + =− ·
[A] ∂S12 S22 v1 − S12 v2 S11 S22 − S12 S21 [B] S11 S22 − S12 S21
1 ∂[A] S12
=
[A] ∂S21 S11 S22 − S12 S21
1 ∂[A] v1 S11 [A] S12
= − = ·
[A] ∂S22 S22 v1 − S12 v2 S11 S22 − S12 S21 [B] S11 S22 − S12 S21
e
1 ∂ [B] v2 S22 [B] S21
= − =
[B] ∂S11 −S21 v1 + S11 v2 S11 S22 − S12 S21 [A] S11 S22 − S12 S21
1 ∂ [B] S21
=
[B] ∂S12 S11 S22 − S12 S21
1 ∂ [B] v1 S21 [B] S11
= − + =−
[B] ∂S21 −S21 v1 + S11 v2 S11 S22 − S12 S21 [A] S11 S22 − S12 S21
1 ∂ [B] S11
= −
[B] ∂S22 S11 S22 − S12 S21
1 1
δ[A] = [20 × 270 × 2% + 20 × 30 × 2% + 30 × 140 × 2% + 30 × 20 × 2%]
[A] 1200
216
= = 0.18 = 18%
1200
1 1
δ[B] = [140 × 270 × 2% + 140 × 30 × 2% + 270 × 140 × 2% + 270 × 20 × 2%]
[B] 1200
426
= = 0.355 = 35.5%
1200
Caso do par 1-3:
270 30
det S = = 53550
15 200
1 1
δ[A] = [200 × 270 × 2% + 200 × 30 × 2% + 30 × 15 × 10% + 30 × 200 × 10%]
[A] 53550
1804,6
= ≈ 0.0337 = 3.37%
52550
1 1
δ[B] = [15 × 270 × 2% + 15 × 30 × 2% + 270 × 15 × 10% + 270 × 200 × 10%]
[B] 53550
5895
= ≈ 0.11 = 11%
53550
Conclusão, apesar de o sensor 3 apresentar uma incerteza cinco vezes maior na
sensibilidade, a escolha do sensor 3 para fazer par ao sensor 1 parece mais ade-
quada. ♦
Exercícios
exp(1/µ)
1 + exp(1/µ)
e
1
exp(−1/µ) + 1
com µ > 0. Verifique que elas são idênticas como funções reais. Teste no compu-
tador cada uma delas para µ = 0,1, µ = 0,01 e µ = 0,001. Qual dessas expressões
é mais adequada quando µ é um número pequeno? Por quê?
de modo que seja possível calcular seus valores para x = 100 utilizando a aritmética
de ponto flutuante ("Double") no computador.
1
Sir Isaac Newton, 1642 - 1727, matemático e físico inglês.
43
44 Cálculo Numérico
f(a)
x∗ b
a x
f(b)
e f (0) = −2 < 0, temos do teorema de Bolzano que existe pelo menos um zero de
f (x) no intervalo (−2, 0). E, portanto, existe pelo menos uma solução da equação
dada no intervalo (−2, 0).
Podemos usar Python para estudarmos esta função. Por exemplo, podemos
definir a função f (x) e computá-la nos extremos do intervalo dado com os seguintes
comandos:
>>> def f(x): return np.exp(x)-x-2
...
>>> f(-2),f(0)
(0.13533528323661281, -1.0)
Alternativamente (e com maior precisão), podemos verificar diretamente o sinal
da função nos pontos desejados com a função numpy.sign:
>>> np.sign(f(-2)*f(0))
-1.0
♦
Quando procuramos aproximações para zeros de funções, é aconselhável isolar
cada raiz em um intervalo. Desta forma, gostaríamos de poder garantir a existência
e a unicidade da raiz dentro de um dado intervalo. A seguinte proposição nos
fornece condições suficientes para tanto.
Em outras palavras, para garantirmos que exista um único zero de uma dada
função diferenciável em um intervalo, é suficiente que ela troque de sinal e seja
monótona neste intervalo.
Exemplo 3.1.2. No exemplo 3.1.1, mostramos que existe pelo menos um zero de
f (x) = ex − x − 2 no intervalo (−2,0), pois f (x) é contínua e f (−2) · f (0) < 0.
Agora, observamos que, além disso, f ′ (x) = ex − 1 e, portanto, f ′ (x) < 0 para
todo x ∈ (−2,0). Logo, da proposição 3.1.1, temos garantida a existência de um
único zero no intervalo dado.
Podemos inspecionar o comportamento da função f (x) = ex − x − 2 e de sua
derivada fazendo seus gráficos no Python. Para tanto, podemos usar o seguinte
código Python:
>>> plt.plot(xx,f(xx))
>>> plt.grid(); plt.show()
Exercícios
E 3.1.2. Mostre que cos x = x tem uma única solução no intervalo [0, π/2].
1
ln(x) + x3 − = 10
x
possui uma única solução positiva.
1
ln(x) + x − =v
x
possui uma solução para cada v real e que esta solução é única.
f (b)
a x(0)
x(1) b x
f (x(0) )
f (a)
a(n) + b(n)
a(n) = a, b(n) = b e x(n) = .
2
Verificamos o critério de parada, isto é, se f (x(n) ) = 0 ou:
|b(n) − a(n) |
< T OL,
2
Licença CC-BY-SA-3.0. Contato: [email protected]
48 Cálculo Numérico
Exemplo 3.2.1. Use o método da bisseção para calcular uma solução de ex = x+2
no intervalo [−2, 0] com precisão T OL = 10−1 .
b−a
|x(n) − x∗ | < , ∀n ≥ 0,
2n+1
Observamos que a hipótese de que f (x) tenha um único zero no intervalo não
é necessária. Se a função tiver mais de um zero no intervalo inicial, as iterações
irão convergir para um dos zeros. Veja o exercício 3.2.3.
Observação 3.2.1. O teorema 3.2.1 nos fornece uma estimativa para a conver-
gência do método da bisseção. Aproximadamente, temos:
1
|x(n+1) − x∗ | . |x(n) − x∗ |.
2
Isto nos leva a concluir que o método da bisseção tem taxa de convergência
linear.
• f - função objetivo
A variável de saída é:
Exercícios
√
E 3.2.1. Considere a equação x = cos(x). Use o método da bisseção com
intervalo inicial [a, b] = [0, 1] e x(1) = (a + b)/2 para calcular a aproximação x(4)
da solução desta equação.
a) Se o método da bisseção for usando com o intervalo inicial [1/2, 3], para qual
raiz as iterações convergem?
√ √
E 3.2.4. O polinômio f (x) = x4 − 4x2 + 4 possui raízes duplas em 2 e − 2.
O método da bisseção pode ser aplicados a f ? Explique.
Use esta estimativa para iniciar o método de bisseção e obtenha o valor da raiz
com pelo menos 6 algarismos significativos para v = 1, 2, 3, 4 e 5.
τ = kθ,
xex = t,
onde t é um número real positivo. Mostre que esta equação possui uma única
solução x∗ que pertence ao intervalo [0, t]. Usando esta estimativa como intervalo
inicial, quantos passos são necessário para obter o valor numérico de x∗ com erro
absoluto inferior a 10−6 quando t = 1, t = 10 e t = 100 através do método da
bisseção? Obtenha esses valores.
que a relação entre a corrente (Id ) e a tensão (vd ) no diodo é dada pela seguinte
expressão:
vd
Id = IR exp −1 ,
vt
onde IR é a corrente de condução reversa e vt , a tensão térmica dada por vt = kTq
com k, a constante de Boltzmann, T a temperatura de operação e q, a carga do
elétron. Aqui IR = 1pA = 10−12 A, T = 300 K. Escreva o problema como uma
equação na incógnita vd e, usando o método da bisseção, resolva este problema
com 3 algarismos significativos para os seguintes casos:
a) V = 30 V e R = 1 kΩ.
b) V = 3 V e R = 1 kΩ.
c) V = 3 V e R = 10 kΩ.
d) V = 300 mV e R = 1 kΩ.
e) V = −300 mV e R = 1 kΩ.
f) V = −30 V e R = 1 kΩ.
g) V = −30 V e R = 10 kΩ.
Dica: V = RId + vd .
y
y=x
x∗
y = g(x)
x∗ x
ex = x + 2 ⇔ ex − x − 2 = 0 ⇔ ex − 2 = x
x(n+1) = g(x(n) ), n ≥ 1,
√ √
• r= A =⇒ A
r
= A;
√ √ √
• r < A =⇒ Ar > A =⇒ A ∈ r, Ar .
√
Ou seja, A sempre está no intervalo entre r e Ar , no qual podemos buscar uma
nova aproximação como, por exemplo, pelo ponto médio:
r+ A
x= r
.
2
Aplicando esse método repetidas vezes, podemos construir a iteração (de ponto
fixo):
x(1) = r
x(n) A
x(n+1) = + (n) , n = 1,2,3,...
2 2x
√
Por exemplo, para obter uma aproximação para 5, podemos iniciar com a
aproximação inicial r = 2 e A = 5. Então, tomamos x(1) = 2 e daí seguem as
aproximações:
2 2,5
x(2) = + = 2,25
2 2
2,25 2,5
x(3) = + = 2,2361111
2 2,25
2,2361111 2,5
x(4) = + = 2,236068
2 2,2361111
2,236068 2,5
x(5) = + = 2,236068
2 2,236068
O método babilônico sugere que a iteração do ponto fixo pode ser uma abor-
dagem eficiente para a solução de equações. Ficam, entretanto, as seguintes per-
guntas:
1. Será que a iteração do ponto fixo é convergente?
2. Caso seja convergente, será que o limite x∗ = limn→∞ x(n) é um ponto fixo?
3. Caso seja convergente, qual é a taxa de convergência?
A segunda pergunta é a mais fácil de ser respondida. No caso de g(x) ser
contínua, se x(n) → x∗ ∈ Dom (g), então:
x∗ = lim x(n) = lim g(x(n−1) ) = g lim x(n−1) = g(x∗ ).
n→∞ n→∞ n→∞
Definição 3.3.1. Uma contração é uma função real g : [a, b] → [a, b] tal que:
• Se |g ′ (x)| < k, 0 < k < 1, para todo x ∈ [a, b], então g(x) é uma contração.
x(n+1) = g(x(n) )
f (a) = a − g(a) ≤ a − a = 0
e
f (b) = b − g(b) ≥ b − b = 0
Se f (a) = a ou f (b) = b, então o ponto fixo existe. Caso contrário, as desigualdades
são estritas e a f (x) muda de sinal no intervalo. Como esta função é contínua, pelo
teorema de Bolzano 3.1.1, existe um ponto x∗ no intervalo (a, b) tal que f (x∗ ) = 0,
ou seja, g(x∗ ) = x∗ . Isto mostra a existência.
Para provar que o ponto fixo é único, observamos que se x∗ e x∗∗ são pontos
fixos, eles devem ser iguais, pois:
A desigualdade |x∗ − x∗∗ | ≤ β|x∗ − x∗∗ | com 0 ≤ β < 1 implica |x∗ − x∗∗ | = 0.
Para demonstrar a convergência da sequência, observamos que:
Daí, temos:
lim |x(n) − x∗ | = 0,
n→∞
|x(n+1) − x∗ | ≤ β|x(n) − x∗ |, n ≥ 1.
−0,85 < − sen (1) ≤ − sen (x) ≤ − sen (1/2) < −0,47.
x(1) = 0,7
x(n+1) = cos(x(n) ), n ≥ 1.
n x(n) ǫn := |x(n) − x∗ |
1 0,70000 3,9E−02
2 0,76484 2,6E−02
3 0,72149 1,8E−02
4 0,75082 1,2E−02
5 0,73113 8,0E−03
6 0,74442 5,3E−03
7 0,73548 3,6E−03
#est. da solucao
xe = sci.optimize.fixed_point(g, 0.7)
#aprox. inicial
x0 = 0.7
eps = np.fabs(x0-xe)
print("%1.5f %1.1e\n" % (x0, eps))
for i in np.arange(7):
x = g(x0);
eps = np.fabs(x-xe);
print("%1.5f %1.1e\n" % (x, eps))
x0 = x
♦
Exemplo 3.3.5. No exemplo 3.3.3, observamos que a função g1 (x) nos forneceu
uma iteração divergente, enquanto que a função g2 (x) forneceu uma iteração con-
vergente (veja a figura 3.5. A razão destes comportamentos é explicada pelo teste
da convergência. Com efeito, sabemos que o ponto fixo destas funções está no
intervalo [1,6, 1,8] e temos:
2.0 1.80
y = g2 (x) y=x
1.9 y = g1 (x)
1.75 x(2)
(1)
1.8 x x(1) x∗
x∗
1.7 1.70
1.6
y=x 1.65
1.5 x (2)
1.4 1.60
1.60 1.65 1.70 1.75 1.80 1.60 1.65 1.70 1.75 1.80
Figura 3.5: Ilustração das iterações do ponto fixo para: (esquerda) y = g1 (x) e
(direita) y = g2 (x). Veja exemplo 3.3.5.
enquanto:
|g2′ (x)| = |1 − 0,05(x + 1)ex | < 0,962, ∀x ∈ [1,6, 1,8].
• Se |g ′ (x∗ )| > 1, então, a distância de x(n) até o ponto fixo x∗ está aumentando
a cada passo.
com base nos valores calculados x(n) . Uma abordagem frequente é analisar a evo-
lução da diferença entre dois elementos da sequência:
∆n = x(n+1) − x(n)
A pergunta natural é: Será que o erro ǫn = x(n) − x∗ é pequeno quando
∆n = x(n+1) − x(n) for pequeno?
x∗ = lim x(n)
n→∞
portanto:
x∗ − x(N ) = x(N +1) − x(N ) + x(N +2) − x(N +1) + x(N +3) − x(N +2) + . . .
∞
X
= x(N +k+1) − x(N +k)
k=0
Portanto:
k
x(N +k+1) − x(N +k) ≈ (x(N +1) − x(N ) ) (g ′ (x∗ ))
E temos:
∞
X
x∗ − x(N ) = x(N +k+1) − x(N +k)
k=0
X∞
k
≈ (x(N +1) − x(N ) ) (g ′ (x∗ ))
k=0
1
= (x(N +1) − x(N ) ) , |g ′ (x∗ )| < 1
1 − g ′ (x∗ )
Observação 3.3.4. Tendo em mente a relação x(n+1) −x(n) ≈ (x(n) −x(n−1) )g ′ (x∗ ),
concluímos:
• Quando 0 < g ′ (x∗ ) < 1, o esquema é monótono e 1−g1′ (x∗ ) > 1, pelo que o
erro ǫN é maior que a diferença ∆N . A relação será tão mais importante
quando mais próximo da unidade for g ′ (x∗ ), ou seja, quando mais lenta for
a convergência. Para estimar o erro em função da diferença ∆N , observamos
(n+1) (n)
(n−1) e
que g ′ (x∗ ) ≈ xx(n) −x−x
∆n
|g ′ (x∗ )| ≈
∆n−1
e portanto
∆N
ǫN ≈ .
1 − ∆∆n−1
n
Exercícios
xex −10
b) g(x) = x − 15
xex −10
c) g(x) = x − 10+ex
a) g(x) = cos(x)
a) ex = x + 2 no intervalo (−2,0).
Use o teorema do ponto fixo para verificar que cada um desses processos con-
verge para a solução da equação x∗ de cos(x) = x. Observe o comportamento
numérico dessas sequências. Qual estabiliza mais rápido com cinco casas deci-
mais? Discuta.
Dica: Verifique que cos([0.5,1]) ⊆ [0.5,1] e depois a mesma identidade para a
função f (x) = 0,4x + 0,6 cos(x).
com erro absoluto inferior a 10−3 usando um método iterativo. Estime o erro
associado ao valor de v = 180 − 100x = 0.052 senh −1 (1013 x), usando cada uma
dessas expressões. Discuta sucintamente o resultado obtido. Dica: Este caso é
semelhante ao problema 3.2.8.
xn+1 = xn − β (xn − x∗ )
xn − x∗ = (1 − β)n−1 (x1 − x∗ ).
x(n+1) = xn + q n ,
x(0) = 0,
onde q = 1 − 10−6 .
a) Calcule o limite
lim x(n)
x∞ = n→∞
analiticamente.
c) Qual deve ser a tolerância especificada para obter o resultado com erro rela-
tivo inferior a 10−2 ?
com x(0) = 10−2 . Prove que {x(n) } é sequência de número reais positivos conver-
gindo para zero. Verifique que são necessários mais de mil passos para que x(n) se
torne menor que 0.9x(0) .
a) Use o teorema do ponto fixo para mostrar que a função g(x) = 1 − sen (x)
possui um único ponto fixo estável o intervalo [ 101
,1]. Construa um método
iterativo x(n+1) = g(x(n) ) para encontrar esse ponto fixo. Use o computador
para encontrar o valor numérico do ponto fixo.
b) Considere a função
g(x) = 10 exp(−x)
e função composta ψ(x) = g ◦ g = g (g(x)). Mostre que ψ possui dois pontos
fixos que não são pontos fixos de g.
Em vista do problema anterior, qual valor de γ você escolheria para que a sequência
x(n) convirja rapidamente para x∗ .
e
h i
V −RI (n)
I (n+1) = IR exp vt
− 1 ,n > 0
B
I (0) = 0
onde α(x) é uma função arbitrária que queremos escolher de forma que a iteração
do ponto fixo tenha ótima taxa de convergência.
Do Teorema do ponto fixo temos que a taxa de convergência é dada em
função do valor absoluto da derivada de g(x). Calculando a derivada temos:
No ponto x = x∗ , temos:
Sabemos que o processo iterativo converge tão mais rápido quanto menor for
|g ′ (x)| nas vizinhanças de x∗ . Isto nos leva a escolher:
g ′ (x∗ ) = 0,
e, então, temos:
1
α(x∗ ) = − ,
f ′ (x∗ )
se f ′ (x∗ ) 6= 0.
A discussão acima nos motiva a introduzir o método de Newton, cujas iterações
são dada por:
(n)
f x
x(n+1) = x(n) − ′ n , n ≥ 1,
f (x )
sendo x(1) uma aproximação inicial dada.
f (x(1) )
x(2) = x(1) − .
f ′ (x(1) )
f(x(1) )
f(x(2) )
x∗
x(3) x(2) x(1) x
Assim, a interseção desta reta com o eixo das abscissas ocorre quando (y = 0):
f (x(1) )
f ′ (x(1) )(x − x(1) ) + f (x(1) ) = 0 ⇒ x = x(1) − .
f ′ (x(1) )
g(x∗ ) = x∗
f ′ (x∗ )f ′ (x∗ ) − f (x∗ )f ′′ (x∗ )
g ′ (x∗ ) = 1 − =0
(f ′ (x∗ ))2
Portanto:
g ′′ (x∗ )
g(x) = x∗ + (x − x∗ )2 + O (x − x∗ )3
2
Com isso, temos:
g ′′ (x∗ ) (n)
x(n+1) = g(x(n) ) = x∗ + (x − x∗ )2 + O (x − x∗ )3 ,
2
ou seja: 2
(n+1)
− x∗ ≤ C x(n) − x∗ ,
x
com constante C = |g ′′ (x∗ )/2|. Isto mostra que o método de Newton tem taxa de
convergência quadrática. Mais precisamente, temos o seguinte teorema.
Teorema 3.4.1 (Método de Newton). Sejam f ∈ C 2 ([a, b]) com x∗ ∈ (a, b) tal
que f (x∗ ) = 0 e:
fornece uma sequência x(n) que converge para x∗ , isto é, x(n) → x∗ quando n → ∞.
Além disso, temos a seguinte estimativa de erro a priori:
2m (2n−1 )
|x(n) − x∗ | ≤ q , n ≥ 2,
M
e a seguinte estimativa de erro a posteriori:
M (n)
|x(n) − x∗ | ≤ |x − x(n−1) |2 , n ≥ 2.
2m
Demonstração. Para n ∈ N, n ≥ 2, temos:
f (x(n) ) 1 h i
xn+1 − x∗ = x(n) − − x ∗
= − f (x (n)
) + (x ∗
− x (n)
)f ′ (n)
(x . (3.1)
f ′ (x(n) ) f (x(n) )
Agora, para estimar o lado direito desta equação, usamos o polinômio de Taylor
de grau 1 da função f (x) em torno de x = x(n) , isto é:
Z x∗
f (x∗ ) = f (x(n) ) + (x∗ − x(n) )f ′ (x(n) ) + f ′′ (t)(x∗ − t) dt.
x(n)
Exemplo 3.4.1. Estime o raio ρ da bacia de atração Kρ (x∗ ) para a função f (x) =
cos(x) − x restrita ao intervalo [0, π/2].
2m
ρ<
M
onde m := min |f ′ (x)| e M := max |f ′′ (x)| com o mínimo e o máximo tomados
em um intervalo [a, b] que contenha o zero da função f (x). Aqui, por exemplo,
podemos tomar [a, b] = [0, π/2]. Como, neste caso, f ′ (x) = − sen (x) − 1, temos
que m = 1. Também, como f ′′ (x) = − cos x, temos M = 1. Assim, concluímos
que ρ < 2 (lembrando que Kρ (x∗ ) ⊂ [0, π/2]). Ou seja, neste caso as iterações
de Newton convergem para o zero de f (x) para qualquer escolha da aproximação
inicial x(1) ∈ [0, π/2]. ♦
Exercícios
tg (x) = 2x2 .
a) Use o método gráfico para isolar as duas primeiras raízes positivas em peque-
nos intervalos. Use a teoria para argumentar quanto à existência e unicidade
das raízes dentro intervalos escolhidos.
b) Calcule cada uma das raízes pelo método de Newton com oito dígitos signi-
ficativos e discuta a convergência comparando com o item b).
trace o gráfico com auxílio do computador e verifique que ela possui uma raiz
positiva. Encontre uma aproximação para esta razão pelo gráfico e use este valor
para inicializar o método de Newton e obtenha uma aproximação para a raiz com
8 dígitos significativos.
a) (1.5) Use o teorema do ponto fixo para provar que se x(0) pertence ao intervalo
[1,3], então a sequência dada iterativamente por
f (x) − f (x0 )
f ′ (x) ≈ , x ≈ x0 .
x − x0
Mais precisamente, o método de Newton é uma iteração de ponto fixo da forma:
onde x(1) é uma aproximação inicial dada e α(x(n) ) = 1/f ′ (x(n) ). Usando a apro-
ximação da derivada acima, com x = x(n) e x0 = x(n−1) , temos:
1 x(n) − x(n−1)
α(x(n) ) = ≈ .
f ′ (x(n) ) f (x(n) ) − f (x(n−1) )
Isto nos motiva a introduzir a iteração do método das secantes dada por:
y
f(x (1) )
f(x (2) )
f(x (3) )
x∗
x (4)x (3) x (2)x (1) x
x(2) − x(1)
x(3) = x(2) − f (x(2) ) .
f (x(2) ) − f (x(1) )
De fato, x(3) é o ponto de interseção da reta secante ao gráfico de f (x) pelos pontos
x(1) e x(2) com o eixo das abscissas. Com efeito, a equação desta reta secante é:
f (x(2) ) − f (x(1) )
y= (x − x(2) ) + f (x(2) ).
x(2) − x(1)
Esta reta intercepta o eixo das abscissas no ponto x tal que y = 0, isto é:
Teorema 3.5.1 (Método das secantes). Seja f ∈ C 2 ([a, b]) uma função com x∗ ∈
(a, b) tal que f (x∗ ) = 0. Sejam, também:
onde {γn }n∈N é a sequência de Fibonacci67 , bem como vale a estimativa a poste-
riori:
M (n)
|x(n) − x∗ | ≤ |x − x(n−1) ||x(n−1) − x(n−2) |, n ≥ 3.
2m
Demonstração. Sejam n ∈ N com n ≥ 2 e x(n) , x(n−1) ∈ Kρ (x∗ ), tal que x(n) 6=
x(n−1) , x(n) 6= x∗ e x(n−1) 6= x∗ . Seja, também:
x(n) − x(n−1)
g(x(n) ,x(n−1) ) := x(n) − f (x(n) ) .
f (x(n) ) − f (x(n−1) )
Com isso, temos:
x(n) − x(n−1)
g(x(n) ,x(n−1) ) − x∗ = x(n) − f (x(n) ) − x∗
f (x ) − f (x
(n) (n−1) )
( )
∗ f (x ) − f (x )
(n) (n−1) (n) (n−1)
x −x (n) (n) ∗
= (x − x ) − f (x ) + f (x ) .
f (x(n) ) − f (x(n−1) ) x(n) − x(n−1)
Logo, temos:
f (x(n) ) − f (x(n−1) ) f (x(n) ) − f (x∗ )
− =
x(n) − x(n−1) x(n) − x∗ (3.5)
Z 1h i
′ (n) (n−1) (n) ′ (n) ∗ (n)
f (x + r(x −x )) − f (x + r(x − x )) dr.
0
6
Leonardo Fibonacci, c. 1170 - c. 1250, matemático italiano.
7
A sequência de Fibonacci {γn }n∈N é definida por γ0 = γ1 = 1 e γn+1 = γn − γn−1 , n ≥ 1.
|∆n |
Iteração Linear < erro
′
ǫn+1 ≈ |φ (x )|εn ∗ 1− ∆∆n
n−1
linear (p = 1) ∆n < ∆n−1
Quadrática 1 f ′′ (x∗ )
Newton ǫn+1 ≈ ′ ∗ ε2n |∆n | < erro
(p = 2) 2 f (x )
√
5+1 f ′′ (x∗ )
p = εn+1 ≈ ′ ∗ εn εn−1
Secante 2 f (x ) |∆n | < erro
≈ 1,618 ≈ M εφn
condição seja satisfeita por alguns poucos passos consecutivos. Outros critérios
podem ser usados. No métodos das secantes, deve-se ter o cuidado de evitar
divisões por zero quando xn+1 − xn muito pequeno em relação à resolução do
sistema de numeração.
Exercícios
E 3.6.1. Refaça as questões 3.4.3, 3.4.4, 3.4.5 e 3.4.6, usando o método das
secantes.
E 3.6.2. Dê uma interpretação geométrica ao método das secantes. Qual a
vantagem do método das secantes sobre o método de Newton?
E 3.6.5. Seja uma função f (x) dada duas vezes continuamente diferenciável.
Faça uma análise assintótica para mostrar que as iterações do método das secantes
satisfazem:
|x(n+1) − x∗ | ≈ C|x(n) − x∗ ||x(n−1) − x∗ |,
para aproximações iniciais x(1) e x(2) suficientemente próximas de x∗ , onde f (x∗ ) =
0.
E 3.7.3. A equação
cos(πx) = e−2x
tem infinitas raízes. Usando métodos numéricos encontre as primeiras raízes dessa
equação. Verifique a j-ésima raiz (zj ) pode ser aproximada por j − 1/2 para j
grande. Use o método de Newton para encontrar uma aproximação melhor para
zj .
a) R = 0Ω
b) R = 10Ω
c) R = 50Ω
d) R = 100Ω
E) R = 500Ω
c(t) = Ate−λt
Ano população
1960 70992343
1970 94508583
1980 121150573
1991 146917459
Use esses parâmetros para calcular a população em 1980 e compare com o valor
do censo. Dica: considere PP (10)−P
(31)−P (0)
(0)
e reduza o sistema a uma equação apenas na
variável λ.
E 3.7.7. (Fluidos) Uma boia esférica flutua na água. Sabendo que a boia tem
10ℓ de volume e 2Kg de massa. Calcule a altura da porção molhada da boia.
E 3.7.8. (Fluidos) Uma boia cilíndrica tem secção transversal circular de raio
10cm e comprimento 2m e pesa 10Kg. Sabendo que a boia flutua sobre água com
o eixo do cilindro na posição horizontal, calcule a altura da parte molhada da boia.
x2
E 3.7.13. Encontre os pontos onde a elipse que satisfaz 3
+ y 2 = 1 intersepta
a parábola y = x2 − 2.
Hidrocarboneto A B C
N-pentano 9.2131 2477.07 -39.94
N-heptano 9.2535 2911.32 -56.51
87
88 Cálculo Numérico
x+y+z = 1
4x + 4y + 2z = 2
2x + y − z = 0
Fixamos, agora, na segunda linha. Dividimos esta linha pelo valor do elemento
em sua diagonal, isto nos fornece
1 1 0 0
0 1 0 −1
0 0 1 1
x+y+z = 1
2x + y − z = 0
2x + 2y + z = 1
#L2 <-> L1
aux = np.copy(E[1,:])
E[1,:] = np.copy(E[0,:])
E[0,:] = np.copy(aux)
print(E)
#zera E[1:2,0]
E[1,:] = E[1,:] - (E[1,0]/E[0,0])*E[0,:]
E[2,:] = E[2,:] - (E[2,0]/E[0,0])*E[0,:]
print(E)
#zera E[2,1]
E[2,:] = E[2,:] - (E[2,1]/E[1,1])*E[1,:]
print(E)
#sub. regressiva
x = np.zeros(3)
x[2] = E[2,3]/E[2,2];
x[1] = (E[1,3] - E[1,2]*x[2])/E[1,1];
x[0] = (E[0,3] - E[0,2]*x[2] - E[0,1]*x[1])/E[0,0]
print(x)
♦
Exemplo 4.1.4 (Problema com elementos com grande diferença de escala). Re-
solva o seguinte sistema usando eliminação gaussiana sem e com pivotamento par-
cial. Discuta, em cada caso, o resultado frente a aritmética de ponto flutuante
quando 0 < |ǫ| ≪ 1.
ε 2 x 4
=
1 ε y 3
Temos
3 − 4/ε
y=
ε − 2/ε
e
4 − 2y
x=
ε
Observe que a expressão obtida para y se aproximada de 2 quando ε é pequeno:
3 − 4/ε 3ε − 4 −4
y= = 2 −→ = 2, quando ε → 0.
ε − 2/ε ε −2 −2
Exercícios Resolvidos
ER 4.1.1. Resolva o seguinte sistema por eliminação gaussiana com pivotamento
parcial.
2y + 2z = 8
x + 2y + z = 9
x+y+z = 6
Portanto x = 2, y = 3 e z = 1.
♦
Exercícios
x+y+z = 0
x + 10z = −48
10y + z = 25
Usando eliminação gaussiana com pivotamento parcial (não use o computador para
resolver essa questão).
x+y+z = 0
x + 10z = −48
10y + z = 25
Usando eliminação gaussiana com pivotamento parcial (não use o computador para
resolver essa questão).
flutuante (flops) feitas pelo algoritmo (observe que o tempo total não depende
apenas disso, mas também de outros fatores como memória, taxas de transferências
de dados da memória para o cpu, redes,...). Entretanto vamos nos concentrar na
contagem do número de operações (flops) para realizar determinada tarefa.
No passado (antes dos anos 80), os computadores demoravam mais tempo
para realizar operações como multiplicação e divisão, se comparados a adição ou
subtração. Assim, em livros clássicos eram contados apenas o custo das operações
× e /. Nos computadores atuais as quatro operações básicas levam o mesmo
tempo. Entretanto, na maioria dos algoritmos de álgebra linear o custo associado
as multiplicações e divisões é proporcional ao custo das somas e subtrações (pois
a maioria dessas operações podem ser escritas como a combinação de produtos
internos). Dessa forma, na maior parte deste material levaremos em conta somente
multiplicações e divisões, a não ser que mencionado o contrário.
Teremos em mente que a ideia é estimar o custo quando lidamos com vetores
e matrizes muito grande, isto é, o custo quando estas dimensões crescem infinita-
mente.
x = [a × x1 , a × x2 , ...,a × xn ]
ax (4.1)
C = n flops. (4.2)
Exemplo 4.2.2 (Produto vetor-vetor). Qual o custo para calcular o produto in-
terno x · y ?
x · y = x1 × y1 + x2 × y2 + ... + xn × yn (4.3)
Note que esse produto tem o custo do produto vetor-vetor, ou seja, 2n − 1. Como
temos n × n elementos em D, o custo total para multiplicar duas matrizes é2
4.4 Fatoração LU
Considere um sistema linear Ax = b, onde a matriz A é densa3 . A fim de
resolver o sistema, podemos fatorar a matriz A como o produto de uma matriz L
triangular inferior e uma matriz U triangular superior, ou seja, A = LU .
Sendo assim, o sistema pode ser reescrito da seguinte forma:
Ax = b
(LU )x = b
L(U x) = b
Ly = b e Ux = y
e fazemos
A2,: ⇐ A2,: − L21 A1,:
Note que denotamos Ai,: para nos referenciarmos a linha i de A. Da mesma
forma, se necessário usaremos A:,j para nos referenciarmos a coluna j de A.
Para zerar o primeiro elemento da terceira linha de A, temos
e fazemos
A3,: ⇐ A3,: − L31 A1,:
até chegarmos ao último elemento da primeira coluna de A.
Repetimos o processo para as próximas colunas, escalonando a matriz A e
coletando os elementos Lij abaixo da diagonal5 .
3
Diferentemente de uma matriz esparsa, uma matriz densa possui a maioria dos elementos
diferentes de zero.
4
Não vamos usar pivotamento nesse primeiro exemplo.
5
Perceba que a partir da segunda coluna para calcular Lij não usamos os elementos de A,
mas os elementos da matriz A em processo de escalonamento
x1 + x2 + x3 = −2
2x1 + x2 − x3 = 1
2x1 − x2 + x3 = 3
y1 = −2
2y1 + y2 = 1
2y1 + 3y2 + y3 = 3
x1 + x2 + x3 = −2
−x2 − 3x3 = 5
8x3 = −8
2n2 flops.
Somando esses 3 custos, temos que o custo para resolver um sistema linear
usando fatoração LU é
2n3 3n2 n
+ − flops.
3 2 6
Quando n cresce, prevalessem os termos de mais alta ordem, ou seja,
2n3 3n2 n 2n3 3n2 2n3
O( + − ) = O( + ) = O( )
3 2 6 3 2 3
Licença CC-BY-SA-3.0. Contato: [email protected]
4.4. FATORAÇÃO LU 103
8n3 n2 n
− − flops.
3 2 6
xi = e i ,
Ax i=1:n
a1 = cn = 0.
xn = d′n (4.21)
xi = d′i − c′i xi+1 , i = n − 1, n − 2, . . . , 1. (4.22)
1 c′1 ′
x1 d1
a2 b2 c2 x2 d2
.
b3 . .
a3 x
3 = d ,
3
... ... .. ..
.
cn−1 .
an bn xn dn
d2 −a2 d′1
onde c′2 = c2
b2 −a2 c′1
e d′2 = b2 −a2 c′1
.
O próximo passo consiste em substituir a terceira linha por ela mesma subs-
a = (0, 1, 1, 1, 1)
b = (2, 2, 2, 2, 2)
c = (1, 1, 1, 1, 0)
d = (4, 4, 0, 0, 2)
x5 = d′5 = 1
4 4
x4 = d′4 − c′4 · x5 = − ·1=0
5 5
3
x3 = d′3 − c′3 · x4 = −1 − · 0 = −1
4
4 2
x2 = d′2 − c′2 · x3 = − · (−1) = 2
3 3
1
x1 = d′1 − c′1 · x2 = 2 − · 2 = 1
2
E assim, obtemos o vetor x = [1, 2, −1, 0, 1].
def TDMA(a,b,c,d):
#preliminares
a = a.astype('double')
b = b.astype('double')
c = c.astype('double')
d = d.astype('double')
#calcula cl e dl
cl[0]=c[0]/b[0]
for i in np.arange(1,n-1,1):
cl[i]=c[i]/(b[i]-a[i]*cl[i-1])
dl[0]=d[0]/b[0]
for i in np.arange(1,n,1):
dl[i]=(d[i]-a[i]*dl[i-1])/(b[i]-a[i]*cl[i-1])
return x
import numpy as np
def TDMA(a,b,c,d):
#preliminares
a = a.astype('double')
b = b.astype('double')
c = c.astype('double')
d = d.astype('double')
#inicializa vetor x
x=np.zeros(n)
#calcula cl e dl sobrescrevendo-os em c e d
c[0]=c[0]/b[0]
for i in np.arange(1,n-1,1):
c[i]=c[i]/(b[i]-a[i]*c[i-1])
d[0]=d[0]/b[0]
for i in np.arange(1,n,1):
d[i]=(d[i]-a[i]*d[i-1])/(b[i]-a[i]*c[i-1])
for i in np.arange(n-2,-1,-1):
x[i]=d[i]-c[i]*x[i+1]
return x
>>>a=np.array([1,1,1,1,1])
>>>b=np.array([2,2,2,2,2])
>>>c=np.array([1,1,1,1,1])
>>>d=np.array([4,4,0,0,2])
>>>TDMA(a,b,c,d)
x1 − x2 = 0
−xj−1 + 5xj − xj+1 = cos(j/10), 2 ≤ j ≤ 10
x11 = x10 /2 (4.25)
Pequenas variações nos coeficientes das matrizes fazem as soluções ficarem bem
distintas, isto é, pequenas variações nos dados de entrada acarretaram em grandes
variações na solução do sistema. Quando isso acontece, dizemos que o problema é
mal-condicionado.
Precisamos uma maneira de medir essas variações. Como os dados de entrada
e os dados de saída são vetores (ou matrizes), precisamos introduzir as definições
de norma de vetores e matrizes.
c) Em construção ...
d) Em construção ...
Solução.
kvk1 = 1 + 2 + 3 + 0 = 6
√ √
kvk2 = 1 + 22 + 32 + 02 = 14
kvk∞ = max{1,2,3,0} = 3
♦
Solução.
>>> A = np.array([[3,-5,7],
... [1,-2,4],
... [-8,1,-7]], dtype='double')
>>> np.linalg.norm(A,1)
18
>>> np.linalg.norm(A,np.inf)
16
>>> np.linalg.norm(A,2)
13.986577820518308
A(x + δx ) = y + δy
e, portanto,
Aδx = δy .
A = np.array([[71,41],[51,30]],dtype='double')
print(np.linalg.cond(A,p=1))
print(np.linalg.cond(A,p=2))
print(np.linalg.cond(A,p=np.inf))
e, análogo para a matriz A2 . ♦
Exercícios
(j/10)3
xj = sen (j/10), yj = j/10 zj = j/10 − , j = 1, . . . ,10
6
Use o Pythonpara construir os seguintes vetores de erro:
|xj − yj | |xj − zj |
ej = fj =
|xj | xj
Calcule as normas 1, 2 e ∞ de e e f
def jacobi(A,b,x0,tol,N):
#preliminares
A = A.astype('double')
b = b.astype('double')
x0 = x0.astype('double')
n=np.shape(A)[0]
x = np.zeros(n)
it = 0
#iteracoes
while (it < N):
it = it+1
#iteracao de Jacobi
for i in np.arange(n):
x[i] = b[i]
for j in np.concatenate((np.arange(0,i),np.arange(i+1,n))):
x[i] -= A[i,j]*x0[j]
x[i] /= A[i,i]
#tolerancia
if (np.linalg.norm(x-x0,np.inf) < tol):
return x
#prepara nova iteracao
x0 = np.copy(x)
raise NameError('num. max. de iteracoes excedido.')
10x + y = 23
x + 8y = 26
23 − y (k)
x(k+1) =
10
26 − x(k+1)
y (k+1) =
8
23 − y (1)
x(2) = = 2,3
10
26 − x(2)
y (2) = = 2,9625
8
23 − y (2)
x(3) = = 2,00375
10
26 − x(3)
y (3) = = 2,9995312
8
def gauss_seidel(A,b,x0,tol,N):
#preliminares
A = A.astype('double')
b = b.astype('double')
x0 = x0.astype('double')
n=np.shape(A)[0]
x = np.copy(x0)
it = 0
#iteracoes
while (it < N):
it = it+1
#iteracao de Jacobi
for i in np.arange(n):
x[i] = b[i]
for j in np.concatenate((np.arange(0,i),np.arange(i+1,n))):
x[i] -= A[i,j]*x[j]
x[i] /= A[i,i]
print(x[i],A[i,i])
#tolerancia
if (np.linalg.norm(x-x0,np.inf) < tol):
return x
#prepara nova iteracao
x0 = np.copy(x)
raise NameError('num. max. de iteracoes excedido.')
>>> A = np.array([[3,1,-1],
... [-1,-4,1],
... [1,-2,5]],
... dtype='double')
>>> D = np.diag(np.diag(A))
>>> L = np.tril(A)-D
>>> U=np.triu(A)-D
Iteração de Jacobi
Vamos, agora, usar a decomposição discutida acima para construir a matriz de
iteração TJ e o vetor de iteração cJ associado ao método de Jacobi. Neste caso,
temos:
Ax = b ⇔ (L + D + U )x = b
⇔ Dx = −(L + U )x + b
−1
⇔ x = −D−1 (L + U ) x + |D{z b} .
| {z }
=:cJ
=:TJ
x(k+1) = TJ x(k) + cJ , k ≥ 1,
Iteração de Gauss-Seidel
A forma matricial da iteração do método de Gauss-Seidel também pode ser
construída com base na decomposição A = L + D + U . Para tando, fazemos:
Ax = b ⇔ (L + D + U )x = b
⇔ (L + D)x = −U x + b
⇔ x = −(L + D)−1 U x + (L + D)−1 b
| {z } | {z }
=:TG =:cG
x(k+1) = TG x(k) + cG , k ≥ 1,
com x(1) uma aproximação inicial dada, sendo TG := −(L + D)−1 U a matriz de
iteração e cJ = (L + D)−1 b o vetor da iteração.
Condições de convergência
Aqui, vamos discutir condições necessárias e suficientes para a convergência de
métodos iterativos. Isto é, dado um sistema Ax = b e uma iteração:
x(k+1) = T x(k) + c, k ≥ 1,
Agora, suponhamos que ρ(T ) < 1 e seja 0 < ǫ < 1 − ρ(T ). Então, existe
1 ≤ p ≤ ∞ tal que (veja [8], Teorema 3, página 12):
kT kp ≤ ρ(T ) + ǫ < 1.
Assim, temos:
lim kT k kp ≤ lim kT km
p = 0.
k→∞ k→∞
∞
X
−1
o que mostra que (I − T ) = T k.
k=0
x(k+1) = T x(k) + c
x∗ − x(k+1) = (T x∗ + c) − (T x(k) + c)
= T (x∗ − x(k) )
..
.
= T (k) (x∗ − x(1) ) = T (k) y.
Condição suficiente
Uma condição suficiente porém não necessária para que os métodos de Gauss-
Seidel e Jacobi convirjam é a que a matriz seja estritamente diagonal domi-
nante.
e para ao menos um i, aii é estritamente maior que a soma dos elementos fora da
diagonal.
Teorema 4.7.3. O erro relativo e o resíduo estão relacionados como (veja [3])
kx − x(k) k krk
≤ κ(A)
kxk kbk
Exercícios
x1 − x2 = 1
−x1 + 2x2 − x3 = 1
−x2 + (2 + ε)x3 − x4 = 1
−x3 + 2x4 − x5 = 1
x4 − x5 = 1
x1 − x2 = 0
−xj−1 + 5xj − xj+1 = cos(j/10), 2 ≤ j ≤ 10
x11 = x10 /2 (4.28)
x1 + 10x2 + 3x3 = 27
4x1 + x3 = 6
2x1 + x2 + 4x3 = 12
Suponha que um dos autovalores seja estritamente maior que os demais, isto é,
k
desde que ξ1 6= 0. Como λλn1 ≤ λn−1
λ1
≤ · · · λλ31 ≤ λλ12 < 1, então λλ1j tende a 0 para
cada j, 2 ≤ j ≤ n. Devido a normalização a cada passo da sequência,
Ak x(1)
x(k+1) =
kAk x(1) k2
converge para ±e1 , e1 = (1, 0, 0, . . . , 0). Também, a sequência
T
λ(k) = x(k) Ax(k)
Considere a mesma hipótese sobre os autovalores: |λ1 | > |λ2 | ≥ |λ3 | ≥ · · · ≥ |λn |.
Então
Ak = (U DU −1 )(U DU −1 )(U DU −1 ) · · · (U DU −1 ) = U Dk U −1
visto que U U −1 = I. Suponha que o chute inicial x(1) pode ser escrito da forma
com ξ1 6= 0. Então
ξ1
k
λ1 ξ1 1
k
k ξ2 λ2
ξ2 λ2 ξ2 ξ λ
1 1
k
Ak x(1) = (U Dk U −1 )U ξ
3 =U λk ξ
3 3 = λk1 ξ1 U ξ3 λ3 .
ξ1 λ1
.. .. ..
. . .
k
ξn λn
ξn λkn ξn ξ1 λ1
Observamos que se x está na forma (4.29), então Ak x pode ser escrito como
!k
n n n
k
X
k
X X βj λj
A x= βj A vj = βj λkj vj = β1 λk1 v 1 + vj .
j=1 j=1 j=2 β1 λ1
Como λλ1j < 1 para todo j ≥ 2, temos
n
!k
X βj λj
vj → 0.
j=2 β1 λ1
Assim,
k
Ak x β1 λk1 λ
2
= v 1 + O . (4.30)
kAk xk2 kAk xk λ1
Ak x
Como a norma de é igual a um, temos
kAk xk
β λk
1 1
→1
v
kAk xk
1
e, portanto,
β λk
1 1 1
→ .
kAk xk kv1 k
β1 λk1
Ou seja, se definimos α(k) = , então
kAk xk
|α(k) | → 1.
Ak x
− α(k) v1 → 0.
kAk xk
Ak x(1)
x(k) = .
kAk x(1) k
Como λλji −σ
−σ
< 1, o último vetor converge para ±ei e
(A − σI)−k x0
xk =
k(A − σI)−k x0 k2
Exercícios
E 4.9.1. O circuito linear da figura 4.9.1 pode ser modelado pelo sistema dado
a seguir. Escreva esse sistema na forma matricial sendo as tensões V1 , V2 , V3 , V4 e
V5 as cinco incógnitas. Resolva esse problema quando V = 127 e
a) R1 = R2 = R3 = R4 = 2 e R5 = R6 = R7 = 100 e R8 = 50
b) R1 = R2 = R3 = R4 = 2 e R5 = 50 e R6 = R7 = R8 = 100
V1 = V
V1 − V2 V3 − V2 V2
+ − = 0
R1 R2 R5
V2 − V3 V4 − V3 V3
+ − = 0
R2 R3 R6
V3 − V4 V5 − V4 V4
+ − = 0
R3 R4 R7
V4 − V5 V5
− = 0
R4 R8
V1 V2 V3 V4 V5
R1 R2 R3 R4
V R5 R6 R7 R8
Caso V1 V2 V3 V4 V5
a
b
Então, refaça este problema reduzindo o sistema para apenas 4 incógnitas (V2 ,
V3 , V4 e V5 ).
a) Encontre o polinômio P (x) = ax2 + bx + c que passa pelos pontos (−1, − 3),
(1, − 1) e (2,9).
Podemos escrever este problema na forma vetorial definindo o vetor x = [x1 ,x2 , . . . ,xn ]T
138
139
e a função vetorial
f1 (x1 ,x2 , . . . ,xn )
f2 (x1 ,x2 , . . . ,xn )
F (x) =
..
.
fn (x1 ,x2 , . . . ,xn )
Portanto,
∂f1 ∂f1 ∂f1
∂x1 ∂x2
··· ∂xn
(0)
x1 − x1
∂f2 ∂f2 ∂f2
∂x1 ∂x2
··· ∂xn
(0)
x2 − x2
F (x) − F (x(0) ) ≈
.. .. ... ..
, (5.1)
. . .
..
.
∂fn ∂fn ∂fn
∂x1 ∂x2
··· ∂xn
xn − x(0)
n
• A aproximação
x(k) é definida como o ponto x em que a linearização F (x(k) )+
JF x(k) x − x(k) é nula, ou seja:
F (x(k) ) + JF x(k) x(k+1) − x(k) = 0
Assim, ∆(k) = −JF−1 x(k) F (x(k) ), ou seja, ∆(k) resolve o problema linear:
JF x(k) ∆(k) = −F (x(k) )
x21
+ x22 = 1
3
x2
x21 + 2 = 1
4
Para tal, definimos a função F (x):
x21
+ x22 − 1
3
F (x) =
2 x22
x1 + −1
4
cuja jacobiana é:
2x1
2x 2
3
JF =
2x1 x22
Faremos a implementação numérica em Python. Para tal definimos as funções
que implementarão F (x) e a JF (x)
>>> def F(x):
... y = np.zeros(2)
... y[0] = x[0]**2/3 + x[1]**2 - 1
... y[1] = x[0]**2 + x[1]**2/4 - 1
... return y
...
>>> def JF(x):
... y = np.zeros((2,2))
... y[0,0] = 2*x[0]/3
... y[0,1] = 2*x[1]
... y[1,0] = 2*x[0]
... y[1,1] = x[1]/2
... return y
...
Ou simplesmente
>>> x = x - np.linalg.inv(JF(x)).dot(F(x))
q q
Observe que as soluções exatas desse sistema são ± 9
11
, ± 8
11
.
x21 = cos(x1 x2 ) + 1
sen (x2 ) = 2 cos(x1 )
Solução. Vamos, aqui, dar as principais ideias para se obter a solução usando o
método de Newton. Começamos definindo nossa aproximação inicial por x(1) =
(1,5, 0,5). Então iteramos:
onde
2
x1 − cos(x1 x2 ) − 1
F (x) =
sen (x2 ) − 2 cos(x1 )
e sua jacobiana é
2x1 + x2 sen (x1 x2 ) x1 sen (x1 x2 )
JF (x) =
2 sen (x1 ) cos(x2 )
return y
def JF(x):
y = np.zeros((2,2))
y[1,0] = 2*np.sin(x[0])
y[1,1] = np.cos(x[1])
return y
>>> x = np.array([1.5,0.5])
>>> x=x-np.linalg.inv(JF(x)).dot(F(x))
def newton(F,JF,x0,TOL,N):
#preliminares
x = np.copy(x0).astype('double')
k=0
#iteracoes
while (k < N):
k += 1
#iteracao Newton
delta = -np.linalg.inv(JF(x)).dot(F(x))
x = x + delta
#criterio de parada
if (np.linalg.norm(delta,np.inf) < TOL):
return x
Exercícios
f (x,y) = x2 y + cos(xy) − 4
2x1 − x2 = cos(x1 )
−x1 + 2x2 − x3 = cos(x2 )
−x2 + x3 = cos(x3 )
a) Faça um esboço das duas curvas e entenda o problema. Verifique que existem
dois pontos de intersecção, um no primeiro quadrante e outro no segundo
quadrante do plano xy.
b) A partir de seu esboço, encontre aproximações para x e y em cada ponto.
x 0
c) Escreva o problema na forma F =
y 0
d) Encontre a jacobiana JF .
e) Construa a iteração do método de Newton.
f) Implemente no computador.
g) Resolva o sistema analiticamente e compare as respostas.
d ) Encontre a jacobiana JF .
e ) Construa a iteração do método de Newton.
f ) Implemente no Python.
g ) Transforme o sistema em um problema de uma única variável e compare
com a resposta do problema 3.4.1.
E 5.1.5. Encontre uma aproximação com erro inferior a 10−5 em cada incóg-
nita para a solução próxima da origem do sistema
6x − 2y + ez = 2
sen (x) − y + z = 0
sen (x) + 2y + 3z = 1
a) k = 1000 Nm/rad
b) k = 500 Nm/rad
c) k = 100 Nm/rad
d) k = 10 Nm/rad
Obs:Você deve escolher valores para iniciar o método. Como você interpretaria
fisicamente a solução para produzir palpites iniciais satisfatórios? O que se altera
entre o caso a e o caso d?
E 5.1.9. (estática - problemas de três variáveis) Considere, agora, o sistema
mecânico semelhante ao do problema 5.1.8, porém constituído de três segmentos
de mesmo comprimento L presos entre si e a uma parede por articulações.
O momento em cada articulação é proporcional à deflexão com constante de
proporcionalidade k. Os segmentos são feitos de material homogêneo de peso P .
A condição de equilíbrio pode ser expressa em termos dos ângulos θ1 , θ2 e θ3
conforme:
5P L
kθ1 = cos θ1 + k (θ2 − θ1 )
2
3P L
k (θ2 − θ1 ) = cos θ2 + k (θ3 − θ2 )
2
PL
k (θ3 − θ2 ) = cos θ3
2
Considere P = 10N, L = 1m e calcule os ângulos θ1 , θ2 e θ3 quando:
a) k = 1000Nm/rad
b) k = 100Nm/rad
c) k = 10Nm/rad
tan−1 (x) + x = y + y 3
Com base no gráfico, encontre soluções aproximadas para o problema e use-as para
iniciar o método de Newton-Raphson. Encontre as raízes com erro inferior a 10−5 .
restrita à condição
E 5.1.15. A função f (x,y,z) = sen (x) + sen (2y) + sen (3z) possui um máximo
quando x = π/2, y = π/4 e z = π/6. Calcule numericamente este ponto.
5.2.1 Gradiente
Considere primeiramente uma função f : Rn → R, ou seja, uma função que
mapeia n variáveis reais em um único real, por exemplo:
d d
g ′ (h) = g(h) = f (x(0) + hv).
dh dh
Pela regra da cadeia temos:
n
d X ∂f dxj
f (x(0) + hv) = .
dh j=1 ∂xj dh
(0)
Observamos que xj = xj + hvj , portanto
dxj
= vj
dh
Licença CC-BY-SA-3.0. Contato: [email protected]
5.2. LINEARIZAÇÃO DE UMA FUNÇÃO DE VÁRIAS VARIÁVEIS 153
Assim:
n
d (0)
X ∂f
f (x + hv) = vj .
dh j=1 ∂xj
Observamos que esta expressão pode ser vista como o produto interno entre o
gradiente de f e o vetor v:
∂f
∂x1 v1
∂f
v2
∂x2
∇f = .
v= .
.. ..
∂f
∂xn
vn
f1 x (0)
+ ∇ f1 (x ) x − x
T (0) (0)
+ O(kx − x k )
(0) 2
f2 x(0) + ∇Tf2 (x(0) ) x − x(0) + O(kx − x(0) k2 )
F (x) =
..
.
fn x(0) + ∇Tfn (x(0) ) x − x(0) + O(kx − x(0) k2 )
| {z }
Vetor coluna
ou, equivalentemente:
∇ f1 (x )
(0) T (0)
f1 x
f2 x (0)
∇Tf2 (x(0) )
F (x) = + x − x(0) +O(kx − x(0) k2 )
|
.. .. {z }
.
. Vetor coluna
fn x(0) ∇Tfn (x(0) )
| {z } | {z }
Vetor coluna Matriz jacobiana
2
F (x) = F x(0) + JF (x(0) ) x − x(0) + O
x − x(0)
fj , ou seja:
∂f1 ∂f1 ∂f1
∂x1 ∂x2
··· ∂xn
∂f2 ∂f2 ∂f2
∂x1 ∂x2
··· ∂xn
∂(f1 ,f2 , . . . ,fn )
JF = =
∂(x1 ,x2 , . . . ,xn )
.. .. .. ..
. . . .
∂fn ∂fn ∂fn
∂x1 ∂x2
··· ∂xn
∂fi
(JF )ij =
∂xj
Interpolação
f (xi ) = yi , i = 1, 2, . . . , n.
a + bx1 = y1 a+b=1
isto é
a + bx2 = y2 a + 2b = 2
156
6.1. INTERPOLAÇÃO POLINOMIAL 157
Figura 6.1: Exemplo de interpolação de dois pontos por uma reta, veja o exem-
plo 6.0.1.
Figura 6.2: Polinômio interpolador do conjunto de pontos {(0, 1), (1, 6), (2, 5),
(3, −8)}. Veja o exemplo 6.1.1.
Teorema 6.1.2. Seja {(xi ,yi )}ni=1 um conjunto de n pares ordenados de números
reais tais que xi 6= xj se i 6= j, então existe um único polinômio p(x) de grau n − 1
ou inferior que passa por todos os pontos dados, isto é, p(xi ) = yi , i = 1, . . . , n.
É fácil ver que se as abscissas são diferentes dois a dois, então o determinante é
não nulo. Disto decorre que a matriz envolvida é inversível e, portanto, o sistema
possui uma solução que é única.
Esta abordagem direta que usamos no exemplo 6.1.1 e na demonstração do
teorema 6.1.2 se mostra ineficiente quando o número de pontos é grande e quando
existe grande variação nas abscissas. Neste caso, a matriz de Vandermonde é mal
condicionada (ver [6]), o que acarreta um aumento dos erros de arredondamento
na solução do sistema.
Uma maneira de resolver este problema é escrever o polinômio em uma base
que produza um sistema bem condicionado.
Exercícios
b) Não existe reta que interpola os pontos {(1, 1), (2, 2,1), (3, 3)}.
c) Não existe parábola de equação y = a0 +a1 x+a2 x2 que interpola dois pontos
dados {(x1 , y1 ), (x1 , y2 )}, com y1 6= y2 . Mas, existem infinitas parábolas de
equação x = a0 + a1 y + a2 y 2 que interpolam estes pontos.
p(x) = a1 + a2 (x − x1 ) + a3 (x − x1 )(x − x2 ) + · · ·
+ an (x − x1 )(x − x2 ) · · · (x − xn−1 ).
a1 = y1
a1 + a2 (x2 − x1 ) = y2
a1 + a2 (x3 − x1 ) + a3 (x3 − x1 )(x3 − x2 ) = y3
..
.
a1 + a2 (xn − x1 ) + · · · + an (xn − x1 ) · · · (xn − xn−1 ) = yn
a1 = y1
y2 − a1 y2 − y1
a2 = =
x2 − x1 x2 − x1
y3 −y2 y2 −y1
y3 − a2 (x3 − x1 ) − a1 (x3 −x2 )
− (x2 −x1 )
a3 = =
(x3 − x1 )(x3 − x2 ) (x3 − x1 )
...
Note que os coeficientes são obtidos por diferenças das ordenadas divididas
por diferenças das abscissas dos pontos dados. Para vermos isso mais claramente,
Tabela 6.1: Esquema de diferenças divididas para um conjunto com três pontos
{(xi , yi )}3i=1 .
j xj f [xj ] f [xj−1 ,xj ] f [xj−2 ,xj−1 ,xj ]
1 x1 f [x1 ] = y1
f [x2 ] − f [x1 ]
f [x1 ,x2 ] =
x2 − x1
f [x2 ,x3 ] − f [x1 ,x2 ]
2 x2 f [x2 ] = y2 f [x1 ,x2 ,x3 ] =
x3 − x1
f [x3 ] − f [x2 ]
f [x2 ,x3 ] =
x3 − x2
3 x3 f [x3 ] = y3
f [xj ] := yj
f [xj+1 ] − f [xj ]
f [xj , xj+1 ] :=
xj+1 − xj
f [xj+1 , xj+2 ] − f [xj , xj+1 ]
f [xj , xj+1 , xj+2 ] :=
xj+2 − xj
..
.
f [xj+1 , xj+2 , . . . , xj+k ] − f [xj , xj+1 , . . . , xj+k−1 ]
f [xj , xj+1 , . . . , xj+k ] :=
xj+k − xj
Chamamos f [xj ] de diferença dividida de ordem zero (ou primeira diferença divi-
dida), f [xi ,xj +1] de diferença dividida de ordem 1 (ou segunda diferença dividida)
e assim por diante.
Uma inspeção cuidadosa dos coeficientes obtidos em (6.2) nos mostra que
j xj f [xj ] f [xj−1 ,xj ] f [xj−2 ,xj−1 ,xj ] f [xj−3 ,xj−2 ,xj−1 ,xj ]
1 −1 3
1−3
= −2
0 − (−1)
2 − (−2)
2 0 1 =2
1 − (−1)
3−1 6−2
=2 =1
1−0 3 − (−1)
20 − 2
3 1 3 =6
3−0
43 − 3
= 20
3−1
4 3 43
Portanto, o polinômio interpolador do conjunto de pontos dados é
Assim, o polinômio p(x) de grau n − 1 que interpola os pontos dados, i.e. p(xj ) =
yj , j = 1, . . . ,n é dado por
n
X
p(x) = y1 L1 (x) + y2 L2 (x) + · · · + yn Ln (x) = yk Lk (x).
k=1
Portanto,
n
Y (x − xj )
Lk (x) =
j=1 (xk − xj )
j6=i
R
Exemplo 6.4.2. Considere o problema de aproximar o valor da integral 01 f (x)dx
pelo valor da integral do polinômio P (x) que coincide com f (x) nos pontos x0 = 0,
x1 = 12 e x2 = 1. Use a fórmula de Lagrange para encontrar P (x). Obtenha o
R
valor de 01 f (x)dx e encontre uma expressão para o erro de truncamento.
O polinômio interpolador de f (x) é
2 3 3 2 1
1 4 3 1
Z 1
2
P (x)dx = f (0) x − x + x + f − x + 2x
0 3 2 0 2 3 0
2 3 1 2 1
+ f (1) x − x
3 2 0
2 3 1 4 2 1
= f (0) − +1 +f − + 2 + f (1) −
3 2 2 3 3 2
1 2 1 1
= f (0) + f + f (1)
6 3 2 6
Para fazer a estimativa de erro usando o teorema 6.4.1 e temos
Z 1 Z 1 Z 1
f (x)dx − P (x)dx =
f (x) − P (x)dx
0 0 0
Z 1
≤ |f (x) − P (x)|dx
0
1
M Z 1
≤ x x − (x − 1) dx
6 "0
2
1
M Z 1/2
= x x− (x − 1)dx
6 0 2
#
1
Z 1
− x x− (x − 1)dx
1/2 2
M 1 1
M
= − − = .
6 64 64 192
Lembramos que M = maxx∈[0,1] |f ′′′ (x)|.
Observação 6.4.1. Existem estimativas melhores para o erro de truncamento
para este esquema de integração numérica. Veremos com mais detalhes tais esque-
mas na teoria de integração numérica.
Exemplo 6.4.3. Use o resultado do exemplo anterior para aproximar o valor das
seguintes integrais:
Z 1
a) ln(x + 1)dx
0
Z 1
2
b) e−x dx
0
Exercícios
(xi+1 − x) (x − xi )
Pi (x) = yi + yi+1
(xi+1 − xi ) (xi+1 − xi )
Exemplo 6.5.1. Construa uma função linear por partes que interpola os pontos
(0,0), (1,4), (2,3), (3,0), (4,2), (5,0).
A função procurada pode ser construída da seguinte forma:
0 x−1 + 1 x−0 ,0 ≤ x < 1
0−1 1−0
f (x) = 4 x−2
1−2
+ 3 x−1
2−1
,1 ≤ x < 2
3 x−3
2−3
+ 0 x−2
3−2
,2 ≤ x ≤ 3
Simplificando, obtemos:
,0 ≤ x < 1
x
f (x) = −x + 5 ,1 ≤ x < 2
−3x + 9 , 2 ≤ x ≤ 3
A figura 6.3 é um esboço da função f (x) obtida.
3.5
2.5
1.5
0.5
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
sj (x) = aj + bj (x − xj ) + cj (x − xj )2 + dj (x − xj )3 ,
aj = yj e bj = yj′ .
As duas últimas são obtidas pela solução do sistema de equações formado pelas
condições de contorno em xj+1 :
yj+1 − yj ′
yj+1 + 2yj′ yj+1 − yj ′
yj+1 + yj′
cj = 3 − e dj = −2 +
(xj+1 − xj )2 xj+1 − xj (xj+1 − xj )3 (xj+1 − xj )2
Esta relação entre o conjunto de valores para a derivada de um spline cúbico
{yj′ }j = 1n nos pontos de interpolação I e os coeficientes dos polinômios em cada
intervalo de interpolação pode ser resumida na seguinte proposição:
Proposição 6.6.1. Seja s um spline cúbico que interpola o conjunto de pontos
I = {(xj ,yj )}nj=1 ⊂ R2 tais que xj+1 > xj . Se {yj′ }nj=1 é o conjunto dos valores da
derivada de s em xj , então em cada intervalo [xj ,xj+1 ) (fechado também à direita
quando j = n − 1) o spline é igual a sj :
sj (x) = aj + bj (x − xj ) + cj (x − xj )2 + dj (x − xj )3 , (6.2)
onde
yj+1 − yj ′
yj+1 + 2yj′
aj = yj , cj = 3 − ,
h2j hj
′
+ yj′ (6.3)
yj+1 − yj yj+1
bj = yj , dj = −2
′
+
h3j h2j
e
hj = xj+1 − xj , j = 1,2, . . . ,n − 1 (6.4)
é a distância entre as abscissas de dois pontos de interpolação consecutivos.
De acordo com a proposição anterior, toda informação sobre um spline cúbico é
armazenada no conjunto {(xj ,yj ,yj′ )}nj=1 . Por construção, uma função s definida a
partir de (6.2), (6.3) e (6.4) com um conjunto {(xj ,yj ,yj′ )}nj=1 ⊂ R3 , onde xj+1 > xj
é de classe C 1 mas não necessariamente um spline cúbico. Para ser um spline
cúbico, os valores do conjunto {yj′ }nj=1 devem garantir a continuidade da derivada
segunda de s em todo intervalo (x1 ,xn ). Ou seja, devemos ter
Proposição 6.6.2. Dado o conjunto de pontos I = {(xj ,yj )}nj=1 ⊂ R2 tais que
xj+1 > xj , as derivadas de um spline cúbico que interpola os pontos I, yj′ , j =
1,2, . . . ,n satisfazem o sistema de equações algébricas lineares
!
′ yj − yj−1 yj+1 − yj
hj yj−1 + 2(hj−1 + hj )yj′ + ′
hj−1 yj+1 = 3 hj + hj−1 , (6.5)
hj−1 hj
ou seja,
y2 − y1
2y1′ + y2′ = 3
h1
. (6.7)
yn − yn−1
′
yn−1 + 2yn′ = 3
hn−1
y2 −y1
h1
y1′
h2 y2h−y
1
1
+ h1 y3h−y
2
2
y2′
h3 y3h−y 2
+ h2 y4h−y 3
y′ = e z= 3 (6.9)
..
2
..
3
.
. .
yn′
hn−1 yn−1hn−2
−yn−2
+ hn−2 ynh−y
n−1
n−1
yn −yn−1
hn−1
Exemplo 6.6.1. Construa um spline cúbico natural que passe pelos pontos (2, 4,5),
(5, − 1,9), (9, 0,5) e (12, − 0,5).
′
y1
1
3
(−1,9
− 4,5) −6,4
3 (−1,9 − 4,5) + 4 (0,5 − (−1,9))
′ 4 3
y2 −20,2
y=
′ e z= 3
3
=
.
y (0,5 − (−1,9)) + 4 (−0.5 − (0,5)) 1,4
3 4 3
y4′ 1
3
(−0,5 − (0.5)) −1
A solução é y1′ = −2,83̄, y2′ = −0,73̄, y3′ = 0,46̄ e y4′ = −0,73̄. Calculamos os
coeficientes usando as expressões (6.3):
c1 = 0 , d1 = 0,07̄,
c2 = 0,7 , d2 = −0,0916̄,
c3 = −0,4 , d3 = 0,04̄.
Portanto:
4,5 − 2,83̄(x − 2) + 0,07̄(x − 2)3 ,2 ≤ x < 5
S(x) =
−1,9 − 0,73̄(x − 5) + 0,7(x − 5)2 − 0,0916̄(x − 5)3 , 5 ≤ x < 9 .
0,5 + 0,46̄(x − 9) − 0,4(x − 9)2 + 0,04̄(x − 9)3
, 9 ≤ x ≤ 12
h2 y2h−y 1
+ h1 y3h−y 2
− h2 y1′
y2′ 1 2
h3 y3h−y 2
+ h2 y4h−y 3
y3′ 2
..
3
y′ = z= 3
.. e . .
.
hn−2 yn−2 −yn−3
hn−3 + hn−3 yn−1 −yn−2
hn−2
′
yn−1
hn−1 yn−1 −yn−2
hn−2 + hn−2 ynh−y
n−1
n−1
− hn−2 yn′
d1 = d2 e dn−2 = dn−1 ,
ou seja,
y2 − y1 y3 − y2
h22 y1′ + (h22 h21 )y2′ h21 y3′ =2 h22 − h21
− −
h1 h2
!
.
yn−1 − yn−2 yn − yn−1
h2n−1 yn−2
′
+ (h2n−1 − h2n−2 )yn−1
′
− h2n−2 yn′ = 2 h2n−1 − h2n−2
hn−2 hn−1
y2 − y1 y3 − y2
2 h22 − h21
h 1 h2
y1′
y2 −y1 y3 −y2
3 h2 h 1 + h 1 h 2
y2′
..
y′ =
z=
.
.. e
.
.
yn−1 −yn−2 yn −yn−1
3 hn−1 hn−2 + hn−2 hn−1
yn′
y n−1 − y n−2 y n − y n−1
2
2 hn−1 2 − hn−2
hn−2 hn−1
Se reduzirmos esse sistema pela eliminação das incógnitas e yn′ , o sistema re- y1′
sultante possui uma matriz de coeficientes diagonal dominante estrita, portanto,
a solução é única.
O termo not-a-knot (não nó) relaciona-se à nomenclatura dos splines. O termo
nó é utilizado para os pontos interpolados. Neles, a derivada terceira da função
spline é descontínua, portanto, quando impomos a continuidade dessa derivada em
x2 e xn−1 é como se esses pontos deixassem de ser nós.
0
y1′
h2 y2h−y 1
+ h1 y3h−y 2
y2′ 1
..
2
y′ = z= 3
.. e . .
. y −y
h n−1
n−1 hn−2
n−2
+ hn−2 ynh−y n−1
yn′ n−1
hn−1 y2h−y
1
1
+ h1 ynh−y
n−1
n−1
Neste caso também, se reduzirmos esse sistema pela eliminação das incógnitas y1′
e yn′ , o sistema resultante possui uma matriz de coeficientes diagonal dominante
estrita, portanto, a solução é única.
Ajuste de curvas
Figura 7.1: Exemplo de um problema de ajuste de uma reta entre três pontos,
veja o exemplo 7.0.1.
N
X
min (f (xj ) − yj )2 ,
f ∈F
j=1
177
178 Cálculo Numérico
Exemplo 7.0.1. Dado o conjunto de pontos {(1, 1,2), (1,5, 1,3), (2, 2,3)} e a famí-
lia de retas f (x) = a + bx, podemos mostrar que f (x) = −0,05 + 1,1x é a reta que
melhor aproxima os pontos dados no sentido de mínimos quadrados. Os pontos e
a reta ajustada e são esboçados na figura 7.1.
seja mínimo.
Para tal, primeiro observamos que f (xj ) = a1 + a2 xj e, portanto, o resíduo
pode ser escrito explicitamente como uma função de a1 e a2 conforme a seguinte
expressão:
N
X
R(a1 ,a2 ) = (a1 + a2 xj − yj )2 .
j=1
Observamos que R(a1 ,a2 ) é uma forma quadrática e que seu mínimo ocorre
quando suas derivadas parciais primeiras são iguais a zero, isto é,
N
∂R ∂ X
= (a1 + a2 xj − yj )2 = 0,
∂a1 ∂a1 j=1
N
∂R ∂ X
= (a1 + a2 xj − yj )2 = 0.
∂a2 ∂a2 j=1
Ou seja,
N
X
2 (a1 + a2 xj − yj ) · 1 = 0,
j=1
N
X
2 (a1 + a2 xj − yj ) · xj = 0,
j=1
P
Observando que N j=1 1 = N , o sistema linear acima pode ser escrito na forma
matricial M a = w, isto é,
PN PN
N j=1 xj a1 j=1 yj
= P (7.1)
P PN .
N 2 N
j=1 xj j=1 xj a2 j=1 xj yj
| {z } | {z } | {z }
M a w
Este sistema linear de duas equações e duas incógnitas admite uma única so-
lução quando o determinante da matriz dos coeficientes for não nulo, isto é,
2
N
X N
X
N x2j − xj 6= 0
j=1 j=1
Por fim, observamos que o sistema M a = w descrito na equação (7.1) pode ser
reescrito na forma V T V a = V T y, onde V := [1 x] é a matriz dos coeficientes do
seguinte sistema linear sobre determinado:
a1 + a2 x1 = y1
a1 + a2 x2 = y2
.. (7.3)
.
a1 + a2 xN = yN
Se os pontos dados não são colineares, este sistema não têm solução. Mas, sempre
que pelo menos duas abscissas foram diferentes, M = V T V é uma matriz invertível
e (veja o exercício 7.1.1), então
−1
a = V TV V T y, (7.4)
nos fornece a chamada solução por mínimos quadrados do sistema (7.3). Note que
esta é uma forma de se obter os coeficientes a = (a1 , a2 ) equivalente àquela dada
em (7.2).
Exemplo 7.1.1. Retornemos ao exemplo 7.0.1. Isto é, dado o conjunto de pontos
{(1, 1,2), (1,5, 1,3), (2, 2,3)}, encontrar a função do tipo f (x) = a1 +a2 x que melhor
se ajusta os pontos dados no sentido de mínimos quadrados.
Solução. Usando as fórmulas em (7.2), obtemos
7,25 · 4,8 − 4,5 · 7,75
a1 = = −0,05,
3 · 7,25 − 20,25
3 · 7,75 − 4,5 · 4,8
a2 = = 1,1.
3 · 7,25 − 20,25
Ou seja, verificamos que, de fato, a função f (x) = −0,05 + 1,1x corresponde à reta
que melhor ajusta os pontos dados no sentido de mínimos quadrados. Os pontos
e a reta ajustada estão esboçados na figura 7.1.
Exercício resolvido
ER 7.1.1. a) Mostre que o sistema linear M a = w descrito na equação 7.1 pode
ser reescrito na forma V T V a = V T y, onde V = [1 x].
b) Mostre que V , como definido no item a), tem posto igual a 2 quando pelo menos
duas abscissas do conjunto de pontos {(xj , yj )}N
j=1 são diferentes. E, portanto,
M = V V é uma matriz invertível.
T
e
y1
PN
1 1··· 1 y2 j=1 yj
V Ty = = P = w.
.
x1 x2 ··· .
xN . N
xj yj
j=1
yN
Exercícios
E 7.1.2. Seja dado o conjunto de pontos {(−0,35, 0,2), (0,15, −0,5), (0,23, 0,54),
(0,35, 0,7)}. Encontre a função f (x) = a1 + a2 x que melhor se ajusta no sentido
de mínimos quadrados aos pontos dados. Faça, então, um gráfico com os pontos e
o esboço da função ajustada.
E 7.1.3. Seja dado o conjunto de pontos {(−1,94, 1,02), (−1,44, 0,59), (0,93, −0,28),
(1,39, −1,04)}. Encontre a função f (x) = a1 + a2 x que melhor se ajusta no sentido
de mínimos quadrados aos pontos dados. Então, responda cada item:
a) Encontre o valor de f (1).
b) Encontre o valor de f (0,93).
c) Encontre o valor de |f (0,93) − (−0,28)|.
PN
d) Encontre o valor do resíduo R = j=1 (f (xj ) − yj )2 .
Forneça os valores calculados com 7 dígitos significativo por arredondamento.
minimiza o resíduo
n
X
R= [f (xi ) − yi ]2 .
i=1
Pm
Do fato que f (xi ) = j=1 aj fj (xi ), temos que cada resíduo pode ser escrito como
2
m
X
Ri = aj fj (xi ) − yi .
j=1
E os vetores a e w, por:
n
P
i=1 f1 (xi )yi
a1 n
P
f2 (xi )yi
i=1
a2
n
a= e w=
P
..
f3 (xi )yi
. i=1
..
.
am
n
P
fm (xi )yi
i=1
Observação 7.2.1. Este problema é equivalente a resolver pelo métodos dos mí-
nimos quadrados o seguinte sistema linear:
f1 (x1 ) f2 (x1 ) · · · fm (x1 ) y1
a1
f1 (x2 ) f2 (x2 ) · · · fm (x2 )
y2
a2
f (x )
1 3 f2 (x3 ) · · · fm (x3 )
.
= y
3
..
.. .. ... .. ..
. . . .
am
f1 (xn ) f2 (xn ) · · · fm (xn ) yn
Exemplo 7.2.1. Encontre a reta que melhor se ajusta aos pontos dados na se-
guinte tabela:
i 1 2 3 4 5
xi 0,01 1,02 2,04 2,95 3,55
yi 1,99 4,55 7,20 9,51 10,82
Exemplo 7.2.2. Encontre a função f (x) = a1 sen (πx) + a2 cos(πx) que melhor se
ajusta pelo critérios dos mínimos quadrados aos seguintes pontos dados
i 1 2 3 4 5
xi 0,00 0,25 0,50 0,75 1,00
yi −153 64 242 284 175
>>> xi = np.array([0,0.25,0.5,0.75,1])
>>> yi = np.array([-153,64,242,284,175])
>>> V = np.array([np.sin(np.pi*xi),np.cos(np.pi*xi)]).transpose()
>>> a = ((np.linalg.inv((V.transpose()).dot(V))).dot(V.transpose())).dot(yi)
>>> x = np.linalg.inv(A).dot(b)
>>> x = np.linalg.solve(A,b)
>>> np.linalg.lstsq(A,b)
p(x) = a1 + a2 x + · · · + am xm−1 .
Neste caso, a matriz V associada ao ajuste dos pontos {(x1 , y1 ), (x2 , y2 ), (x3 , y3 ),
. . ., (xn ,yn )} é dada por:
1 x1 x21 ··· x1m−1
1 x2 x22 · · · x2m−1
V = 1
x3 x23 · · · x3m−1
.. .. .. .
. . . ..
1 xn x2n · · · xnm−1
Exemplo 7.2.3. Entre o polinômio de grau 2 que melhor se ajusta aos pontos
dados na seguinte tabela:
i 1 2 3 4 5
xi 0,00 0,25 0,50 0,75 1,00
yi −153 64 242 284 175
Solução. Um polinômio de grau 2 pode ser escrito na seguinte forma:
p(x) = a1 + a2 x + a3 x2 .
Assim, o problema se resume em encontrarmos a solução por mínimos quadrados
do seguinte sistema linear:
a1 + a2 x1 + a3 x21 = y1
a2 + a2 x2 + a3 x22 = y2
a3 + a2 x3 + a3 x23 = y3
a4 + a2 x4 + a3 x24 = y4
a5 + a2 x5 + a3 x25 = y5
Ou, escrita na forma matricial, V a = y, onde:
1 x1 x21
1 x2 x22
V = 1
2
x3 x3
1 x4 x24
1 x5 x25
>>> xi = np.array([0,0.25,0.5,0.75,1])
>>> yi = np.array([-153,64,242,284,175])
>>> V = np.array([xi**2,xi**1,xi**0]).transpose()
>>> a = ((np.linalg.inv((V.transpose()).dot(V))).dot(V.transpose())).dot(yi)
Para fazermos o gráfico do polinômio e dos pontos, digitamos:
>>> xx = np.linspace(-0.25,1.25)
>>> plt.plot(xi,yi,'ro',xx,np.polyval(a,xx),'b-')
>>> plt.grid();plt.show()
♦
Exercícios
i 1 2 3 4
xi −1,50 −0,50 1,25 1,50
yi 1,15 −0,37 0,17 0,94
i 1 2 3 4 5
xi 0,01 1,02 2,04 2,95 3,55
yi 1,99 4,55 7,20 9,51 10,82
xi 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0
yi 31 35 37 33 28 20 16 15 18 23 31
∂R
= 2(Aex1 b − y1 )ex1 b + 2(Aex2 b − y2 )ex2 b + 2(Aex3 b − y3 )ex3 b = 0
∂A
∂R
= 2Ax1 (Aex1 b − y1 )ex1 b + 2Ax2 (Aex2 b − y2 )ex2 b
∂b
+ 2Ax3 (Aex3 b − y3 )ex3 b = 0
que é não linear em A e b. Esse sistema pode ser resolvido pelo método de Newton-
Raphson, o que pode se tornar custoso, ou mesmo inviável quando não dispomos
de uma boa aproximação da solução para inicializar o método.
Felizmente, algumas famílias de curvas admitem uma transformação que nos
leva a um problema linear. No caso da curva y = Aebx , observe que ln y = ln A+bx.
Assim, em vez de ajustar a curva original y = Aebx a tabela de pontos, ajustamos
a curva submetida a transformação logarítmica
ỹ := a1 + a2 x̃ = ln A + bx.
onde
1 x1
A=
1 x2
1 x3
2
A soma do quadrado dos resíduos.
Temos
1 1
A=
1 2
1 3
xi yi
0,0 9,12
0,1 1,42
0,2 - 7,76
0,3 - 11,13
0,4 - 11,6
0,5 - 6,44
0,6 1,41
0,7 11,01
0,8 14,73
0,9 13,22
1,0 9,93
Solução. Usando o fato que y = A cos(2πx + φ) = a cos(2πx) − b sen (2πx), onde
a = A cos(φ) e b = A sen (φ), z = [ a b ]T é solução do problema
B T Bz = B T y,
onde
1. 0.
0,8090170
−0,5877853
0,3090170 −0,9510565
−0,3090170 −0,9510565
cos(2πx0 ) − sen (2πx0 )
−0,8090170 −0,5877853
cos(2πx1 ) − sen (2πx1 )
B =
..
=
−1,0000000 0,0000000 .
.
−0,8090170 0,5877853
cos(2πx10 ) − sen (2πx10 )
−0,3090170 0,9510565
0,3090170 0,9510565
0,8090170 0,5877853
1,0000000 0,0000000
Observe que
A2 = 7,96147042 + 11,4057212
11,405721
sen (φ) = = 0,8199923
13,909546
φ = 0,9613976
AT Az = AT Y,
onde
1 x1
1 0,0
1 x 1 0,2
2
1 x3 1 0,4
A=
=
1 1 0,6
x4
1 x5
1 0,8
1 x6 1 1,0
e
1/y1
0,0099010
1/y 0,0117647
2
1/y3 0,0133333
Y =
=
1/y4 0,0151515
1/y5
0,0166667
1/y6 0,0181818
Derivação Numérica
197
198 Cálculo Numérico
′ 2f (ξ)
′′
f (x0 + h) = f (x0 ) + hf (x0 ) + h , h > 0, ξ ∈ (x0 ,x0 + h).
2
1
Uma função suave é uma função infinitamente continuamente diferenciável, isto é, f ∈
∞
C (R). Uma análise mais cuidadosa, revela que hipóteses mais fracas podem ser assumidas.
f (x0 + h) − f (x0 )
D+,h f (x0 ) :=
h
é de ordem h.
f ′′ (ξ)
f (x0 − h) = f (x0 ) − hf ′ (x0 ) + h2 , h > 0, ξ ∈ (x0 , x0 + h).
2
Isolando f ′ (x0 ), obtemos
f (x0 ) − f (x0 − h)
D−,h f (x0 ) := ,
h
que possui erro de truncamento de ordem h.
f (x0 + h) − f (x0 − h)
D0,h f (x0 ) := ,
2h
é uma aproximação para f ′ (x0 ) com erro de truncamento de ordem h2 , ou sim-
plesmente ordem 2.
1
Exemplo 8.1.2. Calcule a derivada numérica da função f (x) = e 2 x no ponto
x = 2 usando a diferença progressiva, diferença regressiva e diferença central com
h = 10−1 , h = 10−2 e h = 10−4 . Também, calcule o erro absoluto da aproximação
obtida em cada caso.
f (x + h) − f (x)
D+,h f (x) = D+,h f (x)(1 + ε(x,h)) = (1 + ε(x,h)). (8.3)
h
Também, consideremos
|f (x + h) − f (x + h)| = δ(x,h) ≤ δ
e
|f (x) − f (x)| = δ(x,0) ≤ δ,
onde f (x + h) e f (x) são as representações em ponto flutuante dos números f (x+h)
e f (x), respectivamente.
Então, da equação (8.3), a diferença do valor da derivada e sua aproximação
representada em ponto flutuante pode ser estimada por:
′
′ f (x + h) − f (x)
f (x) − D+,h f (x) = f (x) − (1 + ε(x,h)) .
h
Exercícios Resolvidos
ER 8.1.1. Aproxime a derivada de f (x) = sen (2x) − x2 no ponto x = 2 usando a
fórmula de diferenças finitas progressiva de ordem 1 com: a) h = 0,1 e b) h = 0,01.
Compute, também, o erro absoluto de cada aproximação computada.
#derivada analitica
def fl(x):
return 2*np.cos(2*x) - 2*x
#h=0.1
dy = dp1(f,2)
print("D.F. Progressiva de ordem 1 com h = %f" % 1e-1)
print("Df = %f" % dy)
print("Erro = %1.2e" % np.abs(fl(2)-dy))
#h=0.01
dy = dp1(f,2,1e-2)
print("D.F. Progressiva de ordem 1 com h = %f" % 1e-2)
print("Df = %f" % dy)
print("Erro = %1.2e" % np.abs(fl(2)-dy))
♦
Exercícios
0,0 0,50 1,00 1,50 2,00 2,50 3,00 3,50 4,00 4,50 5,00
0,0 1,05 1,83 2,69 3,83 4,56 5,49 6,56 6,11 7,06 8,29
∂vo
.
∂vi
Calcule o ganho quando vi = 1 e vi = 4.5 usando as seguintes técnicas:
Caso a b c d
vi = 1
vi = 4.5
2
E 8.1.5. Estude o comportamento da derivada de f (x) = e−x no ponto
x = 1,5 quando h fica pequeno.
h2 ′′ h3
f (x0 + h) = f (x0 ) + hf ′ (x0 ) + f (x0 ) + f ′′′ (x0 ) + O(h4 )
2 6
h2 ′′ h3
f (x0 − h) = f (x0 ) − hf ′ (x0 ) + f (x0 ) − f ′′′ (x0 ) + O(h4 ).
2 6
Somando as duas expressões, temos:
Exercícios
(x − x1 )(x − x2 ) (x − x0 )(x − x2 )
f (x) = f (x0 ) + f (x1 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 )
(x − x0 )(x − x1 ) f (ξ(x))
′′′
+ f (x2 ) + (x − x0 )(x − x1 )(x − x2 ).
(x2 − x0 )(x2 − x1 ) 6
2x − x1 − x2 2x − x0 − x2
f ′ (x) = f (x0 ) + f (x1 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 )
2x − x0 − x1
+ f (x2 )
(x2 − x0 )(x2 − x1 )
(8.5)
f (ξ(x))
′′′
+ ((x − x1 )(x − x2 ) + (x − x0 )(2x − x1 − x2 ))
6 !
f ′′′ (ξ(x))
+ Dx (x − x0 )(x − x1 )(x − x2 ).
6
2x0 − x1 − x2 2x0 − x0 − x2
f ′ (x0 ) = f (x0 ) + f (x1 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 )
2x0 − x0 − x1
+ f (x2 )
(x2 − x0 )(x2 − x1 )
f (ξ(x0 ))
′′′
+ ((x0 − x1 )(x0 − x2 ) + (x0 − x0 )(2x0 − x1 − x2 ))
6 !
f ′′′ (ξ(x0 ))
+ Dx (x0 − x0 )(x0 − x1 )(x0 − x2 ).
6
−3h −2h
f ′ (x0 ) = f (x0 ) + f (x1 )
(−h)(−2h) (h)(−h)
−h f (ξ(x0 ))
′′′
+ f (x2 ) + ((−h)(−2h))
(2h)(h) 6
1 3 1 f ′′′ (ξ(x0 ))
= − f (x0 ) + 2f (x1 ) − f (x2 ) + h2
h 2 2 3
1 1 1 f ′′′ (ξ(x1 ))
′
f (x1 ) = − f (x0 ) + f (x2 ) + h2
h 2 2 6
1 1 3 f ′′′ (ξ(x2 ))
f ′ (x2 ) = f (x0 ) − 2f (x1 ) + f (x2 ) + h2
h 2 2 3
1 3 1 f ′′′ (ξ(x0 ))
′
f (x0 ) = − f (x0 ) + 2f (x0 + h) − f (x0 + 2h) + h2
h 2 2 3
1 1 1 ′′′
(ξ(x +
f 0 h))
f ′ (x0 + h) = − f (x0 ) + f (x0 + 2h) + h2
h 2 2 6
1 1 3 f ′′′ (ξ(x0 + 2h))
′
f (x0 + 2h) = f (x0 ) − 2f (x0 + h) + f (x0 + 2h) + h2
h 2 2 3
ou ainda
1 f ′′′ (ξ(x0 ))
f ′ (x0 ) = [−3f (x0 ) + 4f (x0 + h) − f (x0 + 2h)] + h2 (8.6)
2h 3
1 f ′′′
(ξ(x 0 ))
f ′ (x0 ) = [f (x0 + h) − f (x0 − h)] + h2 (8.7)
2h 6
1 f ′′′ (ξ(x0 ))
f ′ (x0 ) = [f (x0 − 2h) − 4f (x0 − h) + 3f (x0 )] + h2 (8.8)
2h 3
Observe que uma das fórmulas é exatamente as diferenças centrais obtida anteri-
ormente.
Analogamente, para construir as fórmulas de cinco pontos tomamos o polinômio
de Lagrange para cinco pontos e chegamos a cinco fórmulas, sendo uma delas a
seguinte:
′ 1 h4 (5)
f (x0 ) = [f (x0 − 2h) − 8f (x0 − h) + 8f (x0 + h) − f (x0 + 2h)]+ f (ξ(x0 ))
12h 30
(8.9)
2
Exemplo 8.3.1. Calcule a derivada numérica de f (x) = e−x em x = 1,5 pelas
fórmulas de três e cinco pontos para h = 0,1, h = 0,01 e h = 0,001.
Exercícios
Em construção ... Gostaria de colaborar na escrita deste livro? Veja como em:
http://www.ufrgs.br/numerico
(x)′x∗ = 1 = c1 x1 + c2 x2 + . . . + cn xn (8.14)
(x2 )′x∗ = 2x∗ = c1 x21 + c2 x22 + . . . + cn x2n (8.15)
(x3 )′x∗ = 3(x∗ )2 = c1 x31 + c2 x32 + . . . + cn x3n (8.16)
.. .
. = .. (8.17)
(xn−1 )′x∗ = (n − 1)(x ) ∗ n−2
= c1 xn−1
1 + c1 xn−1
1 + . . . + cn xn−1
n (8.18)
Considere a base polinomial {φ1 (x), φ2 (x), φ3 (x)} = {1, x, x2 } e substitua f (x) por
φk (x) obtendo
(1)′x=0 = 0 = c1 (1) + c2 (1) + c3 (1) (8.20)
′
(x)x=0 = 1 = c1 (−h) + c2 (0) + c3 (h) (8.21)
(x2 )′x=0 = 0 = c1 (−h)2 + c2 (0)2 + c3 (h)2 (8.22)
que pode ser escrito na forma matricial
1 1 1 c1 0
−h 0
h = 1
(8.23)
c2
h2 0 h2 c3 0
Resolvendo o sistema obtemos {c1 , c2 , c3 } = {− 2h
1
, 0, 2h
1
} fornecendo a regra
1 1
f ′ |x=x1 ≈ − f (x1 ) + f (x3 ) (8.24)
2h 2h
f (x3 ) − f (x1 )
≈ (8.25)
2h
♦
Exercícios resolvidos
Em construção ... Gostaria de colaborar na escrita deste livro? Veja como em:
http://www.ufrgs.br/numerico
Exercícios
E 8.4.4. Seja [x0 ,x1 , . . . ,x4 ] = [−2h, − h,0,h,2h] e x∗ = 0, obtenha uma regra
de diferenciação para aproximar f ′′ (x∗ ).
x y
0 1,95
1 1,67
2 3,71
3 3,37
4 5,12
5 5,79
6 7,50
7 7,55
8 9,33
9 9,41
10 11,48
12
11
10
1
0 1 2 3 4 5 6 7 8 9 10
Exercícios resolvidos
Em construção ... Gostaria de colaborar na escrita deste livro? Veja como em:
http://www.ufrgs.br/numerico
Exercícios
Em construção ... Gostaria de colaborar na escrita deste livro? Veja como em:
http://www.ufrgs.br/numerico
Integração Numérica
n−1
X
I≈S= ∆Si .
i=1
216
217
4.5
3.5
2.5
1.5
0.5
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Z 2
(x2 + 1) dx
0
h1 = 2 h1 f (1) = 4
h2 = 1 h2 f (0,5) + h2 f (1,5) = 4,5
h3 = 0,5 4,625
h4 = 0,25 4,65625
Observe que:
" #2
Z 2
2 x3 8
(x + 1) dx = +x = + 2 = 4,6666667.
0 3 0
3
Nos códigos Python apresentados ao longo deste capítulo, estaremos assumindo:
>>> from __future__ import division
>>> import numpy as np
1
Utilizaremos neste capítulo a notação fi para indicar f (xi ).
onde Z b
Ai = Li (x) dx. (9.6)
a
∆Si ≈ f (xi )h
tal que a área total será aproximada pelas somas de Riemann à esquerda
n−1
X n−1
X
S= ∆Si = f (xi )h
i=1 i=1
Uma terceira opção é utilizar o ponto médio do intervalo [xi ,xi+1 ] o qual fornece
a regra do ponto médio
n−1
X xi + xi+1
S= f (ξi )h, ξi = . (9.9)
i=1 2
onde
" # x2
Z b
x − x1 (x − x1 )2
A1 = dx =
a x2 − x1 2h x1
(x2 − x1 ) h 1 2 2
= = = h
2h 2h 2
Da mesma forma,
Z b
(x − x2 ) 1
A2 = dx = h
a (x1 − x2 ) 2
1 1
Z b
f (x) dx ≈ f1 + f2 h. (9.11)
a 2 2
a+b
x1 := a, x2 := e x3 := b
2
onde
Z b
Ai = Li (x) dx (9.15)
a
1 4 1
Z b
f (x) dx = f (x1 ) + f (x2 ) + f (x3 ) h.
a 6 6 6
+ (x − x2 )
24 x1
1 Z x3 (4)
+ f (ξ(x))(x − x2 )4 dx,
24 x1
Exercícios
R 2
E 9.1.3. Calcule numericamente o valor de 25 e4−x dx usando os métodos
compostos do ponto médio, trapézio e Simpson. Obtenha os resultados utilizando,
em cada quadratura, o número de pontos indicado.
Z b
f (t) dt ≈ w1 f (t1 ) + w2 f (t2 ) + . . . + wn f (tn ) (9.16)
a
2
Por exemplo, se n = 2, então a regra é exata para retas.
1 4 1
Resolvendo o sistema obtemos [w1 ,w2 ,w3 ] = h[ , , ] fornecendo a regra de
6 6 6
Simpson Z h
h 4h h
f (t) dt ≈ f0 + f1 + f2 (9.31)
0 6 6 6
♦
Z b Ni
hX
f (x) dx ≈ [f (xk ) + f (xk+1 )]
a 2 k=1
h
= [f (x1 ) + 2f (x2 ) + 2f (x3 ) + · · · + 2f (xNi ) + f (xNi +1 )]
2
Ni
h X
= [f (x1 ) + f (xNi +1 )] + h f (xi )
2 i=2
Z b n Z xk+1
X
f (x) dx = f (x) dx
a k=1 xk
n
xk+1 + xk
X xx+1 − xk
≈ f (xk ) + 4f + f (xk+1 )
k=1 6 2
Z b " n−1 n
#
h X X
f (x) dx ≈ f (x1 ) + 2 f (x2i+1 ) + 4 f (x2i ) + f (x2n+1 ) + O(h5 )
a 3 i=1 i=1
Z 2
2
x2 ex dx
0
Exercícios
R1
E 9.3.2. O valor exato da integral imprópria 0 x ln(x) dx é dado por
!1
Z 1
x2 x2
x ln(x) dx = ln x −
= −1/4
0 2 4
0
b) Use a identidade
Z 1 Z 1 Z 1
ln(x) sen (x) dx = ln(x)x dx + ln(x) [sen (x) − x] dx
0 0 0
!1
x2 x2
Z 1
= ln x − + ln(x) [sen (x) − x] dx
2 0
4
0
1 Z1
= − + ln(x) [sen (x) − x] dx
4 0
R
e aproxime a integral 01 ln(x) [sen (x) − x] dx numericamente via Gauss-Legendre
com n = 2, n = 3, n = 4, n = 5, n = 6 e n = 7.
Defina I(h) a aproximação desta integral pelo método dos trapézios composto com
malha de largura constante igual a h. Aqui h = b−a
Ni
para algum Ni inteiro, i.e.:
Ni
h X b−a
I(h) = f (a) + 2 f (xj ) + f (b) , Ni =
2 j=2 h
Teorema 9.4.1. Se f (x) é uma função analítica no intervalo (a,b), então a função
I(h) admite uma representação na forma
I(h) = I0 + I2 h2 + I4 h4 + I6 h6 + . . .
I(h) = I0 + I2 h2 + I4 h4 + I6 h6 + . . .
!
h h2 h4 h6
I = I0 + I2 + I4 + I6 + . . .
2 4 16 64
4I(h/2) − I(h) h h
= [f (a) + 2f (c) + f (b)] − [f (a) + f (b)]
3 3 6
h
= [f (a) + 4f (c) + f (b)]
6
Observe que esquema coincide com o método de Simpson.
A partir de agora, usaremos a seguinte notação
R1,1 = I(h)
R2,1 = I(h/2)
R3,1 = I(h/4)
..
.
Rn,1 = I(h/2n−1 )
Observamos que os pontos envolvidos na quadratura Rk,1 são os mesmos pontos
envolvidos na quadratura R(k − 1,1) acrescidos dos pontos centrais, assim, temos
a seguinte fórmula de recorrência:
k−2 !
1 h 2X h
Rk,1 = Rk−1,1 + k−1 f a + (2i − 1) k−1
2 2 i=1 2
Definimos Rk,2 para k ≥ 2 como o esquema de ordem quatro obtido da fórmula
do exemplo 9.4.1:
4Rk,1 − Rk−1,1
Rk,2 =
3
Os valores Rk,2 representam então os valores obtidos pelo método de Simpson
composto aplicado a uma malha composta de 2k−1 + 1 pontos.
Similarmente os valores de Rk,j são os valores obtidos pela quadratura de ordem
2j obtida via extrapolação de Richardson. Pode-se mostrar que
Rk,j−1 − Rk−1,j−1
Rk,j = Rk,j−1 + .
4j−1 − 1
Exemplo
R 2 −x2
9.4.2. Construa o esquema de Romberg para aproximar o valor de
0 e dx com erro de ordem 8.
55,59815 0,000000 0,000000 0,000000
30,517357 22,157092 0,000000 0,000000
O que nos fornece os seguintes resultados:
20,644559 17,353626 17,033395 0,000000
17,565086 16,538595 16,484259 16,475543
Ou seja, temos: Z 2
2
ex dx ≈ 16,475543
0
usando uma aproximação de ordem 8.
Exemplo
R 2 2 x2
9.4.3. Construa o esquema de Romberg para aproximar o valor de
0 x e dx com erro de ordem 12.
218,3926
111,91458 76,421909
66,791497 51,750469 50,105706
O que nos fornece:
51,892538 46,926218 46,604601 46,549028
47,782846 46,412949 46,378731 46,375146 46,374464
46,72661 46,374531 46,37197 46,371863 46,37185 46,371847
Ou seja, temos: Z 2
2
x2 ex dx ≈ 46,371847
0
com uma aproximação de ordem 12.
Exercícios
Rπ
E 9.4.2. Calcule os valores da quadratura de Romberg de R1,1 até R4,4 para
0 sen (x) dx. Não use rotinas prontas neste problema.
R2q
b) 0 2 − cos(x) dx
R2
c) 0
√ 1
dx
2−cos(x)
c) Método de Simpson
" ! #
Z b
a+b b−a
f (x) dx ≈ f (a) + 4f + f (b)
a 2 6
!
b−a 2(b − a) a+b b−a
= f (a) + f + f (b)
6 3 2 6
3
X
:= wj f (xj )
j=1
f (x) = a0 + a1 x + a2 x2 + . . . + an xn + Rn (x)
ou, equivalentemente:
n
X Z b
bk+1 − ak+1
wj xkj = xk dx = , k = 0,1, . . . n
j=1 a k+1
onde wj = b−a
2
, x1 = a e x2 = b.
Pn
(k = 0) : j=1 wj = b − a
Pn b2 −a2
(k = 1) : j=1 wj xj = (a + b) b−a
2
= 2
Pn b3 −a3
(k = 2) : j=1 wj x2j = (a2 + b) 2 =
2 b−a
6 3
onde w1 = w3 = b−a
6
,w2 = 4 b−a
6
, x1 = a, x2 = a+b
2
e x3 = b
Pn
(k = 0) : j=1 wj = (1 + 4 + 1) b−a
6
=b−a
Pn b2 −a2
(k = 1) : j=1 wj xj = (a + 4 a+b
2
+ b) b−a
6
= (a + b) b−a
2
= 2
2
Pn b3 −a3
(k = 2) : j=1 wj xj = (a + 4
2 2 a+b
2
+ b2 ) b−a
6
= 3
3
Pn b4 −a4
(k = 3) : j=1 wj xj = (a + 4
3 3 a+b
2
+ b3 ) b−a
6
= 4
4
Pn b5 −a5
(k = 4) : j=1 wj xj = (a + 4
4 4 a+b
2
+ b4 ) b−a
6
6= 4
w1 + w2 = 2
x1 w1 + x2 w2 = 0
2
x21 w1 + x22 w2 =
3
x31 w1 + x32 w2 = 0
√ √
e−1 + 4e0 + e1 − 3 3
e − e−1 e−1 + e e− 3 +e 3
ex 3
≈ 2,35040 ≈ 3,08616 ≈ 2,36205 ≈ 2,34270
16 4
√
2
p
9 − 9 2
x 3+ x3 3,41421 1,13807 1,15411
≈ 1,14924
3
e−e−1
x2 ex 3 ≈ 0,78347 3,08616 1,02872 0,67905
Exercícios
• Pode-se mostrar que este problema sempre tem solução e que a solução é
única se t1 < t2 < . . . < tn
n tj wj
1 0 2
√
3
2 ± 1
3
8
0
3 s
9
3 5
±
5 9
s
q √
± 3 − 2 6/5 /7 18+ 30
36
4 s
q √
± 3 + 2 6/5 /7 18− 30
36
128
0
v 225
u s √
5 1u 10 322 + 13 70
5−2
t
±
3 7 900
v
u s √
1u 10 322 − 13 70
5+2
t
±
3 7 900
n I
2 2,3094011
3 2,2943456
4 2,2957234
5 2,2955705
Em Python, temos:
def f(x):
return np.sqrt(1+x^2)
#G-L n=2
x2 = sqrt(3)/3
w2 = 1
I2 = w2[1]*f(x2[1]) + w2[1]*f(-x2[1])
print("%1.7f\n" % I2)
#G-L n=3
x3 = [0 -sqrt(3/5) sqrt(3/5)]
w3 = [8/9 5/9 5/9]
I3 = w3[1]*f(x3[1]) + w3[2]*f(x3[2]) + w3[2]*f(-x3[2]);
print("%1.7f\n" % I3)
#G-L n=4
x4 = [sqrt((3-2*sqrt(6/5))/7) sqrt((3+2*sqrt(6/5))/7)]
w4 = [(18+sqrt(30))/36 (18-sqrt(30))/36]
I4 = w4(1)*f(x4(1)) + w4(1)*f(-x4(1)) ...
+ w4(2)*f(x4(2)) + w4(2)*f(-x4(2))
print("%1.7f\n" % I4)
#G-L n=5
x5 = [0 1/3*sqrt(5-2*sqrt(10/7)) 1/3*sqrt(5+2*sqrt(10/7))]
w5 = [128/225 (322+13*sqrt(70))/900 (322-13*sqrt(70))/900]
I5 = w5(1)*f(x5(1)) + w5(2)*f(x5(2)) + w5(2)*f(-x5(2)) ...
+ w5(3)*f(x5(3)) + w5(3)*f(-x5(3))
print("%1.7f\n" % I5)
Mudança de intervalo
Os coeficientes da quadratura de Gauss-Legendre forma obtidos no intervalo
[−1,1]. Para aproximar a integral de f (x) no intervalo [a,b] devemos fazer a mu-
dança de variável
x̄i = αti + β, α = (b − a)/2, β = (b + a)/2
tal que
Z b n
X
f (x) dx ≈ wi f (x̄i )(b − a)/2
a i=1
Quando subdividimos o intervalo inicial [a,b] em N intervalos com extremos
[xi ,xi+1 ] a transformação torna-se
x̄i = αti + β, α = (xi+1 − xi )/2, β = (xi+1 + xi )/2
e Z xi+1 n
X
f (x) dx ≈ wi f (x̄i )(xi+1 − xi )/2
xi i=1
Exemplo 9.6.2. Aproximar
Z 1√
I= 1 + x2 dx
0
pelo método de Gauss-Legendre com 3 pontos.
Solução. Para tanto, fazemos a mudança de variáveis u = 2x − 1:
Z 1√
I= 1 + x2 dx
0
s
1Z 1 u+1 2
= 1+ du.
2 −1 2
E, então aplicamos a quadratura gaussiana nesta última integral, o que nos fornece
I ≈ 1,1478011. No GNU Octave, podemos computar estas aproximações com o
seguinte código:
def f(u):
return sqrt(1+(u+1)^2/4)/2
x3 = [0 -sqrt(3/5) sqrt(3/5)]
w3 = [8/9 5/9 5/9]
I3 = f(x3(1))*w3(1) + f(x3(2))*w3(2) + f(-x3(2))*w3(2)
print("%1.7f\n" % I3)
♦
Exercícios
R1
E 9.7.1. Considere o problema de calcular numericamente a integral I =
cos(x)
−1 f (x)dx quando f (x) =
√ .
|x|
R 1 ex R 1 e−x
E 9.7.3. Calcule as integrais 0 |x|1/4 dx e 0 |x|4/5 dx usando procedimentos
analíticos e numéricos.
E 9.7.4. Use a técnica de integração por partes para obter a seguinte identi-
dade envolvendo integrais impróprias:
Z ∞
cos(x) Z ∞
sen (x)
I= dx = dx.
0 1+x 0 (1 + x)2
c) Melhore sua cultura geral: A lei de Dulong-Petit para o calor específico dos só-
lidos precede a teoria de Debye. Verifique que a equação de Debye é consistente
com Dulong-Petit, ou seja:
lim Cv = 3N kB .
T →∞
Exemplo 10.0.2.
du
=u (10.4)
dt
u(0) = a (10.5)
Exemplo 10.0.3.
du
= sen (u2 + sen (t)) (10.6)
dt
u(0) = a (10.7)
248
10.1. TEORIA DE EQUAÇÕES DIFERENCIAIS 249
A solução do segundo exemplo é fácil de ser obtida: u(t) = aet . Porém como
podemos resolver o terceiro problema?
Muitos problemas de valor inicial da forma (10.1) não podem ser resolvidos
exatamente, ou seja, sabe-se que a solução existe e é única, porém não podemos
expressá-la em termos de funções elementares. Por isso é necessário calcular apro-
ximações numéricas. Diversos métodos completamente diferentes estão disponíveis
para aproximar uma função real.
Existem várias maneiras de obter aproximações para a solução deste problema.
Nos limitaremos a estudar métodos que aproximam u(t) em um conjunto finito de
valores de t chamado malha que será denotado por {ti }N i=1 = {t1 , t2 , t3 , . . . , tN }.
Desta forma, aproximamos a solução u(ti ) por ui em cada ponto da malha usando
diferentes esquemas numéricos.
Nos códgos em Python apresentados neste capítulo, estaremos assumindo que
as seguintes bibliotecas e módulos estão importados:
• A solução é única?
u′ (t) = 2u(t)
u(0) = 1
cuja solução é u(t) = e2t . O método de Euler aplicado a este problema produz o
esquema:
Suponha que queremos calcular o valor aproximado de u(1) com h = 0,2. Então
os pontos t(1) = 0, t(2) = 0,2, t(3) = 0,4, t(4) = 0,6, t(5) = 0,8 e t(6) = 1,0 formam
os seis pontos da malha. As aproximações para a solução nos pontos da malha
usando o método de Euler são:
u(0) ≈ u(1) =1
u(0,2) ≈ u(2) = (1 + 2h)u(1) = 1,4u(1) = 1,4
u(0,4) ≈ u(3) = 1,4u(2) = 1,96
u(0,6) ≈ u(4) = 1,4u(3) = 2,744
u(0,8) ≈ u(5) = 1,4u(4) = 3,8416
u(1,0) ≈ u(6) = 1,4u(5) = 5,37824
#define f(t,u)
def f(t,u):
return 2*u
#cria vetor t e u
t = np.empty(N)
u = np.copy(t)
#C.I.
t[0] = 0
u[0] = 1
#iteracoes
for i in np.arange(N-1):
t[i+1] = t[i] + h
u[i+1] = u[i] + h*f(t[i],u[i])
#imprime
for i,tt in enumerate(t):
print("%1.1f %1.4f" % (t[i],u[i]))
Exemplo 10.2.2. Aproxime a solução do PVI
du
= −0.5u + 2 + t (10.23)
dt
u(0) = 8 (10.24)
Itere a fórmula
%---------------------------
function [u,t]=euler(h,Tmax)
u(1)= 8;
t(1)= 0;
itmax = Tmax/h;
for n=1:itmax
t(n+1)= t(n) + h;
u(n+1)= u(n) + h*(-0.5*u(n)+2+t(n));
end
plot(t,u,'g*-');
%---------------------------
Tabela 10.1: Tabela comparativa entre método de Euler e solução exata para
problema 10.2.3.
t Exato Euler h = 0,1 Euler h = 0,01
0 1/2 0,5 0,5
e1/2
1/2 1+e1/2
≈ 0,6224593 0,6231476 0,6225316
1 e
1+e
≈ 0,7310586 0,7334030 0,7312946
e2
2 1+e2
≈ 0,8807971 0,8854273 0,8812533
3
3 e
1+e3
≈ 0,9525741 0,9564754 0,9529609
e
u = (1 − u)et
Colocando o termo u em evidência, encontramos:
(1 + et )u = et (10.27)
t
E, finalmente, encontramos a solução exata dada por u(t) = 1+e
e
t.
uk+1 = uk + huk (1 − uk ),
u1 = 1/2.
u′ (t) = −u(t) + t
u(0) = 1,
uk+1 = uk + h(−uk + tk )
u1 = 1
Comparação
Exercícios
Se u(tn+1 ) for aproximado por un+1 com erro da ordem O(hp+1 ) dizemos que o
método tem ordem de precisão p.
Queremos obter a ordem de precisão do método de Euler. Para isso, substituí-
mos a EDO u′ = f (t,u) na expansão em série de Taylor
u(tn+1 ) = u(tn ) + hu′ (tn ) + h2 u′′ (tn )/2 + O(h3 ) (10.28)
e obtemos
u(tn+1 ) = u(tn ) + hf (tn ,u(tn )) + h2 u′′ (tn )/2 + O(h3 ) (10.29)
Subtraindo (10.29) do método de Euler
un+1 = un + h f (tn ,un ) (10.30)
obtemos
en+1 = un+1 − u(tn+1 ) (10.31)
= un − u(tn ) + h(f (tn ,u(tn ) + en ) − f (tn ,u(tn ))) + (10.32)
h2 ′′
+ u + O(h3 ) (10.33)
2 n
Defina o erro numérico como en = un − u(tn ) onde u(tn ) é a solução exata e un
é a solução aproximada. Assim
h2 ′′
en+1 = en + h(f (tn ,u(tn ) + en ) − f (tn ,u(tn ))) + u
2 n
+ O(h3 ) (10.34)
Usando a condição de Lipschitz em f temos
h2 ′′
|en+1 | ≤ |en | + h|f (tn ,u(tn ) + en ) − f (tn ,u(tn ))| + |u | + O(h3 )(10.35)
2 n
h2 ′′
≤ |en | + hL|u(tn ) + en − u(tn )| + |u | + O(h3 ) (10.36)
2 n
h2 ′′
≤ |en | + hL|en | + |u | + O(h3 ) (10.37)
2 n
h2
≤ (1 + hL)|en | + |u′′n | + O(h3 ) (10.38)
2
ET G = nET L (10.39)
= n[h2 /2|u′′ | + O(h3 )] (10.40)
= T h/2|u′′ | + O(h2 ) (10.41)
ou seja
ET Gn+1
Euler = O(h)
10.3.1 Convergência
Um método é dito convergente se para toda EDO com f Lipschitz e todo
t > 0 temos que
lim |un − u(tn )| = 0, ∀n
h→0
Convergência significa que a solução numérica tende a solução do PVI.
Teorema 10.3.1. O método de Euler é convergente.
De fato, se f Lipschitz e |e0 | = 0, temos que
10.3.2 Consistência
Definição 10.3.1. Dizemos que um método numérico Rh (un ) = f é consistente
com o PVI u′ (t) = f se para qualquer u(t)
Isto é equivalente a
ET L
lim =0 (10.44)
h→0 h
Licença CC-BY-SA-3.0. Contato: [email protected]
258 Cálculo Numérico
10.3.3 Estabilidade
Definição 10.3.2. Um método numérico é estável se
|un − vn | ≤ C1 |u1 − v1 |, ∀n
Isto significa que dadas duas condições iniciais u1 e v1 , teremos que as soluções
un e vn estarão a uma distância limitada por uma constante C1 vezes |u1 − v1 |. Se
u1 e v1 estiverem próximas então un e vn estão também próximas dependendo da
constante C1 (obviamente C1 depende da função f ).
Considere o PVI linear bem-posto
Para que o método de Euler seja estável, é necessário que h seja escolhido tal que
|1 + hλ| < 1. Ou seja, hλ deve estar em DEuler onde
Note que este método é implícito (a equação é implícita) pois depende de un+1
dos dois lados da equação. Se a função f for simples o suficiente, podemos resolver
a equação isolando o termo un+1 . Se isso não for possível, devemos usar um dos
métodos vistos anteriormente para calcular as raízes da equação (por exemplo,
método da bissecção e método de Newton).
Pode ser mostrado que o erro de truncamento local é
ET Ln+1 2
EulImp = O(h ).
ET Gn+1
EulImp = O(h).
DEulImp onde
{z ∈ C : ℜz < 0} ⊆ D
Z t2
u(t2 ) = u(t1 ) + f (t,u(t)) dt (10.68)
t1
1 1
u2 = u1 + (t2 − t1 ) f (t1 ,u1 ) + f (t2 ,u2 ) (10.69)
2 2
motivando o método trapezoidal
h
un+1 = un + (f (tn ,un ) + f (tn+1 ,un+1 )) (10.70)
2
O método trapezoidal é dito implícito, pois para obter un+1 é necessário calcular
f (tn+1 ,un+1 ).
Entretanto, pode ser mostrado que o erro de truncamento local é
ET Ln+1 3
T rap = O(h )
ET Gn+1 2
T rap = O(h )
ET Ln+1 3
Heun = O(h )
ET Gn+1 2
Heun = O(h )
Exercícios
usando passo h = 0,1. Compare os valores da solução exata dada por u(t) = 2et1−1
com os numéricos nos pontos t = 0, t = 0.1, t = 0.2, t = 0.3, t = 0.4, t = 0.5,
t = 0.6, t = 0.7, t = 0.8, t = 0.9, t = 1.0.
h2 ′′ h3
u(t + h) = u(t) + hu′ (t) + u (t) + u′′′ (t) + . . . (10.76)
2! 3!
Utilizando dois termos temos o método de Euler. Utilizando os três primeiros
termos da série e substituindo u′ (t) = f (t,x) e u′′ (t) = ∂f
∂t
(t,x) temos o método
de Taylor de ordem 2
h2 ∂f
un+1 = un + hf (tn ,un ) + (tn ,un ) (10.77)
2! ∂t
O método de Taylor de ordem 3 é
h2 ∂f h3 ∂ 2 f
un+1 = un + hf (tn ,un ) + (tn ,un ) + (tn ,un )
2! ∂t 3! ∂t2
z2 z3 zp
p(z) = 1 + z + + + ... + (10.78)
2! 3! p!
para t ∈ [0,10].
a. Plote a solução para h = 0.16,0.08, 0.04, 0.02, 0.01 para o método de Taylor
de ordem 1, 2 e 3. (Plote todos de ordem 1 no mesmo gráfico, ordem 2 em
outro gráfico e ordem 3 outro gráfico separado.)
c. Fixe agora o valor h = 0.02 e plote no mesmo gráfico uma curva para cada
método.
e teremos
Ps
un+s = un+s−1 + h m=0 bm fn+m (10.87)
Para isso devemos obter [b3 ,b2 ,b1 ,b0 ] tal que o método seja exato para polinômios
até ordem 3. Podemos obter esses coeficientes de maneira análoga a obter os
coeficientes de um método para integração.
Supondo que os nós tk estejam igualmente espaçados, e para facilidade dos
cálculos, como o intervalo de integração é [tn+3 ,tn+4 ], translade tn+3 para a origem
tal que [tn ,tn+1 , . . . ,tn+4 ] = [−3h, − 2h, − h,0,h].
Considere a base [φ0 (t), . . . ,φ3 (t)] = [1,t,t2 ,t3 ] e substitua f (t) por φk (t) obtendo
Z h
1 dt = h = h(b0 (1) + b1 (1) + b2 (1) + b3 (1))
0
Z h h2
t dt = = h(b0 (0) + b1 (−h) + b2 (−2h) + b3 (−3h))
0 2
Z h
h3
t2 dt = = h(b0 (0)2 + b1 (−h)2 + b2 (−2h)2 + b3 (−3h)2 )
0 3
Z h
h4
t3 dt = = h(b0 (0)3 + b1 (−h)3 + b2 (−2h)3 + b3 (−3h)3 )
0 4
un+4 = un+3 + h
24
[55fn+3 − 59fn+2 + 37fn+1 − 9fn ] (10.94)
un+3 = un+2 + h
12
[23fn+2 − 16fn+1 + 5fn ] (10.96)
Para isso devemos obter [b3 ,b2 ,b1 ,b0 ] tal que o método seja exato para polinômios
até ordem 3. Podemos obter esses coeficientes de maneira análoga a obter os
coeficientes de um método para integração.
un+3 = un+2 + h
24
[9fn+3 + 19fn+2 − 5fn+1 + fn ] (10.102)
un un+1
| |
tn tn+1
τ1 τ2 ··· τν
Z tn+1
un+1 = un + f (t,u(t)) dt (10.115)
tn
Z 1
= un + h f (tn + hτ,u(tn + hτ )) dτ (10.116)
0
ν
X
= un + h bj f (tn + cj h,u(tn + cj h)) (10.117)
j=1
Por exemplo, se ν = 3 estágios teremos [τ0 ,τ1 ,τ2 ] = [tn + c0 h,tn + c1 h,tn + c2 h],
Uj ≡ u(τj ) e Fj ≡ f (τj ,Uj ), j = 1,2,3. Inicie com U1 = un (c1 = 0) como a solução
no passo anterior e aproxime U2 ,U3 , com uma combinação linear dos valores de Fj
anteriores, ou seja,
U1 = un (10.118)
U2 = un + ha21 F1 (10.119)
U3 = un + ha31 F1 + ha32 F2 (10.120)
un+1 = un + h[b1 F1 + b2 F2 + b3 F3 ] (10.121)
10.16.2 Método de RK ν = 2
Assumindo suavidade suficiente em f , expanda em série de Taylor
F2 = f (tn + c2 h,U2 ) (10.122)
= f (tn + c2 h,un + a21 hfn ) (10.123)
∂fn ∂fn
= fn + h[c2 + a21 fn ] + O(h2 ) (10.124)
∂t ∂u
fazendo com que (10.121) se torne
un+1 = un + h[b1 F1 + b2 F2 ] (10.125)
∂fn ∂fn
= un + h(b1 + b2 )fn + h2 b2 [c2 + a21 fn ] + O(h3 ) (10.126)
∂t ∂u
Usando a EDO e derivando-a obtemos
ut = f (t,u) (10.127)
utt = ft + fu ut = ft + fu f (10.128)
e expandindo em série de Taylor a solução exata em tn+1 ,
h2
u(tn+1 ) = un + hut + utt + O(h3 ) (10.129)
2
h2
= un + hfn + [ft + fu f ] + O(h3 ) (10.130)
2
e comparando com (10.126) obtemos as condições para ordem p ≥ 2,
1
b1 + b2 = 1, b2 c 2 = a21 = c2 (10.131)
2
O sistema possui mais de uma solução. Algumas escolhas comuns são
0 0 0
1
2
1
2
, 3 3
2 2
e 1 1
0 1 1
4
3
4
1
2
1
2
U1 = un (10.132)
U2 = un + hF1 (10.133)
1 1
un+1 = un + h[ F1 + F2 ] (10.134)
2 2
Note que o método é de ordem p = 2 pois os termos que sobraram são de O(h3 ).
Seguindo um procedimento similar, podemos obter as condições para um mé-
todo com ν = 3 e ordem p = 3, que são
b1 + b2 + b3 = 1, b2 c2 + b3 c3 = 1
2
(10.135)
1
b2 c22 + b3 c23 = , b3 a32 c2 = 16 (10.136)
3
Alguns exemplos de métodos de RK de 3 estágios são o método clássico de
Runge-Kutta
0
1 1
2 2
1 −1 2
1 4 1
6 6 6
e o método de Nystrom
0
2 2
3 3
2
3
0 2
3
2 3 3
8 8 8
Com paciência e a ajuda de um software algébrico (como Maple) é possível
encontrar um método de quarta ordem e ν = 4 estágios como
0
1 1
2 2
1
2
0 1
2
1 0 0 1
1 2 2 1
6 6 6 6
Portanto existe uma faixa hmin < h < hmax onde o método apresenta a ordem
desejada. Essa região depende do método e do PVI estudado.
Mas se estivermos nessa região podemos aproximar a ordem do método da se-
guinte forma: Considere a solução para um determinado t = T ∗ fixo, u(T ∗ ). Consi-
dere também as aproximações das soluções obtidas com espaçamento h, denotada
por uh ; a aproximação obtida com espaçamento dividido por 2, h/2, denotada por
uh/2 ; a aproximação obtida com espaçamento h/4, denotada por uh/4 ,. . . e assim
por diante, todas calculadas em t = T ∗ .
10.17.1 Método 1
Podemos utilizar uma solução bem refinada, por exemplo, uh/16 como sendo
uma boa aproximação da solução exata e supormos que u∗ = uh/16 . Desta forma
podemos aproximar o erro por eh = ku(h) − u∗ k e a ordem do método é estimada
como
log(eh )−log(eh/2 )
p ≈ log(h)−log(h/2)
(10.143)
eh
log
eh/2
≈ log(h/(h/2))
(10.144)
eh
log
eh/2
≈ log(2)
(10.145)
kuh −u∗ k
log
kuh/2 −u∗ k
≈ log(2)
(10.146)
(10.147)
10.17.2 Método 2
Segundo Ferziger/Peric/Roache, podemos também estimar p diretamente de
kuh/2 −uh k
log
kuh/4 −uh/2 k
p ≈ log(2)
(10.148)
(10.149)
Exercícios
y ′′ + y ′ + y = cos(t),
y(0) = 1,
y ′ (0) = 0,
y′ = w
w′ = −w − y + cos(t)
y(0) = 1
w(0) = 0
Exercícios
v ′ = g − αv 2
b) Calcule o valor produzido por cada um desses método para u(1) com passo
h = 0.1, h = 0.05, h = 0.01, h = 0.005 e h = 0.001. Complete a tabela com
os valores para o erro absoluto encontrado.
Problemas de Valores de
Contorno
281
282 Cálculo Numérico
u(a) = ua ⇒ u1 = ua . (11.5)
u(b) = ub ⇒ uN = ub . (11.6)
u1 = ua , (11.7)
1 2 1
ui−1 − 2 ui + 2 ui+1 = −f (xi , ui ), i = 2, . . . , N − 1, (11.8)
h2 h h
uN = ub . (11.9)
N
1 X
E := |u(xi ) − ui | ,
N i=1
xi = (i − 1)h, i = 1, 2, . . . , N,
a = 0
b = 1
N = 11
h = (b-a)/(N-1)
x = np.linspace(a,b,N)
Completamos este sistema com as condições de contorno dadas nas equações (11.11)
e (11.12), donde
u1 = uN = 0.
Ou seja, obtemos o seguinte problema discreto:
u1 = 0, (11.13)
1
− (ui+1 − 2ui + ui+1 ) = 100(xi − 1)2 , i = 2, . . . , N − 1, (11.14)
h2
uN = 0. (11.15)
Observamos que este é um sistema linear N × N , o qual pode ser escrito na forma
matricial Au = b, cujos matriz de coeficientes é
1 0 0 0 0 ··· 0
1 −2 1 0 0 · · · 0
A= 0 1 −2 1 0 ··· 0
,
.. .. .. .. .. .. ..
. . . . . . .
0 0 0 0 0 ··· 1
A[0,0] = 1
b[0] = 0
for i in np.arange(1,N-1):
A[i,i-1] = 1
A[i,i] = -2
A[i,i+1] = 1
b[i] = -100 * h**2 * (x[i]-1)**2
A[N-1,N-1] = 1
b[N-1] = 0
3. Resolução do problema discreto. Neste caso, o problema discreto
consiste no sistema linear Au = b e, portanto, a solução é
u = A−1 b. (11.16)
Figura 11.2: Esboço dos gráficos das soluções analítica (linha) e numérica (pontos)
do PVC dado no exemplo 11.1.1.
(x − 1)3
−uxx = 100(x − 1)2 ⇒ −ux + c1 = 100
3
(x − 1)4
⇒ −u + c2 x + c1 = 100 ,
12
(x − 1)4
ou seja, u(x) = − + c2 x + c1 . As constantes são determinadas pelas
12
condições de contorno dadas pelas equações (11.11) e (11.12), isto é:
100
u(0) = 0 ⇒ c1 = ,
12
100
u(1) = 0 ⇒ c2 = − .
12
Portanto, a solução analítica é:
(x − 1)4 x 100
u(x) = −100 − 100 + (11.17)
12 12 12
Licença CC-BY-SA-3.0. Contato: [email protected]
11.1. MÉTODO DE DIFERENÇAS FINITAS 287
Tabela 11.1: Erro absoluto médio das soluções numéricas com N = 11 e N = 101
do PVC dado no exemplo 11.1.1.
N h E
11 0,1 1,3 × 10−2
101 0,01 1,4 × 10−4
A figura 11.2 mostra o esboço dos gráficos das soluções analítica (11.17) e a da
solução numérica (11.16).
Em Python, podemos fazer o esboço das soluções analítica e numérica da se-
guinte forma:
#grafico
xx = np.linspace(0,1)
yy = np.zeros(50)
for i,xxi in enumerate(xx):
yy[i] = ue(xxi)
plt.plot(x',u,'ro',xx,yy,'b-')
plt.show()
Por fim, computamos o erro absoluto médio das soluções numéricas com N = 11
e N = 101. A tabela 11.1 mostra os resultados obtidos. Observamos, que ao dimi-
nuirmos 10 vezes o tamanho da malha h, o erro absoluto médio diminui aproxima-
damente 100 vezes. Este resultado é esperado, pois o problema discreto (11.13)-
(11.15) aproxima o problema contínuo (11.10)-(11.12) com erro de truncamento
de ordem h2 . Verifique!
Em Python, podemos computar o erro absoluto médio da seguinte forma:
E = 0
for i,xi in enumerate(x):
E += np.abs(ue(xi) - u[i])
E /= N
x u x u
0.50 1.000000 1.00 1.643900
0.60 1.143722 1.10 1.745332
0.70 1.280661 1.20 1.834176
0.80 1.410269 1.30 1.908160
0.90 1.531724 1.40 1.964534
1.00 1.643900 1.50 2.000000
Exercícios Resolvidos
ER 11.1.1. Use o método de diferenças finitas para resolver o seguinte problema
de valor de contorno:
−uxx + u = e−x , 0 < x < 1, (11.18)
u(0,5) = 1, (11.19)
u(1,5) = 2. (11.20)
Para tanto, use a fórmula de diferenças finitas central de ordem 2 para discretizar
a derivada em uma malha uniforme com passo h = 0,1. Faça, então, um esboço
do gráfico da solução computada.
Solução. O passo h é uma malha uniforme com N pontos no domínio [0,5, 1,5]
satisfaz:
(b − a) (b − a)
h= ⇒N = + 1.
N −1 h
Ou seja, a malha deve conter N = 11 pontos igualmente espaçados. Denotamos
os pontos na malha por xi , onde xi = 0,5 + (i − 1)h.
Agora, a equação diferencial dada no i-ésimo ponto da malha é:
−uxx (xi ) + u(xi ) = exi , i = 2, 3, . . . , N − 1.
Denotando ui ≈ u(xi ) e usando a fórmula de diferenças finitas central de ordem
dois para a derivada uxx , obtemos:
ui−1 − 2ui + ui+1
− 2
+ ui = exi ,
h
para i = 2, 3, . . . , N − 1. Rearranjando os termos e aplicando as condições de
contorno, temos o problema discretizado como segue:
u1 = 1
−ui−1 + (2 + h )ui − ui+1 = h2 exi , i = 2, . . . , N − 1,
2
uN = 2.
#malha
a = 0.5
b = 1.5
N = 11
h = (b-a)/(N-1)
x = np.linspace(a,b,N)
#sistema
A = np.zeros((N,N))
b = np.zeros(N)
A[0,0] = 1
b[0] = 1
for i in np.arange(1,N-1):
A[i,i-1] = -1
A[i,i] = 2 + h**2
A[i,i+1] = -1
#solucao
u = np.linalg.solve(A,b)
#grafico
plt.plot(x,u,'b-o')
plt.show()
Exercícios
293
294 Cálculo Numérico
a) No console temos:
>>> quit()
#!/usr/bin/env python
# -*- coding: utf-8 -*-
s = 'Olá'
print(s)
>>> execfile("ola.py")
#!/usr/bin/env python
# -*- coding: utf-8 -*-
s = 'Olá'
print(s)
$ python ola.py
>>> x=1
>>> y = x * 2.0
>>> print(x,y)
(1, 2.0)
>>> type(x), type(y)
(<type 'int'>, <type 'float'>)
and e lógico
or ou lógico
not negação
== igualdade
!= diferente
< menor que
> maior que
<= menor ou igual que
>= maior ou igual que
A.3 Matrizes
Em Python, temos um ótimo suporte para computação científica com o pacote
numpy. Uma matriz A = [ai,j ]m,n
i,j=1 em Python é definida usando-se a seguinte
sintaxe:
>>> print(A)
[[1 2 3]
[4 5 6]]
♦
A seguinte lista contém uma série de funções que geram matrizes particulares:
>>> A = np.ones((3,2))
>>> print(A)
[[ 1. 1.]
[ 1. 1.]
[ 1. 1.]]
>>> nl, nc = np.shape(A)
>>> print(nl,nc)
(3, 2)
A[i,j]
A[i1:i2, j1:j2]
>>> B = np.random.random((4,4))
>>> B
array([[ 0.94313432, 0.72650883, 0.55487089, 0.18753526],
[ 0.02094937, 0.45726099, 0.51925464, 0.8535878 ],
[ 0.75948469, 0.95362926, 0.77942318, 0.06464183],
[ 0.91243198, 0.22775889, 0.04061536, 0.14908227]])
>>> aux = np.copy(B[:,2])
>>> B[:,2] = np.copy(B[:,3])
>>> B[:,3] = np.copy(aux)
>>> B
array([[ 0.94313432, 0.72650883, 0.18753526, 0.55487089],
[ 0.02094937, 0.45726099, 0.8535878 , 0.51925464],
[ 0.75948469, 0.95362926, 0.06464183, 0.77942318],
[ 0.91243198, 0.22775889, 0.14908227, 0.04061536]])
[ 4., 4.]])
>>> A/B
array([[ 0.5, 0.5],
[ 0.5, 0.5]])
#!/usr/bin/env python
# -*- coding: utf-8 -*-
i = 2
if (i == 1):
print("Olá!")
elif (i == 2):
print("Hallo!")
elif (i == 3):
print("Hello!")
else:
print("Ça Va!")
for i in range(6):
print(i)
import numpy as np
for i in np.arange(1,8,2):
print(i)
for k = 10:-3:1
disp(k)
end
import numpy as np
for i in np.arange(10,1,-3):
print(i)
s = 0
i = 1
while (i <= 10):
s = s + i
i = i + 1
A.5 Funções
Além das muitas funções disponíveis em Python (e os tantos muitos pacotes
livers disponíveis), podemos definir nossas próprias funções. Para tanto, existe a
instrução def. Veja os seguintes exemplos:
def f(x):
return x + np.sin(x)
>>> f(np.pi)
def h(x,y):
if (x < y):
return y - x
else:
return x - y
define a função:
y − x ,x < y
h(x,y) =
x − y ,x ≥ y
def J(x):
y = np.zeros((2,2))
y[0,0] = 2*x[0]
y[0,1] = 2*x[1]
y[1,0] = -x[1]*np.sin(x[0]*x[1])
y[1,1] = -x[0]*np.sin(x[0]*x[1])
return y
∂(f1 ,f2 )
define a matriz jacobiana J(x1 ,x2 ) := ∂(x1 ,x2 )
da função:
A.6 Gráficos
Para criar um esboço do gráfico de uma função de uma variável real y = f (x),
podemos usar a biblioteca Python mathplotlib.A função matplotlib.pyplot.plot
faz uma representação gráfica de um conjunto de pontos {(xi , yi )} fornecidos.
Existe uma série de opções para esta função de forma que o usuário pode ajustar
várias questões de visualização. Veja a documentação.
E 2.1.3. (101,1)2 .
E 2.1.4. (11,1C)16 .
¯ 5 ; b) (45,1)6 .
E 2.1.5. a) (12,31)
E 2.2.1.
a) 2,99792458 × 105 b) 6,62607 × 10−34
c) 6,674 × 10 −8
d) 9,80665 × 104
305
306 Cálculo Numérico
E 2.5.1. a) εabs = 5,9 × 10−4 , εrel = 1,9 × 10−2 %; b) εabs = ×10−5 , εrel = ×10−3 %; c) εabs = 1, εrel = 10−5 %.
E 2.5.4. a) 2; b) 2.
E 2.5.6. a) δabs = 3,46 × 10−7 , δrel = 1,10 × 10−7 ; b) δabs = 1,43 × 10−4 , δrel = 1,00 × 10−3 .
E 2.8.1. 2%, deve-se melhorar a medida na variável x, pois, por mais que o erro relativo seja maior para esta variável, a
propagação de erros através desta variáveis é muito menos importante do que para a outra variável.
E 2.8.2. 3,2% pela aproximação ou 3,4% pela segundo método (0,96758 ≤ I ≤ 1,0342).
E 2.9.1. Quando µ é pequeno, e1/µ é um número grande. A primeira expressão produz um ”overflow” (número maior que o
máximo representável) quando µ é pequeno. A segunda expressão, no entanto, reproduz o limite 1 quando µ → 0+.
2 √ √
2 2
E 2.9.2. a) 1
2
+ x4! + O(x4 ); b) x/2 + O(x2 ); c) 5 · 10−4 x + O(x2 ); d) 4
y + O(y 2 ) = 4
x + O(x2 )
E 2.9.5. 4,12451228 × 10−16 J; 0,002%; 0,26654956 × 10−14 J; 0,002%; 4,98497440 × 10−13 J; 0,057%; 1,74927914 × 10−12 J;
0,522%.
E 3.1.1.
Observamos que a equação é equivalente a cos(x) − x = 0. Tomando, então, f (x) = cos(x) − x, temos que f (x) é contínua
em [0, π/2], f (0) = 1 e f (π/2) = −π/2 < 0. Logo, do teorema de Bolzano 3.1.1, concluímos que a equação dada tem pelo menos
uma solução no intervalo (0, π/2).
E 3.1.2.
No exercício 3.1.1, mostramos que a função f (x) = cos(x) − x tem um zero no intervalo [0, π/2]. Agora, observamos que
f ′ (x) = − sen (x)−1. Como 0 < sen x < 1 para todo x ∈ (0, π/2), temos que f ′ (x) < 0 em (0, π/2), isto é, f (x) é monotonicamente
decrescente neste intervalo. Logo, da proposição 3.1.1, temos que existe um único zero da função neste intervalo.
E 3.1.3.
k ≈ 0,161228
E 3.1.5.
Escolhendo o intervalo [a, b] = [−1,841 − 10−3 , −1,841 + 10−3 ], temos f (a) ≈ 5 × 10−4 > 0 e f (b) ≈ −1,2 × 10−3 < 0, isto
é, f (a) · f (b) < 0. Então, o teorema de Bolzano nos garante que o zero exato x∗ de f (x) está no intervalo (a, b). Logo, da escolha
feita, | − 1,841 − x∗ | < 10−3 .
E 3.2.1. 0,6875
E 3.2.2. Intervalo (0,4, 0,5), zero 0,45931. Intervalo (1,7, 1,8), zero 1,7036. Intervalo (2,5, 2,6), zero 2,5582.
E 3.2.6. kθ = lP
2
cos(θ) com θ ∈ (0, π/2); 1,030.
E 3.3.1. −1,8414057
E 3.3.2.
0,7391
E 3.3.3.
Tomemos x(1) = 1 como aproximação inicial para a solução deste problema, iterando a primeira sequência a), obtemos:
(1)
x = 1
(2) 10
x = ln = 2,3025851
1
(3) 10
x = ln = 1,4685526
2,3025851
.
.
.
(21)
x = 1,7455151
(31)
x = 1,745528
(32)
x = 1,745528
(1)
x = 1
(2) −1
x = 10e = 3,6787944
(3) −3,6787944
x = 10e = 0,2525340
(4) −0,2525340
x = 10e = 7,7682979
(5) −7,7682979
x = 10e = 0,0042293
(6) −0,0042293
x = 10e = 9,9577961
Este experimento numérico sugere que a iteração a) converge para 1,745528 e a iteração b) não é convergente.
E 3.3.10.
0.0431266
cos(x)−x2
E 3.4.1. raiz:0,82413, processo iterativo: x(n+1) = x(n) + sen (x)+2x
E 3.4.3. 0,65291864
E 3.4.9.
x0 > 1.
E 3.4.10.
(0)
x = C.I.
(n+1) (n) (n)
x = x 2 − Ax
E 3.4.11.
x0 = C.I.
(n+1) (n) 1 A
x = x 1− +
n nx(n)
E 3.4.12.
x0 = C.I.
E 3.6.5. Seja f (x) ∈ C 2 um função tal que f (x∗ ) = 0 e f ′ (x∗ ) 6= 0. Considere o processo iterativo do método das secantes:
Definimos ǫn = xn − x∗ , equivalente a xn = x∗ + ǫn
∗ ∗ ′ ∗ 2 f ′′ (x∗ )
f (x + ǫ) ≈ f (x ) + ǫf (x ) + ǫ
2
∗ ′ ∗ 2 f ′′ (x∗ )
f (x + ǫ) ≈ ǫf (x ) + ǫ
2
f ′′ (x∗ )
f ′′ (x∗ )
ǫn ǫn−1 f ′ (x∗ ) + ǫ2
n−1 2
− ǫn−1 ǫn f ′ (x∗ ) + ǫ2
n 2
ǫn+1 ≈
f (x∗ + ǫn ) − f (x∗ + ǫn−1 )
f ′′ (x∗ )
2
ǫn ǫ2 2
n−1 − ǫn−1 ǫn
=
f (x∗ + ǫn ) − f (x∗ + ǫn−1 )
1 ′′ ∗ ǫn ǫn−1 ǫn−1 − ǫn
= f (x )
2 f (x∗ + ǫn ) − f (x∗ + ǫn−1 )
∗ ∗
∗ ′ ∗
∗ ′ ∗
f (x + ǫn ) − f (x + ǫn−1 ) ≈ f (x ) + f (x )ǫn − f (x ) + f (x )ǫn−1
(3.6)
′ ∗
= f (x )(ǫn − ǫn−1 )
Portanto:
1 f ′′ (x∗ )
ǫn+1 ≈ ǫn ǫn−1 (3.7)
2 f ′ (x∗ )
ou, equivalentemente:
(n+1) ∗ 1 f ′′ (x∗ ) (n) ∗
(n−1) ∗
x −x ≈ x −x x −x (3.8)
2 f ′ (x∗ )
E 3.7.2.
x > a com a ≈ 0,4193648.
E 3.7.3.
−2j+1
z1 ≈ 0.3252768, z2 ≈ 1.5153738, z3 ≈ 2.497846, z4 ≈ 3.5002901, zj ≈ j − 1/2 − (−1)j e π , j >4
E 3.7.4.
150 W, 133 W, 87 W, 55 W, 6,5 W
E 3.7.5.
a) 42 s e 8 min2 s, b) 14 min56 s.
E 3.7.6.
118940992
E 3.7.7.
7,7 cm
E 3.7.8.
4,32 cm
E 3.7.9.
(0,652919, 0,426303)
E 3.7.10.
7,19% ao mês
E 3.7.11.
4,54% ao mês.
E 3.7.12.
500 K, 700 K em t = 3 ln(2), 26 min, 4 h27 min.
E 3.7.13.
(±1,1101388, −0,7675919), (±1,5602111, 0,342585)
E 3.7.14.
1,5318075
E 3.7.15.
Aproximadamente 2500 reais por hora.
E 3.7.16.
a) 332,74 K b) 359,33 K
E 3.7.17.
1,2285751, 4,76770758, 7,88704085
E 4.1.1.
Escrevemos o sistema na forma matricial e resolvemos:
Portanto x = 2, y = 3, z = −5
E 4.5.1.
a = (0, 1, 2, 1, 2, 1)
b = (5, 3, 4, 2, 3, 2)
c = (4, 1, 1, 1, 2, 0)
d = (13, 10, 20, 16, 35, 17)
x = (1, 2, 3, 4, 5, 6)
E 4.5.2.
a = (0, −1, −1, −1, −1, −1, −1, −1, −1, −1, −1/2)
b = (1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 1)
c = (−1, −1, −1, −1, −1, −1, −1, −1, −1, −1, 0)
d = (0, cos(2/10), cos(3/10), cos(4/10), cos(5/10),
cos(6/10), cos(7/10), cos(8/10), cos(9/10), cos(1),0)
x = (0,324295, 0,324295, 0,317115, 0,305943, 0,291539,
0,274169, 0,253971, 0,230846, 0,20355, 0,165301, 0,082650)
E 4.6.1.
λ = 71×30
41
≈ 51.95122, para λ = 51: k1 = k∞ = 350.4, k2 = 262.1. Para λ = 52: k1 = k∞ = 6888, k2 = 5163.
E 4.6.2.
k1 (A) = 36, k2 (A) = 18,26, K∞ (A) = 20,8.
E 4.6.3.
√
k1 = k∞ = 6888, k2 = 26656567 e k1 = 180, k2 = 128,40972 e k∞ = 210
E 4.6.4.
18 + 3. Quando ε → 0+, a matriz converge para uma matriz singular e o número de condicionamento diverge para +∞.
ε
E 4.6.5.
As soluções são [−0.0000990 0.0000098]T e [0.0098029 0.0990294]T . A grande variação na solução em função de pequena
variação nos dados é devido ao mau condicionamento da matriz (k1 ≈ 1186274.3).
Exemplo de implementação:
E 4.7.4.
0,324295, 0,324295, 0,317115, 0,305943, 0,291539, 0,274169, 0,253971, 0,230846, 0,203551, 0,165301, 0,082650
E 4.7.5.
Permute as linhas 1 e 2.
E 4.8.1.
λ = 86.1785 associado ao autovetor dado por v1 = [0.65968 0.66834 0.34372]T .
E 4.8.3.
158,726
E 4.9.3.
Dica: P (−1) = −3, P (1) = −1 e P (2) = 9 produzem três equações lineares para os coeficientes a, b e c. Resp: a)
P (x) = 3x2 + x − 5, b) A ≈ 2.49 e B ≈ −1.29 c)A1 ≈ 1.2872058, A2 ≈ −4.3033034, B1 ≈ 2.051533 e B2 ≈ −0.9046921.
cos(x) − x sen (x) 1
JF =
−2e−2x+y e−2x+y
p
√ √
E 5.1.3. As curvas possuem dois pontos de intersecção. A posição exata destes pontos de intesecção é dada por 2 3 − 3,2 3 − 2
p
√ √
e − 2 3 − 3,2 3 − 2 . Use a solução exata para comparar com a solução aproximada obtida.
Da condição de tangência, temos que o coeficiente angular da reta, m, deve igual à derivada da função f (x) nos dois pontos
de tangência.
′ ′
m = f (x1 ) = f (x2 )
E sabemos que:
′ cos(x) sen (x)
f (x) = − .
1+x (1 + x)2
Assim, podemos reescrever o problema como
E 5.1.11. A primeira curva trata-se de uma elipse de centro (3,1) e semi-eixos 4 e 6, portanto seus pontos estão contidos no
retângulo −1 ≤ x ≤ 7 e −5 ≤ y ≤ 7.
A soluções são (−0,5384844, −1,7978634) e (2,8441544, 6,9954443).
p 6,7
E 5.1.13. Incialização do método: A(0) = 3,1 e b(0) = 3,1
A ≈ 3.0297384 e b ≈ 1.4835346.
E 5.1.17.
x1 − x2
−x1 + 5(x2 + x3
2 ) − x3 − 10 exp(−2/3)
−x2 + 5(x3 + x3
3 ) − x4 − 10 exp(−3/3)
F (x) = −x3 + 5(x4 + x3
4 ) − x5 − 10 exp(−4/3)
.
.
.
−x9 + 5(x10 + x3
10 ) − x11 − 10 exp(−10/3)
x11 − 1
1 −1 0 0 0 ... 0
−1 5(1 + 3x2
2) −1 0 0 ... 0
0 −1 5(1 + 3x2
3) −1 0 ... 0
JF (x) =
0 0 −1 5(1 + 3x2
4) −1 ... 0
. . . . .
.
. . . . .. .
. . . . .
0 0 0 0 0 ··· 1
Resposta final: 0,80447, 0,80447, 0,68686, 0,57124, 0,46535, 0,37061, 0,28883, 0,22433, 0,19443, 0,28667, 1
E 6.1.3.
2
a1 + a2 x1 + a3 x1 = y1
2
a1 + a2 x2 + a3 x2 = y2
2
a1 + a2 x1 = y1 − a3 x1
,
2
a1 + a2 x2 = y2 − a3 x2
o qual tem solução única, pois x1 6= x2 . Ou seja, para cada a3 ∈ R dado, existem a1 , a2 ∈ R tais que a parábola de
equação y = a1 + a2 x + a3 x2 interpola os pontos dados.
b) Certamente não existem retas de equação x = a que interpolam os pontos dados. Consideremos então retas de equação
y = a1 + a2 x. Para uma tal reta interpolar os pontos dados é necessário que:
a1 + a2 = 1
a1 + 2a2 = 2,1,
a1 + 3a2 = 3
c) Não existe uma parábola de equação y = a1 + a2 x + a3 x2 que interpole os pontos dados, pois tal equação determina uma
função de x em y. Agora, para mostrar que existem infinitas parábolas de equação x = a1 + a2 y + a3 y 2 que interpolam
os pontos dados, basta seguir um raciocínio análogo ao do item a), trocando x por y e y por x.
E 6.4.1.
R1 f (0)+f (1) 1
P (x)dx = 2
, 12 maxx∈[0,1] |f ′′ (x)|
0
E 8.1.3.
1
h h2 h
h
c) f ′ (0) = h +h l − h2 f (−h1 ) + h1
− h1 f (0) + h1 f (h2 )
1 2 1 2 2
E 8.1.4.
Caso a b c d
vi = 1 1.72 1.56 1.64 1.86
vi = 4.5 2.46 1.90 2.18 1.14
E 8.1.5.
Segue a tabela com os valores da derivada para vários valores de h.
Observe que o valor exato é −0,3161977 e o h ótimo é algo entre 10−8 e 10−9 .
E 8.2.1.
E 9.1.2.
1 2
ISimpson = IT rap + IP M
3 3
E 9.1.3.
E 9.3.2.
-0.2310491, -0.2452073, - 0.2478649.
E 9.3.4.
a)-0.2472261, -0.2416451, -0.2404596, -0.2400968, -0.2399563, -0.2398928. b)-0.2393727, -0.2397994, -0.2398104, -0.2398115,
-0.2398117, -0.2398117.
E 9.4.1.
−1 −12 −2 2 −7 3 −3 4
a)I(h) = 4.41041 · 10 − 8.49372 · 10 h − 1.22104 · 10 h − 1.22376 · 10 h + 8.14294 · 10 h
−1 −11 −2 2 −7 3 −6 4
b)I(h) = 7.85398 · 10 − 1.46294 · 10 h − 4.16667 · 10 h − 2.16110 · 10 h + 4.65117 · 10 h
−3 −10 −7 2 −5 3 −4 4
c)I(h) = 1.58730 · 10 − 9.68958 · 10 h + 2.03315 · 10 h − 1.38695 · 10 h + 2.97262 · 10 h
−1 −12 −2 2 −8 3 −4 4
d)I(h) = 4.61917 · 10 + 3.83229 · 10 h + 2.52721 · 10 h + 5.48935 · 10 h + 5.25326 · 10 h
E 9.4.2.
1.5707963 2.0943951
1.8961189 2.0045598 1.9985707
1.9742316 2.0002692 1.9999831 2.0000055
E 9.4.5. R(6,6) = −10.772065, R(7,7) = 5.2677002, R(8,8) = 6.1884951, R(9,9) = 6.0554327, R(10,10) = 6.0574643. O
valor desta integral com oito dígitos corretos é aproximado por 6.0574613.
E 9.5.1.
w1 = 1/6, w2 = 2/3, w3 = 1/6. O esquema construído é o de Simpson e a ordem de exatidão é 3.
E 9.5.2. 3
E 9.5.3. 5
E 9.5.5. 5, 4, 3
E 9.6.1.
E 9.7.1.
n b c d e f
2 2.205508 3.5733599 3.6191866 3.6185185 3.618146
4 2.5973554 3.6107456 3.6181465 3.6180970 3.6180970
6 2.7732372 3.6153069 3.6181044 3.6180970 3.6180970
8 2.880694 3.6166953 3.6180989 3.6180970 3.6180970
temos
∞
1 − cos(x)
X x2n−1/2
n
√ =− (−1) , x≥0
x (2n)!
n=1
Z 1 ∞
X Z 1
cos(x) − 1 n x2n−1/2
I = 4+2 p dx = 4 − 2 (−1) dx
(2n)!
0 |x| 0
n=1
∞
X 1
n
= 4−2 (−1)
(2n)!(2n + 1/2)
n=1
Solução do item f)
Z 1
−1/2 x3/2 x7/2 1 1 977
2 x − + dx = 2 2− + =
2 24 5 54 270
0
Z 1 Z 1 1+u
1+u
cos(x) − P4 (x) √ cos 2
− P4 2
2 √ dx = 2 √ du
x 1+u
0 −1
E 9.7.5. 4,1138
Z 1 √ √
3 3
E 9.7.8. f (x)dx = f − +f
3 3
−1
E 10.2.1.
2
−1 −2 1+e−1
0,4496 com h = 0,1 e 0,4660 com h = 0,01. A solução exata vale u(1) = 1+2e 4 +e = 2
≈ 0,4678
E 10.2.2.
u(2) ≈ 0,430202 e z(2) = 0,617294 com h = 0,2, u(2) ≈ 0,435506 e z(2) = 0,645776 com h = 0,02, u(2) ≈ 0,435805 e
z(2) = 0,648638 com h = 0,002 e u(2) ≈ 0,435832 e z(2) = 0,648925 com h = 0,0002.
E 10.2.3.
u(2) ≈ 1,161793 com h = 0,1, u(2) ≈ 1,139573 com h = 0,01, u(2) ≈ 1,137448 com h = 0,001, u(2) ≈ 1,137237 com
E 10.6.1.
u(1) ≈ 1,317078 quando h = 0,1 e u(1) ≈ 1,317045.
E 10.6.2.
E 10.19.1.
1 ln 9
1 ln (6)
Os valores exatos para os itens e e f são: 10 4
e 10
E 10.19.2.
q p
g
O valor exato é α
1 − e−200α em t = √1
gα
tanh−1 1 − e−200α
E 10.19.8.
E 11.1.1.
1 0 0 0 0
u1
5
−1 2 −1 0 0 u2 2
0 −1 2 −1 0
u3
= 2
0 0 −1 2 −1 u4 2
0 0 0 0 1 u5 10
Solução: [5, 7.375, 9.25, 10.625, 11.5, 11.875, 11.75, 1.125, 10]
E 11.1.2. 120. 133.56 146.22 157.83 168.22 177.21 184.65 190.38 194.28 196.26 196.26 194.26 190.28 184.38 176.65 167.21
E 11.1.3. 391.13 391.13 390.24 388.29 385.12 380.56 374.44 366.61 356.95 345.38 331.82 316.27 298.73 279.27 257.99 234.99
E 11.1.4. 0., 6.57, 12.14, 16.73, 20.4, 23.24, 25.38, 26.93 , 28, 28.7, 29.06, 29.15, 28.95, 28.46, 27.62 , 26.36, 24.59, 22.18,
[3] R.L. Burden and J.D. Faires. Análise Numérica. Cengage Learning, 8 edition,
2013.
[6] Walter Gautschi and Gabriele Inglese. Lower bounds for the condition
number of vandermonde matrices. Numerische Mathematik, 52(3):241–250,
1987/1988.
[8] E. Isaacson and H.B. Keller. Analysis of numerical methods. Dover, Ontário,
1994.
[9] Arieh Iserles. A first course in the numerical analysis of differential equations.
Cambridge university press, 2009.
[10] W.H. Press. Numerical Recipes 3rd Edition: The Art of Scientific Computing.
Cambridge University Press, 2007.
318
REFERÊNCIAS BIBLIOGRÁFICAS 319
Aqui você encontra a lista de colaboradores do livro. Esta lista contém so-
mente aqueles que explicitamente se manifestaram a favor de terem seus nomes
registrados aqui. A lista completa de colaborações pode ser obtida no repositório
GitHub do livro:
https://github.com/livroscolaborativos/CalculoNumerico
Além das colaborações via GitHub, o livro também recebe colaborações via dis-
cussões, sugestões e avisos deixados em nossa lista de e-mails:
Estas colaborações não estão listadas aqui, mas podem ser vistas no site do grupo
de e-mails.
Caso encontre algum equívoco ou veja seu nome listado aqui por engano, por
favor, entre em contato conosco por e-mail:
320
Índice Remissivo
321
322 Cálculo Numérico
sequência de
Fibonacci, 80
simulação
computacional, 1
numérica, 1
sistema de equações
não lineares, 138
sistema de numeração, 3
sistema linear, 87
condicionamento, 111
sistema numérico
de ponto fixo, 16
de ponto flutuante, 17
notação normalizada, 10
sistemas
de equações diferenciais, 275