book CN

Prévia do material em texto

Em
El
ab
or
aç
ão
Notas de Aula de
Cálculo Numérico
Edson Leite Araújo
Lino Marcos da Silva
Sergio Floquet Sales
Em
El
ab
or
aç
ão
Sumário
1 Sistemas Numéricos 1
1.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Sistemas Binário e Decimal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3 Representação de Inteiros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.4 Conversão Binário para Decimal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.5 Conversão Decimal para Binário . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.6 ExercÃcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.7 Representação de Frações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.8 ExercÃcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.9 Aritmética de Ponto Flutuante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.10 Instabilidades Numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.11 ExercÃcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2 Zeros Reais de Funções Reais 17
2.1 Fase I: Isolamento das Raízes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Fase II: Refinamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3 Método da Bissecção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4 Estimativa do Número de Iterações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.5 Método da Posição Falsa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.6 Método do Ponto Fixo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.7 Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.8 Método da Secante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.9 Métodos Estudados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3 Sistemas de Equações Lineares 28
3.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Sistemas Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3 Resolução de Sistemas Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.4 Métodos Diretos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4 Introdução à Resolução de Sistemas Não Lineares 35
4.1 Método de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.2 Método de Newton Modificado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5 Interpolação 40
5.1 Interpolação Polinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.2 Formas para se obter o polinômio interpolador Pn(x) . . . . . . . . . . . . . . . . . . . . . 41
ii
Em
El
ab
or
aç
ão
Sumário iii
5.3 Resolução do Sistema Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.4 Forma de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.5 Forma de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.6 Forma de Newton para pn(x) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.7 Estudo do Erro na Interpolação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.8 Interpolação Inversa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.9 Fenômeno de Runge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.10 Spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.11 Spline Linear Interpolante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.12 Spline Cúbica Interpolante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6 Ajuste de Curvas pelo Método dos Mínimos Quadrados 61
6.1 Método dos Mínimos Quadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
7 Soluções Numéricas para Equações Diferenciais Ordinárias 66
7.1 Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
7.2 Método de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
7.3 Método das Diferenças Finitas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
Índice Remissivo 75
Em
El
ab
or
aç
ão
Em
El
ab
or
aç
ão
1
Sistemas Numéricos
1.1 Introdução
Nos computadores e em máquinas de calcular, em geral os números são representados de maneira limitada.
Dificuldades relacionadas ao armazenamento em memória, velocidade de processamento e padronização,
fazem com que estudos detalhados sobre a representação dos números em ambientes computacionais,
sejam necessários.
Neste capítulo, pretendemos entender de que maneira esta representação é possível, qual a melhor em
cada situação e quais as limitações de cada uma deles.
1.2 Sistemas Binário e Decimal
De um modo geral, os números inteiros são representados no sistema decimal. Por exemplo, o número
428
na verdade, significa
428 = 4 · 102 + 2 · 101 + 8 · 100
Cada digito possui um valor valor posicional e é determinado pelo seu valor absoluto multiplicado por
uma potência de 10, escolhida de acordo com sua posição. Assim, o 4 em 428 está na posição 2, contando
a partir do algarismo das unidades e tendo o 0 como valor inicial, e portanto seu valor é 4 · 102 = 400.
Tudo indica que a escolha da base 10 como sistema mais usado pela humanidade tenha se dado pelo
simples fato de termos dez dedos e isto induzir a base 10 como forma natural de se realizar os cálculos
em geral.
Apesar de toda nossa familiaridade com o sistema decimal, este não se adequa aos limites e imposição
feitas pelos mecanismos e circuitos que compõem o interior de um computador. Vejamos porquê.
As máquinas, em geral, entendem os números através de pulsos de corrente elétrica. Tem-se o 1
quando há e o 0 quando não há corrente elétrica. De forma que o sistema binário é o sistema numérico
natural para se representar os números em máquinas. Uma vez que tal limitação não surgiu apenas com
os números inteiros, tal sistema é usado, na verdade, para se tentar representar todos os números reais.
1.3 Representação de Inteiros
Um dado inteiro N no sistema binário é representado da seguinte forma:
N = an2n + an2n−1 + ... + a222 + a121 + a020
Em
El
ab
or
aç
ão
2 1. Sistemas Numéricos
onde os ai são zeros ou uns, e sua escrita na base 2 é dada por:
anan−1...a1a0
Usaremos o sublinhado para distinguir a representação binária da decimal.
Exemplo 1.3.1
1 = 1
2 = 2
3 = 101
100 = 1100100
1.4 Conversão Binário para Decimal
Como nós, humanos, estamos habituados a pensar os números no sistema decimal e as máquinas os
entendem, de um modo geral, no sistema binário, é necessário que exista uma forma de conversão entre
estes sistemas, caso queiramos nos comunicar com eles.
Considere N sendo um número inteiro não negativo, cuja representação no sistema binário é dada por
N = ana1...a0,
onde ai = 1 ou 0, para i = 0, ...,n, ou seja
N = an2n + an−12n−1 + ... + a121 + a0
Observe que
N = p(2)
sendo
p(x) = anx
n + an−1x
n−1 + ... + a1x+ a0
perceba ainda que o polinômio p pode ser reescrito da seguinte maneira
p(x) =
(
anx
n−1 + an−1x
n−2 + ... + a1
)
x+ a0
=
((
anx
n−2 + an−1x
n−3 + ... + a2
)
x+ a1
)
x+ a0
=



 an︸︷︷︸
b0
x+ an−1
︸ ︷︷ ︸
b1
 x+ an−2

︸ ︷︷ ︸
b2
x+ ... + a2

x+ a1
︸ ︷︷ ︸
bn−1
x+ a0
Esta forma de enxergar o polinômio p, torna possível colocar na forma de algoritmo, o cálculo de p(2).
Em
El
ab
or
aç
ão
1.5. Conversão Decimal para Binário 3
x = 2
b0 = an
b1 = b0 · x+ an−1
b2 = b1 · x+ an−2
...
bn = bn−1 · x+ a0
Algoritmo 1: Binário → Decimal
Entrada: N = a1an−1 · · ·a0
início
x = 2; b0 = an;
para i = 1; i 6 n; i++ faça
bi = bi−1 · x+ an−i;
fim
retorne bn;
fim
1.5 Conversão Decimal para Binário
Para a conversão de inteiros do sistema decimal para o binário convém observarmos que, o último dígito
de um dado inteiro N é zero, se e somente se, N é par. O dígito seguinte, do mesmo modo, será zero
se, e somente se, N−a0
2 é par, e assim sucessivamente.
Desta forma, um algoritmo possível para a conversão de decimal para binário seria:
Algoritmo 2: Decimal → Binário
Entrada: N (na forma decimal)
Saída: a; Obs: a indexação de a não é a mesma utilizada em C/C#;
início
n = N; i = 0;
enquanto n > 1 faça
se n%2 == 0 então
ai = 0;
fim
senão
ai = 1;
fim
n =
(n− ai)
2
;
fim
retorne a;
fim
O algoritmo anterior pode ser melhor compreendido através de um esquema prático que veremos a
seguir, na forma de um exemplo:
Exemplo 1.5.1
Em
El
ab
or
aç
ão
4 1. Sistemas Numéricos
N = 50
50 2
0 25 2
↑ 1 12 2
a5 ↑ 0 6 2
a4 ↑ 0 3 2
a3 ↑ 1 1
a2 ↑ ↑
a1 a0
N = a0a1a2a3a4a5 = 110010
1.6 ExercÃcios
1.6.1 ???
1.7 Representação de Frações
Uma fração binária (aqui estaremos pensando apenas em números menores que 1) é um número que
pode ser escrito na forma
0 · a−1a−2a−3 · · ·a−n = a−12−1 + a−22−2 + ... + a−n2−n
=
n∑
i=1
a−12−i
onde os coeficientes ai ∈ {0, 1}.
Observação 1.7.1 Um fato interessante é que uma dada fração pode ter duas representações
diferentes. Por exemplo
0.5 =
1
2
= 2−1 = 0.1
Por outro lado, temos também
0.5 = 0.0111111...
Em
El
ab
or
aç
ão
1.7. Representação de Frações 5
Uma vez que:
0.0111111... = 2−2 + 2−3 + 2−4 + ... + 2−n + ...
=
1
4(
1 − 1
2
) (Usando a soma de uma P.G)
=
1
4
1
2
=
1
2
= 0.5
Tal característica não é exclusiva da representação binária uma vez que, no sistema decimal, por
exemplo:
1 = 0.9999...
Entre as duas possibilidades, em ambos os sistemas, opta-se pela mais simples para representar uma
dada fração. �
Uma fração binária, dada por
z = 0.a−1a−2 · · ·a−n
pode ser vista como o valor do polinômio
p(x) = a−1x+ a−2x
2 + · · ·+ a−nxn
= a−nx
n + a−n+1x
n−1 + · · ·+ a−2x
2 + a−1x
em x = 1
2 . Assim, procedendo-se de maneira semelhante à dedução do Algoritmo 1, teremos o seguinte
algoritmo:
Algoritmo 3: Binário → Decimal (Frações)
Entrada: N = 0 · a−1a−2 · · ·a−n
início
x = 1/2; b0 = a−n; a0 = 0
para i = 1; i 6 n; i++ faça
bi = bi−1x+ a−n+i;
fim
retorne bn;
fim
Exemplo 1.7.1
???
Em
El
ab
or
aç
ão
6 1. Sistemas Numéricos
Observação 1.7.2 Sendo z = 0.a−1a−2 · · ·a−n uma fração binária, outro modo de convertê-la para
uma fração decimal, consiste em converter o inteiro
2n. z = a−1a−2 · · ·a−n
para decimal, e dividir o resultado por 2n. �
Precisamos agora, pensar em como seria a conversão de uma fração decimal para fração binária.
Para isto, considere z ∈ R, tal que 0 6 z 6 1 e suponha que após convertido, sua representação no sistema
binário seja:
z = 0.a−1a−2a−3 · · ·a−n
tal representação pode também ser vista como:
z = a−12−1 + a−22−2 + a−32−3 + · · ·+ a−n2−n
Observe que
2z = 2
(
a−12−1 + a−22−2 + a−32−3 + ... + a−n2−n
)
= a−120 + a−22−1 + a−32−2 + ... + a−n2−n+1
= a−1.a−2a−3 ...a−n
Ou seja, ao multiplicarmos o número z por 2, fazemos com que sua representação binária se desloque
uma posição para a esquerda. Com isto podemos concluir que:
2z > 1⇒ a−1 = 1
2z 6 1⇒ a−1 = 0
Repetindo este procedimento consecutivamente, tomando o cuidado de remover, a cada repetição, a
parte inteira de z, teremos então um algoritmo para a conversão do número z numa fração binária.
Algoritmo 4: Decimal → Binário (Frações)
Entrada: z ∈ R, 0 6 z < 1
início
α = z; i = 1
enquanto α > 0 faça
α = 2α;
se α > 1 então
a−i = 1;
fim
senão
a−i = 0;
fim
α = α− parteInteira(α);
fim
retorne a[ ];
fim
Em
El
ab
or
aç
ão
1.8. ExercÃcios 7
Exemplo 1.7.2
Converter 0.7 para binário.
Na Tabela 1.1, temos em cada linha, o que seria a execu-
ção de cada iteração do laço existente no algoritmo anterior.
Perceba, que a partir da sexta linha, teremos a repetição das
linhas 2, 3, 4 e 5, ou seja, temos uma dízima periódica.
0.7 = 0 · 101100 01100 · · ·
i α 2α a−i
1 0.7 1.4 1
2 0.4 0.8 0
3 0.8 1.6 1
4 0.6 1.2 1
5 0.2 0.4 0
6 0.4 0.8 0
...
...
...
...
Tabela 1.1: Decimal →
Binário
Observe que o Algoritmo 4 possui uma falha, uma vez que em geral as frações decimais, quando
escritas na forma binária possui uma quantidade infinita de dígitos após a vírgula. São dízimas periódicas
em sua grande maioria.
Há portanto a necessidade, por questões de armazenamento em memória, de se limitar (truncar) esta
representação, determinando para isto um número máximo de dígitos após a vírgula.
Concluimos esta seção observando que, em posse dos quatro algoritmos apresentados, estamos aptos
à conversão de qualquer número real do sistema binário para o decimal e vice-versa. Bastando para isso
apenas uma decomposição deste número em suas partes inteiras e fracionárias.
1.8 ExercÃcios
1.8.1 ???
1.9 Aritmética de Ponto Flutuante
Os números reais, apesar do nome, apresentam-se impossíveis de serem representados no mundo real
dos computadores. Entretanto, muitos métodos têm sido propostos para aproximar o sistema numérico
real usando representações finitas. O método utilizado em quase todos os computadores, hoje em dia, é
conhecido como Aritmética de Ponto Flutuante.
Nesta representação, um número x ∈ R, possui três partes:
- sinal
- mantissa
- expoente
e escreve-se
x = ±︸︷︷︸
sinal
.d1d2 · · ·dt︸ ︷︷ ︸
mantissa
· β
expoente︷︸︸︷
e
onde
• β: base do sistema (2 em geral)
• t: número de dígitos na mantissa (precisão)
• e: expoente
• d′is: dígito da mantissa, sendo d1 6= 0 e 0 6 di 6 β− 1, i = 2, ..., t.
Em
El
ab
or
aç
ão
8 1. Sistemas Numéricos
Por razões de precisão, que ficarão claras mais adiante, é comum que a representação em ponto
flutuante seja normalizada, ou seja, a mantissa é maior ou igual a 1.
Observação 1.9.1 Numa representação normalizada, atenção especial deve ser dada ao número
zero. Em geral usa-se
di = 0, i = 1, ..., t
e = 0 (menor expoente não negativo)
Uma consequência interessante da normalização, para sistemas de base 2, é que o primeiro dígito da
mantissa é sempre "1" e, assim sendo, não há necessidade de armazená-lo. Em muitos sistemas computa-
cionais, isto é o que, de fato, é feito. Este bit1 é conhecido como bit oculto ou bit implícito.
Quantidade de Números
Considere um sistema (uma aritmética) de ponto flutuante na base β com t dígitos na mantissa e expoente
e ∈ [emin, emax]. De agora em diante denotaremos um tal sistema por F(β, t, emin, emax). Um número
x ∈ R neste sistema, supondo representação normalizada terá o seguinte formato:
x = ± d1.d2d3 ...dt · βe
Para sabermos quantos números podem ser representados neste sistema, observe que para o dígito d1
existem β − 1 possibilidades2, para os t − 1 dígitos restantes, tem-se β possibilidades para cada um. Ou
seja, com t dígitos disponíveis na mantissa, podemos representar (β − 1)βt−1 números. Levando-se em
conta o sinal, que pode ser “+” ou “−”, esta quantidade passa a ser
2(β− 1)βt−1
Para o expoente e ∈ [emin, emax] teremos
emax − emin + 1
possíveis números. Portanto, a quantidade de números representáveis numa aritmética F(β, t, emin, emax)
incluindo a representação especial para o número zero, é
2(β− 1)βt−1(emax − emin + 1) + 1
O que nos leva à conclusão já esperada, que o conjunto F é um subconjunto finito dos números reais,
tendo o valor β1−t como menor distância possível entre os números representáveis na mantissa.
Exemplo 1.9.1
Consideremos o sistema F(2, 3,−1, 3). Neste sistema temos β = 2, t = 3, emin = −1, e emax = 3.
Logo, supondo representação normalzada, a mantissa deste sistema terá a forma
d1 ·d2d3 d1 6= 0, d2,d3 ∈ {0, 1}
ou seja, temos as seguintes possibilidades:
1A palavra bit é uma abreviação de binay digit, que significa dígito binário, do inglês
2O zero está excluido pela normalização
Em
El
ab
or
aç
ão
1.9. Aritmética de Ponto Flutuante 9
1.00 1.00
1.01 1.25
1.10 1.50
1.11 1.75
binário decimal
Tabela 1.2: Conversão
Para os expoentes, temos
e = −1, 0, 1, 2, 3
Calculando os números positivos deste sistema teremos
e
−1 0 1 2 3
1.00 · 2e ⇒ 0.5 1 2 4 8
1.25 · 2e ⇒ 0.625 1.25 2.5 5 10
1.5 · 2e ⇒ 0.75 1.5 3 6 12
1.75 · 2e ⇒ 0.875 1.75 3.5 7 14
Tabela 1.3: Números positivos do sistema
Dispondo numa reta,
x
0
0.5
1
1.5
2
2.5
3
3.5
4 5 6 7 8 10 12 14
Observe que, de acordo com o que foi dito antes, a qantidade total de números deste sistema é:
2(β− 1)βt−1(emax − emin + 1) = 2 · 1 · 22 · 5 + 1 = 41
Além disso,
β1−t = 2−2 = 0.25 (menor distância possível)
1.11 · 23 = 1.75 · 8 = 14 (maior número possível)
1.00 · 2−1 = 0.5 (menor número possível)
�
Overflow e Underflow
Como pôde ser visto no exemplo anterior, 1.9.1, a quantidade de números possíveis em F(β, t, emin, emax)
é finito e, diante disto, fica clara a impossibilidade de representar a continuidade dos números reais.
Qualquer número maior que o valor máximo do sistema F, não é representável em F. Esta limitação
leva aos fenômenos conhecidos como overflow. Para muitos propósitos, tal impossibilidade também se
apresenta para números diferentes de zero e menores que o menor valor de F. Este tipo de ocorrência é
também conhecida como underflow.
Estas limitações ocorrem por dois motivos:
i) Os expoentes de F estão num intervalo limitado;
ii) Número finito de dígitos para a mantissa, o que leva a valores para a mesma, no intervalo [β−1,β− β−t]
Em
El
ab
or
aç
ão
10 1. Sistemas Numéricos
Exemplo 1.9.2
Considere o sistema F(10, 4,−5, 5) e sejam x = 875 e y = 3172. Vejamos o que ocorre quando
calcula-se x · y. Observe inicialmente que os números x e y, representados em F, possuem o seguinte
formato:
x = 8.750× 102
y = 3.172× 103
e
xy = 2.776 · 106 (Overflow)
uma vez que o expoente máximo do sistema é 5. �
Exemplo 1.9.3
Considere novamente o sistema F(10, 4,−5, 5) e sejam x = 0.0064 e y = 7312. Desejamos agora,
observar o que ocorre quando calculamos x/y. Escrevendo x e y no sistema F, teremos
x = 6.400× 10−3
y = 7.312× 103
Assim:
x
y
= 8.753 · 10−7 (Underflow)
Pois, o menor expoente possível em F é −5. �
Na ocorrência de overflow e underflow, cada implementação reage de alguma forma. As principais são:
• Overflow
– interrompe-se o cálculo;
– retorna o número que representa o infinito da máquina.
• Underflow
– interrompe-se o cálculo;
– arredonda-se o resultado para zero;
– arredonda-se para um número subnormal3.
As duas maneiras como as máquinas tratam o overflow possuem aspectos indesejáveis. O primeiro
caso não apresenta resposta enquanto que, no segundo caso, apesar de dar uma resposta, não é muito útil,
exceto para interpretações físicas ou aplicações específicas, uma vez que tal resposta não possui valor
numérico.
Na ocorrência de underflow, o primeiro caso, pelo mesmo motivo apresentado para o overflow, é
também indesejado. Os segundo e terceiro casos resultam em respostas úteis embora exista o perigo de,
numa operação subsequente, ocorrer um overflow.
Em geral, tenta-se evitar o overflow e underflow em implementações computacionais e, uma maneira
prática de se fazer isto é utilizando uma técnica conhecida como escalonamento.
Exemplo 1.9.4 (Escalonamento)
3Veremos posteriormente o que são números subnormais.
Em
El
ab
or
aç
ão
1.9. Aritmética de Ponto Flutuante 11
Considere o sistema F(10, 10,−99, 99) e os números a = 1060 e b = 1060. Calcular c =
√
a2 + b2.
Observe que, o maior número deste sistema é
10 − 10−99 . 1099︸ ︷︷ ︸
(β−β−emin) .βemax
= 10100 − 1
e
a2 = b2 = 10120 (overflow)
Porém,
c =
√
a2 + b2 =
√
2 . 10120 =
√
2 . 1060
é representável em F. Para evitar este problema, o escalonamento consiste em observar que
c =
√
a2 + b2
=
√
s2 (a
2 + b2)
s2
= s
√(a
s
)2
+
(
b
s
)2
(1.9.1)
Com isto, se tomarmos
s = |a|+ |b| = 2 · 1060
teremos
a
s
=
b
s
=
1060
2 · 1060 =
1
2
⇒
(
a
s
2
)
=
(
b
s
2
)
= 0.125
Portanto, s,
(
a
s
)2 e
(
b
s
)2 são todos representáveis em F e em consequência disto, o número escrito na
forma 1.9.1, também.
Erros
Chamamos de número de máquina a todo número que pode ser exatamente representado na aritmética
de ponto flutuante em consideração.
Em geral, operações aritméticas envolvendo números da máquina, podem ainda assim, ter como re-
sultados, números que não são exatamente representáveis. Surge portanto a necessidade de se aproximar
tais resultados por números da máquina, gerando com isto, erros de aproximação ou erros de arre-
dondamento.
Existem 4 maneiras para se realizar esta aproximação. Sendo x um número a ser aproximado, pode-se
utilizar os seguintes métodos:
• Truncamento: O valor escolhido, denotado aqui por ∇(x) é o maior número da máquina menor
que x.
• Arredondamento: Consiste em escolher o menor número da máquina maior que x. Denota-se ∆(x).
• Em direção ao 0: Denotado por Z(x), é igual a ∇(x) se x > 0 e igual a ∆(x) se x < 0
• Vizinho mais próximo: Consiste em escolher o número de máquina, denotado por N(x), que está
mais próximo de x. Se x estiver exatamente no meio de dois números da máquina, qualquer um dos
dois pode ser o escolhido.
Em
El
ab
or
aç
ão
12 1. Sistemas Numéricos
x0
∇(x)
Z(x),N(x)
∆(x)
(a)
x
0
∇(x)
Z(x) N(x)
∆(x)
(b)
0
∇(x)
N(x) Z(x)
∆(x)
(c)
Erros Absoluto e Relativo
A exatidão de qualquer método de aproximação é medida em termos do erro cometido. Existem duas
maneiras de se especificar o erro: erro absoluto e erro relativo.
Sendo x um número real e x uma aproximação, define-se erro absoluto, como:
Ea = |x− x|
e o erro relativo como:
Er =
|x− x|
|x|
=
Ea
|x|
O erro absoluto mede a magnitude do erro, enquanto o erro relativo determina o quão grande é o erro
cometido em relação ao valor correto.
Exemplo 1.9.5
Sejam x = 100 e y = 0.0006 e suponha que estes números foram aproximados por x = 100.1 e
y = 0.0004. Observe que
Eax = |x− x| = 0.1
Eay = |y− y| = 0.0002
ou seja, o erro absoluto cometido na aproximação de y é muito menor do que o erro absoluto cometido
na aproximação feita para x. Por outro lado, temos que:
Erx =
Eax
|x|
=
0.1
100
= 0.001
Ery =
Eay
|y|
=
0.0002
0.0006
=
1
3
Desta forma, podemos concluir que a aproximação feita para x é melhor que a aproximação feita para
y. �
Precisão da Máquina
O maior erro relativo possível, cometido quando um número real é apresentado numa aritmética de ponto
flutuante é chamado de precisão da máquina. Como é comum usar o símbolo ε para denotar a precisão
da máquina, é também usual se referir a este número como epsílon da máquina.
Outros autores definem o ε da máquina como o menor número, tal que
1 + ε > 1
Proposição 1.9.1. Seja F(β, t, emin, emax) uma aritmética de ponto flutuante. Então
ε 6
1
2
β1−t
Em
El
ab
or
aç
ão
1.9. Aritmética de Ponto Flutuante 13
Demonstração. Suponha que um dado x ∈ R seja aproximado por
x = a1 · a2a3 · · ·at · βe
onde e ∈ [emin, emax]. O erro absoluto cometido será
|x− x| < 0 · 00 · · · 0︸ ︷︷ ︸
t
β‘ · βe (1.9.2)
onde β‘ é o dígito na base β tal que
β‘ =
β
2
ou seja
|x− x| <
β
2
· β−t · βe = 1
2
β1−t · βe
Observe que, para todo x escrito na forma (1.9.2), tem-se que
βe < |x| < βe+1
e disto segue-se que:
1
βe+1 <
1
|x|
<
1
βe
Logo,
ε =
|x− x|
|x|
<
1
2β
1−t · βe
βe
=
1
2
β1−t
�
Embora o ε da máquina esteja, conforme vimos, precisamente especificado numa aritmética de ponto
flutuante, não é crucial que o usemos de forma exata. É comum que seu valor seja aproximado por um
número calculado através do seguinte código:
Algoritmo 5: Epsilon da Máquina
Entrada: eps = 1; eps1 = 0;
início
repita
eps = 0.5 · eps;
eps1 = eps+ 1;
até (eps1 < 1.0);
retorne eps;
fim
Números Subnormais e Exceções
Em sistemas de ponto flutuante, númerossubnormais são aqueles que estão no espaço vazio em torno do
zero. Portanto, qualquer número diferente de zero, com módulo menor que o menor número normalizado
do sistema é um número subnormal.
Em aritméticas de ponto flutuante normalizadas, o primeiro dígito da mantissa é sempre diferente
de zero. Com o intuito de amenizar o underflow, os números subnormais são implementados com
expoente 0 e mantissa não normalizada.
Também são implementadas representações especiais para algumas exceções, tais como:
Em
El
ab
or
aç
ão
14 1. Sistemas Numéricos
• NaN (Not a Number)4 é o resultado de operações aritméticas inválidas, por exemplo,
√
−5, ∞∞ , · · · ;
• ± ∞ resultante de um overflow ou divisão por zero;
• ± 0 resultante de um underflow, ou divisão por ±∞.
1.10 Instabilidades Numéricas
De um modo informal, podemos dizer que um processo numérico é instável se pequenos erros cometidos
em algum estágio do processo, são ampliados em estágios subsequentes e degradam seriamente o resultado
final.
Como ilustração deste conceito, considere a sequência de números reais definida indutivamente por:
x0 = 1
x1 =
1
3
xn+1 =
13
3
xn −
4
3
xn−1, n > 1
(1.10.1)
Esta relação de recorrência gera a sequência
xn =
(
1
3
)n
se usarmos a relação de recorrência (1.10.1), para calcular os termos desta sequência, numa aritmética de
ponto flutuante F(24, 2,−64, 64)5, teremos os seguintes resultados, já convertidos para o sistema decimal:
x0 = 1.0000000
x1 = 0.3333333 (7 dı́gitos corretos)
x2 = 0.1111112 (6 dı́gitos corretos)
x3 = 0.0370373 (6 dı́gitos corretos)
x4 = 0.0123466 (5 dı́gitos corretos)
x5 = 0.0041187 (5 dı́gitos corretos)
x6 = 0.0013857 (4 dı́gitos corretos)
x7 = 0.0005131 (0 dı́gitos corretos)
...
x15 = 3.6574980 (erro relativo de 108)
Este algoritmo é portanto instável. Qualquer erro presente em xn é multiplicado por 13
3 ao calcular
xn+1. Existe portanto a possibilidade de que o erro cometido no cálculo de x1 seja propagado para x15
com um fator de ampliação de
(13
3
)14.
Um processo é classificado como numericamente instável ou não, em função do erro relativo final,
cometido em sua execução. Assim, mesmo o erro absoluto sendo grande, seu resultado pode ser aceitável,
desde que se trate de um valor alto.
No exemplo anterior, considere x0 = 1 e x1 = 4 em (1.10.1). A relação de recorrência permanece a
mesma e os erros serão propagados e ampliados do mesmo modo. Entretanto, a sequência gerada, agora é
xn = 4n
4Não é um número, em inglês.
5Sistema Binário de 32 bits.
Em
El
ab
or
aç
ão
1.10. Instabilidades Numéricas 15
e os resultados obtidos, quando calculamos os termos dessa sequência, usando (1.10.1), são, por exemplo:
x1 = 4.000006
...
x10 = 1.048576× 106
...
x20 = 1.099512× 1012
E neste caso, em razão da grandeza dos números obtidos, os erros cometidos são relativamente pequenos,
e portanto negligenciáveis.
Condicionamento
O termo condicionamento é em geral, utilizado para indicar a sensibilidade da solução de um problema
em relação à pequenas alterações nos dados de entrada. Um problema é considerado mal condicionado
se pequenas mudanças nos dados resultam em grandes mudanças nas respostas do mesmo.
Em certos problemas, é possível definir um número de condicionamento. Quando este número é
grande, indica um problema mal condicionado.
Como ilustração, considere o problema de se calcular uma função f, de uma variável, em x. Qual é
o efeito de uma pequena perturbação em x, no cálculo de f(x)? se quisermos calcular o erro absoluto,
usando o Teorema do Valor Médio, temos que
f(x+ ∆x) − f(x) = f′(ξ) · ∆x ≈ f′(x) · ∆x
Portanto, se f′(x) não é muito grande, o efeito desta perturbação em f(x) será pequeno. Ao realizar uma
perturbação ∆x em x temos um erro relativo de
∆x
x
. Do mesmo modo, quando f(x) é perturbada e torna-se
f(x+ ∆x), o erro relativo será
f(x+ ∆x) − f(x)
f(x)
≈ ∆x · f
′(x)
f(x)
=
[
xf′(x)
f(x)
]
· ∆x
x
E assim, o fator
xf′(x)
f(x)
funciona como um número de condicionamento para este problema.
Exemplo 1.10.1
Qual é o número de condicionamento no cálculo da função arcsen(x)? Considere
f(x) = arcsen(x)
De acordo com o que vimos anteriormente, o número de condicionamento para este tipo de problema
é o fator:
xf ′(x)
f(x)
=
x
arcsen(x)
√
1 − x2
Assim, para x próximo de 1, por exemplo, temos arcsen(x) ≈ π
2
e o número de condicionamento tende
a +∞ à medida que x tende a 1. Portanto, pequenas perturbações em x produzam grandes alterações
em f(x), quando está próximo de 1.
Em
El
ab
or
aç
ão
16 1. Sistemas Numéricos
1.11 ExercÃcios
1.11.1 ???
Em
El
ab
or
aç
ão
2
Zeros Reais de Funções Reais
Em diversas áreas da ciências ocorre situações que é necessário obter um determinado valor de uma função,
recaindo numa equação do tipo
f(x) = 0
os valores de x que anulam a função são chamados de raízes ou zeros da função.
As equações de 1◦, 2◦, 3◦ e 4◦ grau possuem soluções exatas, para o caso em que o grau > 5◦ não há
solução geral. Nestes casos o problema só pode ser resolvido numericamente.
Exemplo 2.0.1
FIGURA
ξ = Ri+ g(i)
ξ− Ri− g(i)︸ ︷︷ ︸ = 0
f(i) = 0
Meu problema de encontrar a corneta do circuito se reduz a encontrar as raízes da função f(i).
Podemos analisar exemplos envolvendo a bolsa de valores, produção de objetos, e assim por diante.
Nosso trabalho será obter as raízes reais desta função. A partir dos métodos numéricos, podemos
encontrar soluções com a precisão desejada.
Os métodos consistem de duas fases:
Fase I - Localização ou Isolamento das Raízes: encontrar um intervalo que contém a raíz.
Fase II - Refinamento: melhorar sucessivamente a solução até a precisão desejada.
2.1 Fase I: Isolamento das Raízes
Nesta fase devemos determinar a região que iremos procurar a raiz.
Teorema 2.1.1 Seja f(x) uma função contínua em [a,b]. Se f(a) · f(b) < 0 então existe pelo menos um
ponto x = ξ entre a e b que f(ξ) = 0.
Ideia: FIGURA
f(a) · f(b) < 0
17
Em
El
ab
or
aç
ão
18 2. Zeros Reais de Funções Reais

f(a) > 0 e f(b) < 0
ou
f(a) < 0 e f(b) > 0
Como é contínua, em algum momento ela passou pelo zero.
Observação 2.1.1 Se f(a) · f(b) < 0 e f ′(x) preservar o sinal, para todo x ∈ (a,b) então a raiz é
única.
A derivada f ′(x) preservar o sinal implica na função ser unicamente
{
crescente, f ′(x) > 0
decrescente, f ′(x) < 0 no
intervalo.
Como a função só cresce (ou decresce) no intervalo:
FIGURA
Uma forma de se isolar as raízes de f(x) é analisar um conjunto de valores para a função f(x) e sua
derivada f ′(x).
Exemplo 2.1.1
a) f(x) = x3 − 9x+ 3 ; f ′(x) = 3x2 − 9
FIGURA
x −100 −10 −5 −3 −1 0 1 2 3 4 5
f(x)
− − − + + + − − + + +
−907 −77 3 11 3 −5 −7 3 31 83
f ′(x)
+ + + + − − − + + + +
299 16 18 −6 −9 −6 3 18 39 66︸ ︷︷ ︸ ︸ ︷︷ ︸ ︸ ︷︷ ︸
Nessa região a função mudou de sinal, ou seja, há uma raiz,
e a derivada não altera o sinal, então a função é crescente
durante todo o intervalo, ou seja, só há uma raiz.
Como o polinômio é de grau 3, possui no máximo 3 raízes.
b) f(x) =
√
x− 5e−x ; f ′(x) =
1√
x
+ 5e−x
x 0 1 2 3 4
f(x)
− − + + +
−5 −0, 84 0, 74 3, 48 1, 91
f ′(x) Inf
+ + +
2, 34 1, 03 0, 54 0, 34
Só há uma raiz.
FIGURA
Em
El
ab
or
aç
ão
2.2. Fase II: Refinamento 19
Observação 2.1.2 É equivalente analisar as funções f(x) = 0 ou f(x) = g(x) − h(x) = 0⇒ g(x) =
h(x), ou seja
a) f(x) = x3 − 9x+ 3 = 0 ou x3 = 9x− 3 FIGURA
b)
√
x− 5e−x ou
√
x = 5e−x
2.2 Fase II: Refinamento
Após definir a região onde há uma raiz, vamos encontrar o valor a partir do Método Iterativo.
O Método Iterativo consiste num procedimento ou uma sequência de instruções, reduzidas passo a
passo. Cada iteração (ciclo de procedimentos) baseia-se no resultado anterior, para encontrar o próximo.
Os métodos iterativos em geral seguem o seguinte diagrama de fluxo:
DIAGRAMA
O critério de parada pode ser:
i) |x− ξ| < ε
ii) |f(x)| < ε
;

x é a raiz aproximada
ξ é a raiz exata
ε é a precisão desejada
Observação 2.2.1 Apesar da opção i) da a impressão de termos que conhecer a raíz, o que estamos
na verdaderequerindo é que a raíz aproximada esteja em um intervalo de confiança.
FIGURA
Os critérios i) e ii) nem sempre são satisfeitos simultaneamente, podem haver casos como
FIGURA
|f(x)| < ε, satisfeito
|x− ξ| > ε, não satisfeito .
FIGURA
|x− ξ| < ε, satisfeito
|f(x)| > ε, não satisfeito .
2.3 Método da Bissecção
Seja f(x) contínua no intervalo [a,b] e f(a) · f(b) < 0. Vamos supor (por simplicidade) que só existe uma
raiz neste intervalo e que f(a) < 0 e f(b) > 0.
FIGURA
O primeiro método de solução consiste em ir dividindo sucessivamente o intervalo ao meio, até
• Encontrar a raiz (sorte);
• Obter a precisão requerida x, ou seja, x ∈ ]ξ− ε; ξ+ ε[.
x0 =
a+ b
2

f(a0) < 0
f(b0) > 0
f(x0) > 0
;
então a raiz está entre a e x0; ξ ∈ ]a0; x0[
novo intervalo:
{
a1 = a0
b1 = x0
x1 =
a1 + b1
2

f(a1) < 0
f(b1) > 0
f(x1) < 0
;
então a raiz está entre x1 e b1; ξ ∈ ]x1;b1[
novo intervalo:
{
x1 = a2
b1 = b2
Em
El
ab
or
aç
ão
20 2. Zeros Reais de Funções Reais
x2 =
a2 + b2
2

f(a2) < 0
f(x2) < 0
f(b2) > 0
;
ξ ∈ ]x2;b2[
novo intervalo:
{
a3 = x2
b3 = b2
Observação 2.3.1 Se eu quero uma precisão de ε, meu intervalo tem que ter no máximo tamanho
ε, ou seja
bn − an = ε.
Exemplo 2.3.1
Seja f(x) = x logx−1
FIGURA
f(2) = −0, 40 < 0 ; f(3) = 0, 43 > 0. Então ξ ∈ [2; 3], a0 = 2;b0 = 3.
x0 =
a0 + b0
2
=
2 + 3
2
= 2, 5

f(2) = −0, 4
f(2, 5) = −0, 0051
f(3) = 0, 4
a1 = 2, 5
b1 = 3
x ∈ ]2, 5; 3[
x1 =
a1 + b1
2
= 2, 75

f(2, 5) = −0, 0051
f(2, 75) = 0, 21
f(3) = 0, 4
a2 = 2, 5
b2 = 2, 75
x ∈ ]2, 5; 2, 75[
x2 =
a2 + b2
2
= 2, 625

f(2, 5) = −0, 0051
f(2, 625) = 0, 10
f(2, 75) = 0, 4
a3 = 2, 5
b3 = 2, 625
x ∈ ]2, 5; 2, 625[
x3 =
a3 + b3
2
= 2, 5625

f(2, 5) = −0, 0051
f(2, 5625) = 0, 047
f(2, 625) = 0, 10
a4 = 2, 5
b4 = 2, 5625
x ∈ ]a4;b4[
...
Até encontrar a raiz na precisão desejada, isto é,
x ∈ ]an,bn[; tal que bn − an = ε.
2.4 Estimativa do Número de Iterações
Uma pergunta que surge é, quantas iterações serão necessárias para atingir uma certa precisão?
Como vamos dividindo o intervalo sucessivamente ao meio
bk − ak =
bk−1 − ak−1
2
=
b0 − a0
2k
Para termos um precisão de ε, temos que
bk − ak < ε⇒
b0 − a0
2k
< ε
⇒ 2k
b0 − a0
> ε
⇒ 2k >
b0 − a0
ε
Em
El
ab
or
aç
ão
2.5. Método da Posição Falsa 21
aplicando o logarítmo
log 2k > log
(
b0 − a0
ε
)
⇒ k log2 > log(b0−a0)− logε
k >
log(b0−a0)− logε
log(2)
Portanto, no Exemplo 2.3.1, com f(x) = x log x− 1 e x ∈ [2, 3], para termos uma precisão de
• ε = 10−2 :
k >
log(3−2)− log10−2
log2 =
0 + 2 log10
log2
= 6, 64⇒ k > 7
• ε = 10−4 :
k >
3
log2 = 13, 28⇒ k > 14
Observação 2.4.1 É ummétodo simples, que não envolve cálculos laboriosos, mas em contrapartida,
a convergência é muito lenta, já que vamos sempre dividindo o intervalo ao meio.
Exemplo 2.4.1
Se b0 − a0 = 10 e ε = 10−8
k >
log(b0−a0)− logε
log(2) =
log(10)− log10−8
log(2)
=
1 + 8
log(2)
=
9
log(2) = 29, 89⇒ k > 30
2.5 Método da Posição Falsa
Seja novamente f(x) contínua no intervalo [a,b] e f(a) · f(b) < 0, com uma única raiz neste intervalo.
No método da bissecção consideramos xn =
an + bn
2
a cada iteração. Em vez de considerar a média
aritmétrica entre a e b podemos considerar média ponderada, levando em conta a posição relativa dos
pontos.
x =
a|f(b)|+ b|f(a)|
|f(b)|+ |f(a)|
Se f(b) > 0 e f(a) < 0, podemos reescrever a relação como
Em
El
ab
or
aç
ão
22 2. Zeros Reais de Funções Reais
x =
af(b) − bf(a)
f(b) − f(a)
, para f(b) > 0 e f(a) < 0.
Graficamente, esse ponto x é a reta que passa por (a, f(a)) e (b, f(b)) em y = 0.
FIGURA
x =
af(b) − bf(a)
f(b) − f(a)
De fato, seja
y = y0 + α · x
então,
y = 0⇒ x =
−y0
α
Como a reta passa por (a, f(a)) e (b, f(b)),
{
f(a) = y0 + α · a
f(b) = y0 + α · b
• f(b) − f(a) = α(b− a)⇒ α =
f(b) − f(a)
b− a
•

b · f(a) = y0 · b+ αab
⇒
a · f(b) = y0 · a+ αab
a · f(b) − b · f(a) = y0(a− b)
y0 =
a · f(b) − b · f(a)
a− b
Portanto,
x =
a · f(b) − b · f(a)
a− b
f(b) − f(a)
b− a
=
a · f(b) − b · f(a)
f(b) − f(a)
Exemplo 2.5.1
Seja f(x) = x log(x) − 1, com uma raiz entre [2, 3], então a0 = 2, b0 = 3.
f(2) = −0, 40 ; f(3) = 0, 43
x0 =
a · |f(b)|+ b · f|a|
|f(b)|+ |f(a)|
=
2 · 0, 43 + 3 · 0, 40
0, 40 + 0, 43
= 2, 48

f(2) = −0, 40 < 0
f(2, 48) = −0, 022 < 0
f(3) = 0, 43 > 0
a1 = 2, 48
b1 = 3
f(2, 48) = −0, 02
f(b1) = −0, 43
f(a1) = −0, 02
x1 =
a · |f(b)|+ b · f|a|
|f(b)|+ |f(a)|
=
2, 48 · 0, 43 + 3 · 0, 02
0, 43 + 0, 02
= 2, 50

f(2, 48) = −0, 022 < 0
f(2, 50) = −0, 0051 < 0
f(3) = 0, 43 > 0
a2 = 2, 50
b2 = 3
f(2, 50) = −0, 0051
f(b2) = 0, 43
f(a2) = −0, 0051
Em
El
ab
or
aç
ão
2.6. Método do Ponto Fixo 23
x2 =
a · |f(b)|+ b · f|a|
|f(b)|+ |f(a)|
=
2, 50 · 0, 43 + 3 · 0, 0051
0, 43 + 0, 0051
= 2, 51

f(2, 50) = −0, 0051 < 0
f(2, 51) = −0, 0051 > 0
f(3) = 0, 43
{
a3 = 2, 50
a4 = 3
f(2, 51) = 0, 0032
Observação 2.5.1 A aproximação que estamos tomando influência o resultado.
2.6 Método do Ponto Fixo
Os três próximos métodos

Método do Ponto Fixo
Método de Newton
Método das Secantes
são baseados na definição de ponto fixo, ou seja,
f(x) = x, e no seguinte Teorema:
Teorema 2.6.1 (Ponto Fixo) Seja f(x) : [a,b] → [a,b] uma função contínua com |f ′(x)| < 1, ∀ x ∈
[a,b], então
(i) Existe um único ponto fixo x∗ ∈ [a,b], ou seja, f(x∗) = x∗;
(ii) Qualquer sequência da forma xn+1 = f(xn) com x0 ∈ (a,b) converge para o ponto fixo x∗.
Exemplo 2.6.1
Seja f(x) = x2 com x ∈
(
−1
2 , 1
2
)
. f(0) = 0 é um ponto fixo e f ′(x) = 2x < 1, ∀ x ∈
(
−1
2 , 1
2
)
então
x =
1
4
; x1 = f(x0) = f
(
1
4
)
=
1
16
= 0, 0625
x2 = f(x1) = f
(
1
16
)
=
1
256
= 0, 0039
x3 = f(x2) = f
(
1
256
)
=
1
65536
= 0, 000015
Os próximos 3 métodos consistiram em reescrever a função f(x) de forma a explicitar um ponto fixo e
em seguida calcula-se a partir de diversas iterações.
f(x) = 0⇔ ϕ(x) = x
a função ϕ(x) é chamada de função de iteração.
Exemplo 2.6.2
Seja f(x) = x2 + x− 6 = 0, podemos ter várias funções de iteração.
i) Com x2 + x− 6 = 0⇒ x = 6 − x2, temos a função de iteração ϕ1(x) = 6 − x2
ii) Com x2 + x− 6 = 0⇒ x2 = 6− x⇒ x = ±
√
6 − x, temos a função de iteração ϕ2(x) = ±
√
6 − x
iii) Com x2 + x− 6 = 0⇒ x2 = 6− x⇒ x = 6−x
x = 6
x − 1, temos a função de iteração ϕ3(x) =
6
x − 1
iv) Com x2 + x− 6 = 0⇒ x(x+ 1) = 6⇒ x = 6
x+1 , temos a função de iteração ϕ4 = 6
x+1
Em
El
ab
or
aç
ão
24 2. Zeros Reais de Funções Reais
De modo geral, as funções de iteração ϕ(x) são da forma ϕ(x) = x + A(x)f(x); ϕ(ξ) = ξ em que
f(ξ) = 0 e A(ξ) 6= 0.
Observação 2.6.1 A segunda condição é importante, pois implica em
f(ξ) = 0⇔ ϕ(ξ) = ξ
⇒ ϕ(ξ) = ξ+A(ξ)f(ξ) = ξ
⇒ A(ξ)f(ξ) = 0
como A(ξ) 6= 0 temos que f(ξ) = 0.
Portanto, satisfazendo estas condições, temos a equivalência.
Observação 2.6.2 Dada uma função f(x), podemos ter infinitas funções de iteração ϕ(x) para a
equação f(x) = 0, mas só utilizaremos as que satisfazem as condições do Teorema do Ponto Fixo.
Graficamente:
• f(x) satisfaz a condição de que |f ′(x)| < 1. A inclinação da curva é menor que a da bissetriz. FIGURA
• f(x) satisfaz a condição de que |f ′(x)| > 1.
FIGURA
A condição de que |f ′(x)| < 1 implica na convergência do valor.
Observação: no exemplo anterior f(x) = x2 + x − 6, sabendo que há uma raiz entre [1, 3], temos que{
f(1) = −4
f(3) = 6
i) ϕ3(x) = 6 − x2 ⇒ ϕ
′
3(x) = −2x;
ii) ϕ2(x) =
√
6 − x⇒ ϕ
′
2(x) =
−1
2
√
6 − x
; |ϕ ′2(x)| < 1, ∀ x ∈ [1, 3];
iii) ϕ3(x) =
6
x − 1⇒ ϕ
′
3(x) =
−6
x2 (não satisfaz a condição para o intervalo [1, 3]) Observação: pode ser
útil para a outra raiz ξ = −3 ∈ [−3, 5; 2, 5];
iv) ϕ4(x) =
6
x+ 1
⇒ ϕ ′(x) =
−6
(x+ 1)2
– Para a região [1, 3], não satisfaz ϕ ′(x) =
−6
4
– Para a região [1, 5, 3] satisfaz ϕ ′(1, 5) =
−6
2, 52 = −6
6,25 = −0, 96
Então ii) ϕ2(x) =
√
6 − x; x ∈ [1, 3]. Escolhendo um número aleatório entre 1 e 3.
x0 = 2, 2
x1 = ϕ(x0) = ϕ(2, 2) = 1, 949
x2 = ϕ(x1) = ϕ(1, 949) = 2, 0127
x3 = ϕ(x2) = 1, 996822
x4 = ϕ(x3) = 2, 000794
Em
El
ab
or
aç
ão
2.7. Método de Newton-Raphson 25
iv) ϕ4(x) =
6
x+ 1
;x ∈ [1, 5; 3] ; x0 = 2, 2
x1 = ϕ(x0) = 1, 875
x2 = ϕ(x1) = 2, 08696
x3 = ϕ(x2) = 1, 94366
x4 = ϕ(x3) = 2, 03828
Observação 2.6.3 Ordem de convergência pode ser obtida analisando que{
ϕ(xk) −ϕ(ξ) = xk+1 − ξ ; já que xk+1 = ϕ(xk) e ϕ(ξ) = ξ
ϕ(xk) −ϕ(ξ) = ϕ
′(ck) · (xk − ξ)
então xk+1−ξ = ϕ(k)−ϕ(ξ) =
ϕ ′(ck)(xk − ξ)⇒
xk+1 − ξ
xk − ξ
= ϕ ′(ck) < 1. Denotando
{
xk+1 − ξ = ck+1
xk − ξ = ck
o erro na iteração k + 1
e na iteração k respectivamente, temos que
ck+1
ck
= ϕ ′(ck) < 1⇒ ck+1 ≈ c · ck a
tomando o limite de k→ +∞
lim
k→+∞ ck+1
ck
= lim
k→+∞ϕ ′(ck) = ϕ ′(ξ) (em torno da raiz)
Portanto, após um número grande de iterações (k → +∞) a iteração anterior é proporcional a
iteração posterior, com a constante de proporcionalidade sendo ϕ ′(ξ).
aA aproximação é dita ser linear já que ck+1 ≈ ccpk com p = 1, onde p é o denominador de ordem de convergência.
2.7 Método de Newton-Raphson
No método a convergência depende de ϕ ′(ξ). O método de Newton consiste em tomar ϕ ′(ξ) = 0, a fim
de acelerar a convergência.
Seja ϕ(x) = x+A(x) + f(x) a função de iteração
ϕ ′(x) = 1 +A ′(x)f(x) +A(x)f ′(x),
tomando x = ξ,
ϕ ′(ξ) = 1 +A ′(ξ)f(ξ) +A(ξ)f ′(ξ)⇒ ϕ ′(ξ) = 1 +A(ξ)f ′(ξ)
Para ϕ ′(ξ) = 0 temos que
0 = 1 +A(ξ)f ′(ξ)⇒ A =
−1
f ′(ξ)
O oposto inverso da derivada de f(x).
Se considerarmos A(x) =
−1
f ′(x)
temos que
ϕ(x) = x−
f(x)
f ′(x)
, com f ′(x) 6= 0, x ∈ [a,b].
Portanto, o método de Newton-Raphson consiste em fixar uma função de iteração em que
xk+1 = xk −
f(xk)
f ′(xk)
Em
El
ab
or
aç
ão
26 2. Zeros Reais de Funções Reais
Exemplo 2.7.1
Seja f(X) = x2 + x− 6; f ′(x) = 2x+ 1. Então,
ϕ(x) = x−
x2 + x− 6
2x+ 1
, seja x0 = 2, 2
x1 = ϕ(x0) = 2, 2 −
4, 84 + 2, 2
4, 4 + 1
= 2, 2 −
1, 04
5, 4
= 2, 2 − 0, 19259 = 2, 0074
x2 = ϕ(x1) = 2, 0074 −
2, 00742 + 2, 0074 − 6
2(2, 0074) + 1
= 2, 0074 −
0, 03705
5, 0148
= 2, 00001092
Com duas iterações, conseguimos uma grande aproximação da raiz real, que é 2.
2.8 Método da Secante
Uma grande desvantagem no Método de Newton-Raphson é a necessidade de se obter f ′(x). Uma forma
de se contornar tal problema é substituir a derivada f ′(x) pelo quociente das diferenças.
f ′(xk) =
f(xk) − f(xk−1
xk − xk−1
A função de iteração é:
ϕ(xk) = xk −
f(xk)
f(xk)−f(xk−1)
xk−xk−1
=
xk(f(xk) − f(xk−1)) − f(xk)(xk − xk−1)
f(xk) − f(xk−1)
=
xk−1f(xk) − f(xk−1)xk
f(xk − f(xk−1)
Observação 2.8.1 São necessários dois pontos iniciais, para se iniciar o método.
Exemplo 2.8.1
Seja f(x) = x2 + x− 6, raízes
{
−3
2 , tomando x0 = 1, 5 e x1 = 1, 7,

f(x0) = −2, 25
f(x1) = −1, 41
f(x2) = 0, 17983
, então
xk+1 =
xk−1f(xk) − f(xk−1)xk
f(xk) − f(xk−1)
x2 =
x0f(x1) − f(x0)x1
f(x1) − f(x0)
=
1, 5x(−1, 41) − 1, 7(−2, 25)
−1, 41 + 2, 25
= 2, 03571
x3 =
x1f(x2) − f(x1)x2
f(x2) − f(x1)
=
1, 7x(0, 17983) − (−1, 41)(2, 03575)
0, 17983 − (−1, 41)
= 1, 99774 ; f(x3) = −0, 01131
Com poucas iterações, temos uma boa convergência!
x4 =
x2f(x3) − f(X2)x3
f(x3) − f(x2)
=
2, 03571(−0, 01131 − 0, 17983(1, 99774)
−0, 01131 − 0, 17983
= 1, 99999
Em
El
ab
or
aç
ão
2.9. Métodos Estudados 27
2.9 Métodos Estudados
• Método da Bissecção: xk =
ak + bk
a
(ζ ∈ [a,b])
• Método da Posição Falsa: xk =
akf(bb) − bkf(ak)
f(bk) − f(ak)
• Método do Ponto Fixo: xk+1 = ϕ(xk) ; com |ϕ ′(x)| < 1 (função de iteração)
• Método de Newton: xk+1 = xk −
f(xk)
f ′(xk)
• Método da Secante: xk+1 = xk −
f(xk)
f(xk)−f(xk−1)
xk−xk−1
=
xk−1f(xk) − xkf(xk−1)
f(xk) − f(xk−1)
Observação 2.9.1 Na implementação dos métodos deve ser levado em conta o esforço computaci-
onal. A escolha do método está diretamente ligado a equação que queremos resolver.
Em
El
ab
or
aç
ão
3
Sistemas de Equações Lineares
3.1 Introdução
A solução de muitos problemas oriundos das Engenharias necessitam da resolução de um ou mais sistemas
lineares. Por exemplo, os problemas de cálculo de estruturas na engenharia civil; o cálculo da intensidade
de corrente elétrica em uma malha, em engenharia elétrica; e o problema de programação linear, em
engenharia de produção, são exemplos de situações nas quais a resolução de pelo menos um sistema linear
pode ser necessária.
Neste capítulo, trataremos da resolução numérica de sistema lineares que possuem o número de equa-
ções igual ao número de incógnitas. Serão abordados três métodos diretos, (Eliminação Gaussiana, Fato-
ração LU e Fatoração de Cholesky) e dois métodos iterativos (o método de Gauss-Jacobi e o método de
GAuss-Sidel).
3.2 Sistemas Lineares
Definição 3.2.1 Um sistema de equações lineares, ou simplismente sistema linear, com n equações
e n incógnitas é um conjunto de equações do tipo:
a11x1 + a12x2 + · · ·+ a1nxn = b1
a21x1 + a22x2 + · · ·+ a2nxn = b2
...
an1x1 + an2x2 + · · ·+ annxn = bn
(3.2.1)
sendo aij números reais e 1 6 i 6 n e 1 6 j 6 n. No caso em que b1 = b2 = · · · = bn = 0, o
sistema linear é chamado de Sistema Linear Homogêneo. Uma solução do sistema 3.2.1 é uma n-upla
(x1, x2, . . . , xn) de números reais que satisfazem as n equações simultaneamente.
O sistema linear 3.2.1 também pode ser representado por meio da seguinte equação matricial
a11 a12 · · · a1n
a21 a22 · · · a2n
...
...
. . .
...
an1 an2 · · · ann

︸ ︷︷ ︸

x1
x2
...
xn

︸ ︷︷ ︸
=

b1
b2
...
bn

︸ ︷︷ ︸
(3.2.2)
A x b
onde A é chamada de matriz dos coeficientes, x a matriz das incógnitas e b a matriz dos termos
independentes. Algumas vezes nos referimos ao sistema linear 3.2 na forma matricial simplificada.
28
Em
El
ab
or
aç
ão
3.2. Sistemas Lineares 29
Ax = b (3.2.3)
A matriz C que obtemos ao escrevermos numa única matriz a matriz A e a matriz b é chamada de
matriz ampliada do sistema linear, escrevemos
C = [Ab]
ou
C =

a11 a12 · · · a1n b1
a21 a22 · · · a2n b2
...
...
. . .
...
...
an1 an2 · · · ann bn
 (3.2.4)
Exemplo 3.2.1
O sistema linear

x1 + 12x2 + 3x3 = 8
2x1 + 5x2 + 4x3 = 8
3x1 + 14x2 + 6x3 = 12
é representado na forma matricial por
 1 2 3
2 5 4
3 4 6
 x1
x2
x3
 =
 8
8
12

e sua matriz ampliada é
C =
 1 2 3 8
2 5 4 8
3 4 8 12

Existência e Unicidade da Solução de um Sistema Linear
Quanto à existência e número de soluções um sistema linear é classificado do seguinte modo:
Impossível. Quando não existe solução
Possível e determinado. Quando existe uma única solução
Possível e indeterminado. Quando existem infinitas soluções.
Transformações Elementares
De um modo geral, os sistemas lineares que cuja matriz dos coeficientes são matrizes triangulares, inferior
ou superior, são mais simples de resolver.
Exemplo 3.2.2
A solução (x1, x2, x3) do sistema linear
x1 + 2x2 + 3x3 = 8
x2 − x3 = −8
−7x3 = −28
(3.2.5)
pode facilmente ser encontrada, substituindo o valor de x3 na equação 3.2 para encontrar o valor de
x2; e em seguida substituindo os valores de x2 e x3 na primeira equação para calcular o valor de x1.
Dessa maneira, ao fim de dois passos a solução do sistema 3.2.2 é calculada. Verifique que a solução
encontrada é (−4, 0, 4).
Em
El
ab
or
aç
ão
30 3. Sistemas de Equações Lineares
Os métodos diretos para resolução de sistemas lineares que estudamos aqui procuram transformar
a matriz dos coeficientes em uma matriz triangular para em seguida resolver o sistema triangular e se
baseiam no fato de que dado um sistema linear Ax = b podemos, através de certas operações, obter outro
sistema linear Ax = b, com A sendo uma matriz triangular e com a mesma solução do sistema Ax = b.
Transformações de matrizes
Seja A uma matriz de ordem m× n. Definimos as transformações, ou operações, elementares nas linhas
de A do seguinte modo
(i) Troca das linhas i e k: Li ⇔ Lk;
(ii) Multiplicação de uma linha Li por um número real c diferente de zero: Li → cLi;
(iii) Substituição de uma linha Li pela soma desta linha com um múltiplo de outra linha: Li → Li+ cLk
Teorema 3.2.1 Seja C = [Ab] a matriz ampliada do sistema linear Ax = b e seja C = [Ab] a
matriz obtida a partir da matriz C através da aplicação sucessiva de um número finito de transformações
elementares,então a solução do sistema linear Ax = b é igual a solução do sistema Ax = b.
Exemplo 3.2.3
O sistema linear do Exemplo 3.2.2 foi obtido do sistema linear do Exemplo 3.2.1 a partir de
transformações elementares nas linhas de sua matriz ampliada, do seguinte modo:
C =
 1 2 3 | 8
2 5 4 | 8
3 4 6 | 12

Multiplicamos a primeira linha por −2 e somamos à segunda linha da matriz C; multiplicamos a
primeira linha por −3 e somamos à terceira linha, assim obtemos: 1 2 3 | 8
0 1 −2 | −8
0 −2 −3 | −12

Agora multiplicamos a 2◦ linha por 2 e somamos à 3◦ linha. Obtemos a matriz 1 2 3 | 8
0 1 −2 | −8
0 0 −7 | −28

que é a matriz ampliada do seguinte sistema linear triangular superior:

x1 + 2x2 + 3x3 = 8
x2 − 2x3 = −8
−7x3 = −28
3.3 Resolução de Sistemas Lineares
Sistema Triangular Inferior
A solução de um sistema triangular inferior como
a11x1 = b1
a21x1 + a22x2 = b2
...
...
an1x1 + an2x2 + · · ·+ annxn = bn
pode ser calculada efetuando-se substituições sucessivas do seguinte modo:
Em
El
ab
or
aç
ão
3.3. Resolução de Sistemas Lineares 31
x1 =
b1
a11
x2 =
b2 − a21x1
a22
x3 =
b3 − (a31x1 + a32x2)
a33
...
xn =
bn − (an1x1 + an2x2 + · · ·+ ann−1xn−1)
ann
De um modo geral, para i = 1, 2, . . ., n
xi =
bi −
i−1∑
j=1
aij
aii
Algoritmo 6: Sistema Triangular Inferior
Entrada: aij, 1 6 i 6 n e i > j
bi, 1 6 i 6 n
Saída: x1, x2, . . . , xn
início
x1 =
b1
a11
;
para i = 2, 3, . . . ,n faça
soma = 0;
para j = 1, 2, . . . , i− 1 faça
soma = soma− aij;
fim
xi =
bi − soma
aii
;
fim
fim
Sistema Triangular Superior
A solução de um sistema triangular superior

a11x1 + a12x2 + · · ·+ a1nxn = b1
a22x2 + · · ·+ a2nxn = b2
...
annxn = bn
pode ser calculada através de substituição retroativas, do seguinte modo:
Em
El
ab
or
aç
ão
32 3. Sistemas de Equações Lineares
xn =
bn
ann
xn−1 =
bn−1 − an−1nxn
an−1an−1
...
x1 =
b1 − (a12x2 + a13x3 + · · ·+ a1nxn)
a11
Generalizando, para i = n, n− 1, . . ., 1
xi =
bi −
n∑
j=i+1
aij
aii
Algoritmo 7: Substituição Retroativa
Entrada: aij, 1 6 i 6 n e i 6 j
bi, 1 6 i 6 n
Saída: x1, x2, . . . , xn
início
xn =
bn
an1
;
para i = n− 1,n− 2, . . . , 1 faça
soma = 0;
para j = i+ 1, . . . ,n faça
soma = soma− aij;
fim
xi =
bi − soma
aii
;
fim
fim
Número de Operações
Tanto o algoritmo de substituições numéricas quanto o de substituições retroativas realizam n2
2 + n
2 − 1
adições, n
2
n − n
2 multiplicações e n divisões, totalizando, cada um n2 + n− 1 operações.
3.4 Métodos Diretos
Métodos diretos são aqueles que resolvem o sistema linear Ax = b realizando um número finito de passos.
Nesta seção estudaremos o método da Eliminação Gaussiana, a Fatoração LU e a Fatoração de Cholesky.
Esses métodos realizam duas fases distintas para calcular a solução de um sistema linear. Primeiro, efetua-
se transformações elementares nas matrizes do sistema de modo a triangularizar a amtriz dos coeficientes;
e em seguida resolve-se um ou mais sistemas triangulares.
Eliminação Gaussiana
Considere um sistema linear de ordem 4 cuja matriz ampliada é dada por
Em
El
ab
or
aç
ão
3.4. Métodos Diretos 33
C =

a11 a12 a13 a14 | b1
a21 a22 a23 a24 | b2
a31 a32 a33 a34 | b3
a41 a42 a43 a44 | b4

A primeira fase da eliminação Gaussiana consiste na triangularização dessa matriz, o que é feito por
erapas.
Etapa 1. Eliminação da variável x1 nas equações 2, 3 e 4.
Para isso calcularesmo os multiplicadores m21, m31 e m41 do seguinte modo:
m21 =
a21
a11
, m31 =
a31
a11
e m41 =
a41
a11
.
Os elementos da 1◦ coluna que estão abaixo de a11 serão zerados com as seguintes transformações
elementares:
L
(1)
2 = L2 −m21L1; L
(1)
3 = L3 −m31L1 e L(1)
4 = L4 −m41L1.
Ao final dessa etapa, obteremos a matriz
C(1) =

a11 a12 a13 a14 | b1
0 a
(1)
22 a
(1)
23 a
(1)
24 | b
(1)
2
0 a
(1)
32 a
(1)
33 a
(1)
34 | b
(1)
3
0 a
(1)
42 a
(1)
43 a
(1)
44 | b
(1)
4

Etapa 2. Nesta etapa o objetivo é a eliminação dos coeficientes x2 nas linhas 3 e 4 da matriz C(1).
Os multiplicadores são m32 =
a
(1)
32
a
(1)
22
e m42 =
a
(1)
42
a
(1)
22
e as transformações elementares são:
L
(2)
3 = L
(1)
3 −m32L
(1)
2 e L(2)
4 = L
(1)
4 −m42L
(1)
2
ao final da etapa 2, teremos a matriz
C(2) =

a11 a12 a13 a14 | b1
0 a
(1)
22 a
(1)
23 a
(1)
24 | b
(1)
2
0 0 a
(2)
33 a
(2)
34 | b
(2)
3
0 0 a
(2)
43 a
(2)
44 | b
(2)
4

Por fim, temos a Etapa 3. O objetivo desta etapa é eliminar o coeficiente de x3 da linha 4 de C(2).
Nesta etapa, o multiplicador é m43 =
a43
a44
e a transformação elementar é L(3)
4 = L
(2)
4 −m43L
(2)
3 .
Dessa forma, ao final da etapa 3 teremos a matriz
C(3) =

a11 a12 a13 a14 | b1
0 a
(1)
22 a
(1)
23 a
(1)
24 | b
(1)
2
0 0 a
(2)
33 a
(2)
34 | b
(2)
3
0 0 0 a
(3)
44 | b
(3)
4

e o processo de triangularização está finalizado.
Em
El
ab
or
aç
ão
34 3. Sistemas de Equações Lineares
Os elementos a11, a
(1)
22 , a(2)
33 e a(3)
44 que aparecem na diagonal principal da matriz C(k), a cada etapa
k, são chamados de pivôs. Se em alguma dessas etapas, o pivô se anular, devemos efetuar uma troca
de linhas, antes de começar a etapa seguinte. Sempre deve-se escolher linhas abaixo da diagonal
principal para não perder a eliminação realizada na etapa anterior.
Ao final do processo de triangularização, o sitema triangular superior obtido

a11 a12 a13 a14
0 a
(1)
22 a
(1)
23 a
(1)
24
0 0 a
(2)
33 a
(2)
34
0 0 0 a
(3)
44


x1
x2
x3
x4
 =

b1
b
(1)
2
b
(2)
3
b
(3)
4

é resolvido através do Algoritmo 7 (distribuição retroativa).
Algoritmo 8: Eliminação Gaussiana sem troca de linhas
Entrada: aij, 1 6 i 6 n e 1 6 j 6 n
bi, 1 6 i 6 n
para k = 1, 2, . . . ,n− 1 faça
para i = k+ 1, . . . ,n faça
m =
aik
akk
(supondo akk 6= 0);
para j = k+ 1,k+ 2, . . . ,n faça
aij = aij −m× akj;
fim
bi = bi −m× bk;
fim
fim
Em
El
ab
or
aç
ão
4
Introdução à Resolução de Sistemas
Não Lineares
Dada uma função não linear Fξ : U→ Rn em que F = (f1, f2, f3, . . . , fn), nossoobjetivo é resolver o sistema
com n funções e n variáveis. 
f1(x1, x2, . . . , xn) = 0
f2(x1, x2, . . . , xn) = 0
...
fn(x1, x2, . . . , xn) = 0
Definição 4.0.1 O vetor gradiente de fi(x) é o vetor formado pelas derivadas parciais de fi(x).
∇fi(x) =
(
∂fi
∂x1
;
∂fi
∂x2
;
∂fn
∂xn
)
.
Definição 4.0.2 A matriz Jacobiana, ou somente o Jacobiano, é a matriz
J(x) =

∇f1(x)
∇f2(x)
∇f3(x)
...
∇fn(x)
 =

∂f1
∂x1
∂f1
∂x2
∂f1
∂x3
· · · ∂f1
∂xn
∂f2
∂x1
∂f2
∂x2
∂f2
∂x3
· · · ∂f2
∂xn
∂f3
∂x1
∂f3
∂x2
∂f3
∂x3
· · · ∂f3
∂xn
...
...
...
...
...
∂fn
∂x1
∂fn
∂x2
∂fn
∂x3
· · · ∂fn
∂xn

Exemplo 4.0.1
Seja F(x) =
{
x3
1 − 3x1x
2
2 + 1 = 0 ⇔ f1(x1, x2) = 0
3x2
1x2 − x
3
2 = 0 ⇔ f2(x1, x2) = 0 . Calcule o Jacobiano.
J(x) =
(
3x2
1 − 3x2
2 −6x1x2
6x1x2 3x1 2 − 6x2
2
)
Para obtermos a solução de sistemas não lineares, utilizaremos métodos iterativos, isto é, a partir de
um ponto inicial x(0), geramos uma sequência de vetores que convergem para a solução.
35
Em
El
ab
or
aç
ão
36 4. Introdução à Resolução de Sistemas Não Lineares
4.1 Método de Newton
Similar ao que foi feito, para se definir a função de iteração no calcular da raiz, pelo método de Newton,
em que f(x) = 0 vamos aproximar o sistema de equações lineares a partir das primeiras derivadas das
funções, então
fi(x) = f1(x
k) +∇fi(ci) + (x− xk) ; F(x) =

f1(xi) = 0
f2(xi) = 0
...
fn(xi) = 0
tomando ∇fi(xi) podemos aproximar a função fi(x) como e escrevendo na forma matricial.
F(x) ≈ Lk(x) = F(xk) + J(xk)(x− xk)
onde Lk(x) é a aproximação da função F(x).
Como F(x) = 0, a aproximação Lk(x) = 0⇒ 0 = F(xk) + J(xk)(x− xk)⇒
J(xk)(x− xk) = −F(xk)
Portanto, temos que resolver esse sistema linear com xk+1 = xk.
Para simplificar a notação podemos denotar sk = x− xk, então
J(xk)sk = −F(xk) e xk+1 = sk + xk
ou seja, resolvemos o sistema linear J(xk)sk = −F(xk) e com o resultado somamos o ponto xk, isto é,
xk+1 =sk + xk.
Observação 4.1.1 Para que seja garantida a convergência da solução, é necessário que a aproxi-
mação inicial seja próxima da solução.
Exemplo 4.1.1
Determine a solução do sistema não linear
{
x1 + x2 − 3 = 0
x2
1 + x
2
2 − 9 = 0 com x0 =
(
1
5
)
.
O Jacobiano do sistema é J(x) =

∂f1
∂x1
∂f1
∂x2
∂f2
∂x1
∂f2
∂x2

J(x) =
(
1 1
2x1 2x2
)
⇒ J(x0) =
(
1 1
2 10
)
F(x0) =
(
3
17
)
Resolvendo o sistema J(x0)s = −F(x0):(
1 1
2 10
)(
s0
1
s0
2
)
=
(
−3
−17
)
Utilizando a eliminação de Gauss (
1 1 | −3
2 1 | −17
)
Pivô é a11 = 1 ; m21 = a21
a11
= 2 ; m21( 1 1 | 0 ) = ( −2 −2 | 6 )
Em
El
ab
or
aç
ão
4.1. Método de Newton 37
(
1 1 | −3
0 8 | −11
)
⇒ 8s(0)
2 = −11 ⇒
s1 + s2 = −3 ⇒
{
s2 = −11/8
s1 = −13/8
Logo
x1 = x0 + s0 =
(
1
5
)
+
(
−13/8
−11/8
)
=
(
−0, 625
3, 625
)
⇒
x1 =
(
−5/8
29/8
)
2◦ Iteração:
J(x1) =
(
1 1
−5/4 29/4
)
; F(x1) =
(
0
290
64
)
; J(x1)s = −F(x1)(
1 1
−5/4 29/4
)(
s
(1)
1
s
(1)
2
)
=
(
0
−290/64
)
Eliminação de Gauss. Pivô a11 = 1, m21 = a21
a11
= −5
4 =
(
5/4 5/4 | 0
)
(
1 1 | 0
−5/4 29/4 | −290/64
)
=
(
1 1 | 0
0 34/4 | −290/64
)
⇒
34
4
s2 =
290
64
s2 = −
145
(16× 17)
≈ −0, 53⇒
s1 + s2 = 0⇒
s1 =
145
(16× 17)
≈ −0, 53
x2 = s(1) + x(1) =
(
+0, 53
−0, 53
)
+
(
−5/8
29/8
)
⇒ x0 =
(
0, 09
3, 098
)
3◦ Iteração:
J(x2) =
(
1 1
0, 18 6, 18
)
; F(x2) =
(
0, 18
0, 56
)
; J(x2)s = −F(x2)
Logo,
(
1 1 | −0, 18
0, 18 6, 18 | −0, 56
)
.
Eliminação de Gauss, pivô a11 = 1, m21 = a21
a11
= 0, 18 ; m21
(
−0, 18 −0, 18 | 0, 01
)
(
1 1 | −0, 18
0 6 | −0, 53
)
⇒ 6s(2)
2 = 0, 53 ⇒ s
(2)
2 = −0, 09
s1 + s2 = 0, 18 ⇒ s1 = −0, 09
x(3) = s2 + x2 =
(
−0, 09
−0, 09
)
+
(
0, 09
3, 09
)
=
(
0, 00
3, 00
)
Na 3◦ iteração, conseguimos encontrar a solução, com erro na segunda casa decimal!
Observação 4.1.2 O Método de Newton converge bem, mas é trabalhoso resolver o sistema linear
a cada iteração.
Em
El
ab
or
aç
ão
38 4. Introdução à Resolução de Sistemas Não Lineares
Observação 4.1.3 Se para resolvermos o sistema linear utilizamos algum método iterativo, o método
recebe o nome de Método de Newton Iterativo.
4.2 Método de Newton Modificado
Para simplificar a solução de um sistema linear, para cada iteração do método de Newton, é feita uma
modificação
J(x0)s = −F(xk)
fixando o Jacobiano J(x0)m da iteração inicial, ao invés de J(xk) para cada ponto. Dessa forma, a matriz
Jacobiana é avaliada somente uma vez. Esse método é denominado de Método de Newton Modificado.
Exemplo 4.2.1
Aplicando o método de Newton Modificado na solução do exemplo anterior, até a 3◦ iteração, com{
x1 + x2 − 3 = 0
x2
1 + x
2
2 − 9 = 0 e x0 =
(
1
5
)
, J(x) =
(
1 1
2x1 2x2
)
A primeira iteração é feita da mesma forma: J(x0) =
(
1 1
2 10
)
e F(x0) =
(
3
17
)
, então(
1 1
2 10
)(
s
(0)
1
s
(0)
2
)(
3
17
)
⇒ s =
(
−13/8
−11/8
)
Com x1 =
(
−5/8
29/8
)
=
(
−0, 625
3, 625
)
2◦ iteração:
f(x1) =
(
0
290/64
)
⇒ J(x0)s = −F(x1)⇒
(
1 1
2 10
)(
s
(1)
1
s
(1)
2
)
=
(
0
−4, 53
)
Eliminação de Gauss, pivô a11 = 1⇒ m21 = a21
a11
= 2, m21( −2 −2 | 0 )(
1 1 | 0
0 8 | −290/64
)
⇒ s2 = −
290
64× 8
= −0, 57⇒ s1 + s2 = 0⇒ s1 = +0, 57
x2 = s1 + x1 =
(
0, 57 −0, 57
)
+
(
−6, 625
3, 625
)
=
(
−6, 06
3, 06
)
3◦ iteração:
f(x2) =
(
0
0, 37
)
⇒ J(x0)s = −F(x2)(
1 1
2 10
)(
s2
1
s2
1
)
=
(
0
−0, 37
)
(
1 1
0 8
)(
s2
1
s2
2
)
=
(
0
−0, 37
)
Eliminação de Gauss, m21 = 2, m21( −2 −2 | 0 )
8s2 = −0, 37⇒ s2
2 = −0, 05 ; s1 + s2 = 0⇒ s1 = −s2 = 0, 05
x3 = s2 + x2 =
(
0, 05
−0, 05
)
⇒ x3 =
(
−0, 01
3, 01
)
Em
El
ab
or
aç
ão
4.2. Método de Newton Modificado 39
O Método de Newton Modificado não converge tão rápido quanto o Método de Newton, mas em
compensação só é necessário calcular uma matriz Jacobiana.
Uma outra modificação feita no método de Newton é aproximar a matriz Jacobiana J(xk) por uma
matriz B(k), muito útil em situações onde calcular a derivada pode ser custoso computacionalmente.
Esses métodos são chamados de Quase-Newton.
Para sistemas nãolineares 2× 2, podemos obter uma expressão mais simples para a inversa da matriz
Jacobiana, já que A =
(
a b
c d
)
, então A−1 = 1
detA
(
d −b
−c a
)
= 1
ad−bc
(
d −b
−c a
)
.
Uma vez calculado o Jacobiano J =
(
∂xf1 ∂yf1
∂xf2 ∂yf2
)
Jsk = −f(xk)⇒ sk = −J−1F(xk)
Logo,
xk+1 = sk + xk = xk − J−1F(x)⇒
xk+1 = xk −
1
∂xf1∂yf2 − ∂yf1∂xf2
(
∂xf1 ∂yf1
∂xf2 ∂yf2
)(
F1(x
k)
F2(xk)
)
xk+1 =
(
xk1
xk2
)
−
1
∂xf1∂yf2 − ∂yf1∂xf2
(
F1(x
k)∂xf1 + F2(x
k)∂yf1
F1(x
k∂xf2 + F2(x
k)∂yf2
)
Em
El
ab
or
aç
ão
5
Interpolação
Interpolar uma função f(x) consiste em aproximar essa função por uma outra, g(x), escolhida entre um
classe de funções definidas previamente. Assim, a análise de f(x) é substituida pela análise de g(x).
x −20 −12 −10 −5 0 2 8 16
f(x) 54 13 12 1 −1 12 14 26
Qual seria o valor de x = −3, x = −2, x = 5?
f(x) ≈ g(x) e a partir do cálculo de

g(−3) ≈ f(−3)
g(−2) ≈ f(−2)
g(5) ≈ f(5)
Sejam (m+ 1) pontos distintos, da forma

(x0, f(x0
0))
(x1, f(x1)
...
(xn, f(xn))
esses pontos são chamados de nós da inter-
polação.
Uma forma de se obter a interpolação é requisitando que g(x) coincida com f(x) nos pontos conhecidos
g(x0) = f(x0)
g(x1) = f(x1)
g(x2) = f(x2)
...
...
g(xn) = f(xn)
FIGURA
Observação 5.0.1 Existem outras formas para se obter critério para interpolação, como por exem-
plo, a série de Taylor.
Ao longo do curso, utilizaremos g(x) como uma função polinomial, mas vale a pena solicitar que a
função g(x) pode ser funções de Hermite, ou funções trigonométricas.
5.1 Interpolação Polinomial
Sejam (x0, f(x0), (x1, f(x1), . . . , (xn, f(xn))m+1 pontos, então efetuar uma interpolação polinomial consiste
em aproximar f(x) por um polinômio de grau menor ou igual a n, com p(x) = a0+a1x+a2x
2+· · ·+anxn,
onde queremos obter os coeficientes ak,k = 0, 1, . . . ,m, então
40
Em
El
ab
or
aç
ão
5.2. Formas para se obter o polinômio interpolador Pn(x) 41

f(xk) = p(xk), para k = 0, 1, . . . ,n logo
f(x0) = a0 + a1x0 + a2x
2
0 + · · ·+ anx2
n
f(x1) = a0 + a1x1 + a2x
2
1 + · · ·+ anx2
1
f(x2) = a0 + a1x2 + a2x
2
2 + · · ·+ anx2
2
...
f(xn) = a0 + a1xm + a2x
2
m + · · ·+ anx2
n
um sistema com n+ 1 equações e n+ 1 variáveis.
Observação 5.1.1 Se aumentamos o grau do polinômio, não temos como obter uma solução para
o sistema acima.
Representando na forma matricial
f(x0)
f(x1)
f(x2)
...
f(xn)
 =

1 x0 x2
0 · · · xn0
1 x1 x2
1 · · · xn1
1 x2 x2
2 · · · xn2
...
...
...
. . .
...
1 xn x2
n · · · xnn

︸ ︷︷ ︸

a0
a1
a2
...
an

A matriz dos coeficientes
Então, determinar os coeficientes (a0,a1, . . . ,an) é equivalente a resolver o sistema acima. A matriz
A é uma matriz de Vondermonde, e desde que xi 6= xj, ∀ i, j = 0, 1, 2, . . . ,n temos que det(A) 6= 0. Com
isso, podemos afirmar que
Teorema 5.1.1 Existe um único polinômio pn(x), de grau 6 n, tal que pn(xk) = f(xk), k = 0, 1, 2, . . . ,n,
desde que xi 6= xj ∀ i, j.
5.2 Formas para se obter o polinômio interpolador Pn(x)
Verem os três formas equivalentes de se obter o polinômio interpolador que são:

·Resolução do Sistema Linear;
·Método de Lagrange;
·Método de Newton.
Observação 5.2.1 Os três métodos conduzem ao mesmo polinômio interpolador.
5.3 Resolução do Sistema Linear
Como o próprio título indica, vamos obter o polinômio interpolador, resolvendo o sistema linear.
Exemplo 5.3.1
A partir dos dados abaixo, responda:
x −1 0 2
f(x) 4 1 −1
a) Qual é o maior grau do polinômio interpolador?
b) Determine o polinômio interpolador;
c) Estime o valor da função f(x), para x = 1.
Em
El
ab
or
aç
ão
42 5. Interpolação
a) Com 3 pontos distintos, o grau de polinômio é 2
p(x) = a0 + a1x+ a2x
2
b) Para
{
x0 = −1
f(x0) = 4 ;
{
x1 = 0
f(x1) = 1 ;
{
x2 = 2
f(x2) = −1 , como f(xk) = p(xk)
4 = a0 − a1 + a2 ⇒
{
3 = −a1 + a2
−2 = 2a1 + 4a2
1 = a0
−1 = a0 + 2a1 + 4a2{
6 = −2a1 + 2a2
−2 = 2a1 + 4a2
⇒ 4 = 6a2 ⇒ a2 =
4
6
=
23
a1 = −3 + a2 = −3 +
2
3
= −
7
3

a0 = 1
a1 = −7/3
a2 = 2/3
Logo, p2(x) = 1 −
7
3
x+
2
3
x2
c) Para x = 1 temos que p2(1) = 1 − 7
3 + 2
3 = 3−7+2
3 = −2
3 . Então f(1) ≈ p2(1) = −2
3 então
f(1) ≈ −2
3 .
5.4 Forma de Lagrange
Sejam (x0, f(x0)), (x1, f(x1)), . . ., (xn, f(xn)), n + 1 pontos distintos. A forma de Lagrange consiste em
escrever o polinômio interpolador de grau n, de forma
pn(x) = f(x0)L0(x) + f(x1)L1(x) + f(x2)L2(x) + · · ·+ f(xn)Ln(x)
onde os polinômios Lk(x) são de grau n. Para cada xi, queremos que p(xi) = f(xi)
p(xi) = f(x0)L0(xi) + f(x1)L1(xi) + · · ·+ L(xn)Ln(xi) = f(xi), ∀ i = 0, 1, 2, . . . ,n.
Uma forma de se obter essa condição é impor que Lk(xi) =
{
0, se k 6= i
1, se k = i
e definirmos
Lk(x) =
(x− x0)(x− x1)(x− x2) . . . (x− xk−1)(x− xk+1) . . . (x− xn)
(xk − x0)(xk − x1) . . . (xk − xk−1)(xk − xk+1) . . . (xk − xn)
um produto com todos os pontos, menos o com índice k
Lk(x) =
n∏
i=0
i 6=k
(x− xi)
(xk − xi)
, então
 ·Lk(xk) =
(xk − x0)(xk − x1) . . . (xk − xk−1)(xk − xk+1) . . . (xk − xn)
(xk − x0)(xk − x1) . . . (xk − xk−1)(xk − xk+1) . . . (xk − xn)
= 1
·Lk(xi) = 0, se i 6= k.
Portanto, o polinômio interpolador a partir da forma de Lagrange é
pn(x) =
n∑
k=0
f(xk)Lk(x), com Lk(x) =
n∏
i=0
(x− xi)
(xk − xi)
Em
El
ab
or
aç
ão
5.4. Forma de Lagrange 43
pn(x) =
n∑
k=0
[
n∏
i=0
f(xk)(x− xi)
(xk − xi)
]
Exemplo 5.4.1 (Interpolação Linear pela Forma de Lagrange)
Sejam (x0, f(x0)) e (x1, f(x1)) dois pontos distintos, então o polinômio interpolador pn(x) deve ser
de grau 1.
p1(x) = f(x0)L0(x) + f(x1)L1(x)
L0(x) =
1∏
i=0
i 6=0
(x− xi)
(x0 − xi)
=
(x− x1)
(x0 − x1)
L1(x) =
1∏
i=0
i 6=o
(x− x1)
(x1 − xi)
=
(x− x0)
(x1 − x0)
Logo
p1(x) = f(x0)
(x− x1)
x0 − x1
+ f(x1)
(x− x0)
x1 − x0
⇒
p1(x) = −
f(x0)(x− x1) + f(x1)(x− x0)
x1 − x0
⇒
p1(x) =
(x1 − x)f(x0) + (x− x0)f(x1)
x1 − x0
que é a equação da reta que passa pelos pontos{
(x0, f(x0))
(x1, f(x1))
Exemplo 5.4.2
Seja
x −1 0 2
f(x) 4 1 −1
determine o polinômio interpolador pela forma de Lagrange.
Como temos 3 pontos, pn(x) é de grau 2.
p2(x) = f(x0)L0(x) + f(x1)L1(x) + f(x2)L2(x)
com f(x0) = 4 ; f(x1) = 1 e f(x2) = −1.
Os polinômios Li(x) são
• L0(x) =
2∏
i=0
i 6=0
(x− xi)
(x0 − xi)
=
(x− x1)(x− x2)
(x0 − x1)(x0 − x2)
=
(x− 0)(x− 2)
(−1 − 0)(−1 − 2)
L0(x) =
x(x− 2)
3
=
x2 − 2x
3
• L1(x) =
2∏
i=0
i 6=0
(x− xi)
(x1 − xi)
=
(x− x0)(x− x2)
(x1 − x0)(x1 − x2)
=
(x+ 1)(x− 2)
(0 + (−1))(0 − 2)
L1(x) =
(x+ 1)(x− 2)
−2
=
x2 − x− 2
−2
Em
El
ab
or
aç
ão
44 5. Interpolação
• L2(x) =
2∏
i=0
i 6=0
(x− xi)
(x2 − xi)
=
(x− x0)(x− x1)
(x2 − x0)(x2 − x1)
=
(x+ 1)(x)
(2 − (−1))(2 − 0)
L2(x) =
x2 + x
6
Logo,
p2(x) = 4L0(x) + 1L1(x) − 1L2(x)⇒
p2(x) = 4
(x2 − 2x)
3
+
x2 − x− 2
−2
+ (−1)
(x2 + x)
6
⇒
p2(x) =
8x2 − 16x− 3x2 + 3x+ 6 − x2 − x
6
⇒
p2(x) =
4x2 − 14x+ 6
6
⇒
p2(x) =
2
3
x2 −
7
3
x+ 1
Igual ao polinômio obtido resolvendo o sistema linear!
5.5 Forma de Newton
A partir de n+ 1 pontos (x0, f(x0)), (x1, f(x1)), . . ., (xn, f(xn)) o polinômio interpolador pn(x) será
pn(x) = d0 + d1(x− x0) + d2(x− x0)(x− x1) + · · ·+ dn(x− x0)(x− x1)(x− x2) . . . (x− xn−1)
onde di são os operadores de diferenças divididas que veremos a seguir.
Operador de Diferenças Divididas
Seja f(x) uma função tabelada em n + 1 pontos distintos: x0, x1, . . ., xn. Definiremos o operador de
diferenças divididas como:
• Ordem Zero:
f[x0] = f(x0)
• Ordem 1:
f[x0, x1] =
f[x1] − f[x0]
x1 − x0
=
f(x1)f(x0)
x1 − x0
• Ordem 2:
f[x0, x1, x2] =
f[x1, x2] − f[x0, x1]
x2 − x0
=
f(x2)−f(x1)
x2−x1
−
f(x1)−f(x0)
x1−x0
x2 − x0
• Ordem 3:
f[x0, x1, x2, x3] =
f[x1, x2, x3] − f[x0, x1, x2]
x3 − x0
=
f[x2, x3] − f[x1, x2]
(x3 − x1)(x3 − x0)
=
f[x1, x2] − f[x0, x1]
(x2 − x0)(x3 − x0)
• Ordem 4:
f[x1, x2, x3, x4] =
f[x0, x1, x2, x3, x4] − f[x0, x1, x2, x3]
x4 − x0
Em
El
ab
or
aç
ão
5.5. Forma de Newton 45
• Ordem n:
f[x0, x1, x2, . . . , xn] =
f[x1, x2, . . . , xn] − f[x0, x1, . . . , xn−1]
xn − x0
Para calcular as diferenças divididas utiliza-se a tabela para facilitar os cálculos.
x Ordem 0 Ordem 1 Ordem 2 Ordem n
x0 f[x0] f[x0, x1] f[x0, x1, x2]
f[x0, x1, . . . , xn]
x1 f[x1] f[x1, x2] f[x1, x2, x2]
x2 f[x2] f[x2, x3] f[x2, x3, x4]
x3 f[x3] f[x3, x4] f[x3, x4, x5]
...
...
...
...
xn f[xn] f[xn−1, xn] f[xn−2, xn−1, xn]
Observação 5.5.1 A diferença dividida são simétricas, ou seja, qualquer troca nos valores de xi
não altera o resultado, isto é,
f[x0, x1] =
f[x1] − f[x0]
x1 − x0
=
f[x0] − f[x1]
x0 − x1
= f[x1, x0]
f[x0, x1, x2] = f[x2, x0, x1] = f[x1, x0, x2]
= f[x0, x2, x1] = f[x1, x2, x0] = f[x2, x1, x0]
válido para qualquer combinação entre xi.
Exemplo 5.5.1
Seja
x −1 0 1 2 3
f(x) 1 1 0 −1 −2
x0 x1 x2 x3 x4
então as diferenças divididas são:
• Ordem 0:
f[x0] = 1 ; f[x1] = 1 ; f[x2] ; f[x3] = −1 ; f[x4] = 0
x0 = −1 ; x1 = 0 ; x2 = 1 ; x3 = 2 ; x4 = 3
• Ordem 1:
f[x0, x1] =
f[x1] − f[x0]
x1 − x0
=
1 − 1
0 − (−1)
= 0
f[x1, x2] =
f[x2] − f[x1]
x2 − x1
=
0 − 1
1 − 0
= −1
f[x2, x3] =
f[x3] − f[x2]
x3 − x2
=
−1 − 0
2 − 1
= −1
Em
El
ab
or
aç
ão
46 5. Interpolação
f[x3, x4] =
f[x4] − f[x3]
x4 − x3
=
−2 − (−1)
3 − 2
= −1
• Ordem 2:
f[x0, x1, x2] =
f[x1, x2] − f[x0, x1]
x2 − x0
=
−1 − 0
1 − (−1)
= −
1
2
f[x1, x2, x3] =
f[x2, x3] − f[x1, x2]
x3 − x1
=
−1 − (−1)
2 − 0
= 0
f[x2, x3, x4] =
f[x3, x4] − f[x2, x3]
x4 − x2
=
−1 − (−1)
3 − 1
= 0
• Ordem 3:
f[x0, x1, x2, x3] =
f[x1, x2, x3] − f[x0, x1, x2]
x3 − x0
=
0 −
(
−1
2
)
2 − (−1)
=
1
6
f[x1, x2, x3, x4] =
f[x2, x3, x4] − f[x1, x2, x3]
x4 − x1
=
0 − 0
3 − 0
= 0
• Ordem 4:
f[x0, x1, x2, x3, x4] =
f[x1, x2, x3, x4] − f[x0, x1], x2, x3]
x4 − x0
=
0 −
(1
6
)
3 − (−1)
= −
1
24
A tabela de diferenças divididas é:
x Ordem 0 Ordem 1 Ordem 2 Ordem 3 Ordem 4
x0 = −1 f[x0] = 1
f[x0, x1] = 0
−1
2
− 1
24
x1 = 0 f[x1] = 1
−1
0
1
6
x2 = 1 f[x2] = 0
−1
0
0
x3 = 2 f[x3] = −1
−1
x4 = 3 f[x4] = −2
5.6 Forma de Newton para pn(x)
Seja f(x) uma função contínua e derivável em [a,b], com n+ 1 pontos
a = x0 < x1 < x2 < . . . < xn = b
ordenados de forma crescente. Recursivamente, temos que
• Para n = 0, o polinômio de grau 0 que interpola f(x) em x = x0, então
p0(x) = f(x0) = f[x0] Diferença dividida de ordem zero.
então para qualquer x ∈ [a,b], x 6= x0
f[x0, x] =
f[x] − f[x0]
x− x0
=
f(x) − f(x0)
x− x0
⇒
f(x) = f(x0)︸ ︷︷ ︸+ [x− x0] · f[x0, x]︸ ︷︷ ︸⇒
p(x0) E0(x) é o erro cometido ao se
aproximar f(x) por p0(x), dado por
Em
El
ab
or
aç
ão
5.6. Forma de Newton para pn(x) 47
E0(x) = f(x) − p(x0) com p0(x) = f(x0).
• Para n = 1, o polinômio de grau 1 que interpola f(x) em x0 e x1, temos que
f[x0, x1, x] = f[x1, x0, x] =
f[x0, x] − f[x1, x0]
x− x1
=
f[x]−f[x0]
(x−x0)
− f[x1, x0]
x− x1
=
f(x) − f[x0] − f[x1, x0](x− x0)
(x− x1)(x− x0)
⇒
f(x) = f[x0] + (x− x0)f[x1, x0]︸ ︷︷ ︸+(x− x0)(x− x1)f[x0, x1, x]︸ ︷︷ ︸
p1(x) polinômio de grau 1 E1(x) é o erro ao se interpolar f(x) por p1(x)
então p1(x) = f(x0)︸ ︷︷ ︸+(x− x0)f[x0, x1]︸ ︷︷ ︸
p0(x) q1(x)
é o polinômio de grau 1 que passa por (x0, f(x0)) e (x1, f(x1) já que
 p1(x0) = f(x0) + (x0 − x0)f[x0, x1] = f(x0)
p1(x1) = f(x0) + (x1 − x0)
f[x1] − f[x0]
x1 − x0
= f(x1)
• Para n = 2, o polinômio de grau 2 que interpola f(x) em x0, x1 e x2, com
f[x0, x1, x2, x] = f[x2, x1, x0, x] =
f[x1, x0, x] − f[x2, x1, x0]
x− x2
=
f[x0,x]−f[x1,x0]
(x−x1)
− f[x2, x1, x0]
x− x2
=
f(x) − f(x0) − f[x1, x0](x− x0) − (x− x0)(x− x1)f[x2, x1, x0]
(x− x0)(x− x1)(x− x0)
⇒
f(x) = f[x0] + f[x1, x0](x− x0) + f[x2, x1, x0](x− x0)(x− x1)︸ ︷︷ ︸+ f[x0, x1, x2, x](x− x0)(x− x1)(x− x2)︸ ︷︷ ︸
p2(x) E2(x) é o erro ao se interpolar f(x) por p2(x)
p2(x) = f[x0] + (x− x0)f[x1, x0]︸ ︷︷ ︸+(x− x0)(x− x1)f[x2, x1, x0]︸ ︷︷ ︸
p1(x) p2(x)
Então, de modo geral o polinômio
pk(x) = pk−1(x) + qk(x)
Aplicando sucessivamente o mesmo raciocínio para um polinômio de grau n, que interpola f(x) em
x0, x1, x2, . . . , xn obtemos que
pn(x) = f(x0) + (x− x0)f[x0, x1] + (x− x0)(x− x1)f[x0, x1, x2]+
(x− x0)(x− x1)(x− x2)f[x0,x1, x2, x3] + · · ·+ (x− x0)(x− x1) . . . (x− xn−1)f[x0, x1, x2, . . . , xn],
onde o erro na aproximação de f(x) por pn(x) é
En = (x− x0)(x− x1) . . . (x− xn)f[x0, x1, x2, . . . , xn, x]
Em
El
ab
or
aç
ão
48 5. Interpolação
Portanto, pn(x) construído pela forma de Newton passa pelos n+1 pontos (x0, f(x0)), (x1, f(x1)), . . . , (xn, f(xn))
já que
f(x) = pn(x) + En(x), então
f(xk) = pn(xk) + En(xk) e como En(xk) = 0, ∀ k = 0, 1, 2, . . . ,n
f(xk) = pn(xk)
Observação 5.6.1 Ao utilizarmos a forma de Newton para o polinômio interpolador, obtemos uma
expressão para o erro obtido na interpolação.
Exemplo 5.6.1
Usando a forma de Newton, obtenha o polinômio p2(x) que interpola f(x) nos pontos
x −1 0 2
f(x) 4 1 −1
o polinômio p0(x) pela forma de Newton é:
p2(x) = f[x0] + (x− x0)f[x0, x1] + (x− x0)(x− x1)f[x0, x1, x2]
Montando a tabela de diferenças divididas:
x Ordem 0 Ordem 1 Ordem 2
x0 = −1 f(x0) = 4
f[x0, x1] =
1−4
0−(−1) = −3
f[x0, x1, x2] =
−1−(−3)
2−(−1) = 2
3x1 = 0 f(x1) = 1
f[x1, x2] =
−1−1
2−0 = −1
x2 = 2 f(x2) = −1
Logo,
p2(x) = 4 + (x− x0)(−3) + (x− x0)(x− x1)
2
3
p2(x) = 4 − 3(x+ 1) +
2
3
(x+ 1)x
= 4 − 3x− 3 +
2
3
x2 +
2
3
x2 ⇒
p2(x) = 1 −
7
3
x+
2
3
x2
o mesmo polinômio obtido pela forma de Lagrange.
5.7 Estudo do Erro na Interpolação
Sempre que aproximamos uma função f(x), com n+ 1 pontos, por um polinômio interpolador de grau n
comete-se um erro
En = f(x) − pn(x)
Em
El
ab
or
aç
ão
5.8. Interpolação Inversa 49
FIGURA
Fizendo os n+ 1 pontos temos um único polinômio pn(x), mas temos uma infinidade de funções que
passam pelos pontos.
Pela forma de Newton obtivemos uma expressão para o erro, dado a partir das diferenças divididas
En = (x− x0)(x− x1) . . . (x− xn)f[x0, x1, . . . , xn, x]
mas para calcular o erro devemos conhecer a função f(x), o que nem sempre é possível. Um outro resultado,
na mesma perspectiva é que:
Teorema 5.7.1 Sejam a = x0 < x1 < . . . xn = b n+ 1 pontos e f(x) contínua e derivável em [a,b], com
polinômio interpolador pn(x), então
En = f(x) − pn(x) = (x− x0)(x− x1) . . . (x− xn)
fn+1(ξ)
(n+ 1)!
com ξ ∈ [a,b] um ponto do intervalo.
Observação 5.7.1 o Teorema 5.7.1 só garante a existência de ξ, não como obtê-lo.
Comparando o erro pela forma de Newton, com o Teorema anterior, obtemos uma relação entre as
diferenças divididas de ordem n+ 1 com a derivada de ordem n+ 1 de f(x)
f[x0, x1, x2, . . . , xn, x] =
f(n+1)(ξ)
(n+ 1)!
, com ξ ∈ [a,b]
Quando a função f(x) é dada na forma de tabela, o valor do erro só pode ser estimado pois não
podemos calcular sua diferença dividida, nem as suas derivadas.
Uma estimativa para o erro consiste em utilizar o maior valor em método da tabela das diferenças
divididas, como aproximação para as derivadas de ordem n+ 1.
Então
|En(x) ≈ (x− x0)(x− x1) . . . (x− xn)| max{diferenças divididas}
Exemplo 5.7.1
No exemplo anterior, uma estimativa para o erro seria
|E2(x)| = (x+ 1) · (x) · (x− 2) · 4︸︷︷︸
f[x0], o maior valor da
tabela de diferenças di-
vididas
5.8 Interpolação Inversa
O problema da interpolação inversa consiste em, a partir de um conjunto de dados:
x x0 x1 x2 · · · xn
f(x) f(x0) f(x1) f(x2) · · · f(xn)
Dado y ∈ [f(x0); f(xn)] obter x tal que f(x) = y.
Podemos resolver este problema de duas formas:
i) Obter o polinômio pn(x) que interpola x0, x1, . . . , xn e determinar o x tal que p(x) = y.
Em
El
ab
or
aç
ão
50 5. Interpolação
Observação 5.8.1 Seja y ∈ [f(xi); f(xi+1)] podemos tomar pn como um polinômio de grau 1
que passa por (xi, f(xi)) e (xi+1, f(xi+1)) ou qualquer outro polinômio.
ii) Se f(x) for crescente ou decrescente no intervalo em análise, podemos supor que f(x) é inversível
nesta região e consequentemente interpolar a função inversa.
f(x) f(x0) f(x1) f(x2) · · · f(xn)
f−1(f(x)) = x x0 x1 x2 · · · xn
Exemplo 5.8.1
A prtir da tabela
x 0, 5 0, 6 0, 7 0, 8 0, 9 1, 0
f(x) 1, 65 1, 82 2, 01 2, 23 2, 46 2, 72
obtenha uma estimativa para x, em que f(x) = 2, 1 utilizando
i) Um polinômio de grau 1;
ii) A interpolação inversa. Forneça uma estimativa para o erro.
i) 2, 1 ∈ (2, 01; 2, 23) com
{
x2 = 0, 7
x3 = 0, 8
o polinômio linear pode ser obtido por qualquer método.
{
p1(x) = f(x0) + (x− x0)f[x0, x1] Newton
p1(x) = f(x0)L0(x) + f(x1)L1(x) Lagrange
⇒ p1(x) = 2, 01 + (x− 0, 7)
f[x1] − f[x0]
x1 − x0
= 0, 7 + (x− 0, 7) · 2, 2 = 2, 2x+ 0, 147
⇒ p1(x) = 2, 01
(x− x1)
x0 − x1
+ 2, 23
(x− x0)
x1 − x0
=
2, 01
−0, 1
(x− 0, 8) + 2, 23
(x− 0, 7)
0, 1
=
(2, 23 − 2, 01)
0, 1
x+
2, 23(x− 0, 7)
0, 1
−
2, 23 · 0, 7
0, 1
= 2, 2x+ 0, 47
Portanto, o x tal que p1(x) = 2, 1 é 2, 1 = 2, 2x+ 0, 47
x =
2, 1 − 0, 47
2, 2
= 0, 741
Por esse método não temos como avaliar o erro ao calcularmos x = 0, 741.
ii) Efetuando a interpolção inversa, já que a função é crescente no intervalo [0, 5; 1, 0], logo
y 1, 65 1, 82 2, 01 2, 23 2, 46 2, 72
g(y) 0, 5 0, 6 0, 7 0, 8 0, 9 1, 0
Em
El
ab
or
aç
ão
5.9. Fenômeno de Runge 51
Obtemos p5(x) pelo Método de Newton, a fim de estimar o erro no cálculo do processo.
x f(x) Ordem 1 Ordem 2 Ordem 3 Ordem 4 Ordem 5
x0 = 1, 65 f(x0) = 0, 5
0, 59
−0, 17
−0, 05
0, 37
−0, 68
x1 = 1, 82 f(x1) = 0, 6
0, 53
−0, 20
0, 25
−0, 36
x2 = 2, 01 f(x2) = 0, 7
0, 45
−0, 04
−0, 08
x3 = 2, 23 f(x3) = 0, 8
0, 43
−0, 1x4 = 2, 46 f(x4) = 0, 9
0, 38
x5 = 2, 72 f(x5) = 1, 0
p5(x) = g(y0) + (y− y0)g[y0,y1] + (y− y0)(y− y1)g[y0,y1,y2]
+ (y− y0)(y− y1)(y− y2)g[y0,y1,y2,y3]
+ (y− y0)(y− y1)(y− y2)(y− y3)g[y0,y1,y2,y3,y4]
+ (y− y0)(y− y1)(y− y2)(y− y3)(y− y4)g[y0,y1,y2,y3,y4,y5]
Logo, para y = 2, 1, temos que
p5(x) = 0, 5 + (y− 1, 65)0, 59 + (y− 1, 65)(y− 1, 82)(−0, 17)
+ (y− 1, 65)(y− 1, 82)(y− 2, 01)(−0, 05)
+ (y− 1, 65)(y− 1, 82)(y− 2, 01)(y− 2, 23)(0, 37)
+ (y− 1, 65)(y− 1, 82)(y− 2, 01)(y− 2, 23)(y− 2, 43)(−0, 68)
Logo
p(2, 1) = 0, 742636735
com uma diferença na 3a casa decimal.
Pelo método de Newton somos capazes de estimar o erro, que é dado por
E(y) ≈ |(y−y0)(y−y1)(y−y2)(y−y3)(y−y4)(y−y5)| (máximo {|diferençasdivididas|}) = 1
E(2, 1) = 3, 2904× 10−4
≈ 4× 10−4 = 0, 0004
5.9 Fenômeno de Runge
Ao efetuarmos uma interpolação de uma função f(x) num intervalo [x0, xn] pode decorrer um erro ao se
calcular os valores nas bordas do intervalo.
Mesmo aumentando o grau do polinômio pn, este erro persiste FIGURA e consiste em estarmos
utilizando polinômios para aproximar uma função.
Em
El
ab
or
aç
ão
52 5. Interpolação
Este comportamento de pn(x) nas bordas aumenta cada vez mais com o grau de pn(x) e é conhecido
como fenômeno de Runge.
Este fenômeno pode ser contornado de três formas:
i) Tomando como base outras funções, diferentes dos polinômios (o que está fora do nosso curso).
{
Funções racionais
Trigonométricas
ii) Escolher os pontos de forma diferente, para tornar o erro mais homogêneo. Como exemplo temos os
pontos de Chebyshev:
xi =
x0 + xn
2
+
xn − x0
2
· cos
(
2i+ 1
2π+ 2
π
)
ii) Usa funç eos Spline, que são um conjunto de funções com um grau menor, e impor condições para
que as aproximações sejam contínuas e com derivadas contínuas.
5.10 Spline
Seja f(x) uma função com
x x0 x1 x2 . . . xn
f(x) f(x0) f(x1) f(x2) . . . f(xn)
construiremos uma função Sp(x) que interpola f(x) num grupo de pontos.
Definição 5.10.1 Seja f(x) uma função com n+ 1 pontos (x0, f(x0), . . . , (xn, f(xn)), a função Sp(x)
é denominada de Spline interpolante se:
i) Para cada [xi, xi+1] Sp(x) é um polinômio de grau p · sp(x);
ii) Sp(x) é contínua e tem derivadas contínuas até a ordem (p− 1) em [x0, xn];
iii) Sp(x) = f(xi) ∀ i ∈ {0, 1, 2, . . . ,n}
FIGURA(caso de spline linear)
Sp =

s1(x), se x ∈ [x0, x1]
s2(x), se x ∈ [x1, x2]
s3(x), se x ∈ [x2, x3]
...
sn(x), se x ∈ [xn−1, xn]
Dividimos os pontos em grupos e calculamos os polinômios de grau k nessas regiões pedindo que{
Sk(x) seja contínua e com derivadas
contínuas até ordem (k− 1).
Em
El
ab
or
aç
ão
5.11. Spline Linear Interpolante 53
5.11 Spline Linear Interpolante
Consideraremos o caso do Spline ser uma função linear, ou seja, de grau 1 agrupandoos pontos em pares.
Então
S1(x) =
si(x) = f(xi − 1)
(x− xi)
xi−1 − xi
+ f(xi)
(x− xi−1)
xi − xi−1︸ ︷︷ ︸, ∀ x ∈ [xi−1, xi]
Polinômio de Lagrande que passa por xi−1 e xi
Portanto, agrupamos os pontos em pares de calculamos a reta que passa por xi−1 e xi para cada
intervalo.
FIGURA
Exemplo 5.11.1
Determine a função spline linear que interpola a função
x0 x1 x2 x3
x 1 2 5 7
f(x) 1 2 3 2, 5
Basta determinar as 3 retas que passam por
(x0, f(x0)); (x1, f(x1)) 7→ s1(x)
(x1, f(x1)); (x2, f(x2)) 7→ s2(x)
(x2, f(x2)); (x3, f(x3)) 7→ s3(x)
(x3, f(x3)); (x4, f(x4)) 7→ s4(x)
dados pela Forma de Lagrange de grau 1.
• si(x) = f(xi)
x− xi
xi−1 − xi
+ f(xi)
x− xi−1
xi − xi−1
s1(x) = 1
(x− 2)
1 − 2
+ 2
(x− 1)
2 − 1
= −x+ 2 + 2x− 2 = x
logo, s1(x) para x ∈ [1, 2].
• s2(x) = 2
(x− 5)
2 − 5
+ 3
(x− 2)
5 − 2
=
2
3
(5 − x) + x− 2 =
1
3
(x+ 4)⇒
s2(x) =
1
3
x+
4
3
se x ∈ [2, 5].
• s3(x) = 3
(x− 7)
5 − 7
+ 2, 5
(x− 5)
7 − 5
=
1
2
(−0, 5x+ 8, 5)⇒
s3(x) = −
1
4
x+
17
4
, se x ∈ [5, 7].
Observação 5.11.1 Podemos ter splines interpolante de diversos graus

grau 1
grau 2
grau 3
grau n
. Nos ateremos
Em
El
ab
or
aç
ão
54 5. Interpolação
ao grau 1, que é a linear.
5.12 Spline Cúbica Interpolante
• A Spline Linear tem a desvantagem de 1◦ derivada ser descontínua nos pontos xi, i ∈ {1, . . . ,n}.
• A Spline Quadrática tem a desvantagem da 2◦ derivada (concavidade de função) ser descontínua
nos pontos xi, i ∈ {1, . . . ,n}.
FIGURA
Para que não haja descontinuidade na 1◦ derivada (crescente ou decrescente) ou na 2◦ derivada
(concavidade
{
∪
∩
)
analisaremos a Spline Cúbica, S3(x).
A Spline Cúbica Interpolante S3(x) é:
S3(x) =

s1(x); x ∈ [x0, x1]
s2(x); x ∈ [x1, x2]
...
sn(x); x ∈ [xn−1, xn]
Observação 5.12.1 Vamos tomar uma função para cada subintervalo e utilizar as informações
adicionais para fixar a continuidade das derivadas.
Observação 5.12.2
• Spline
{
função definida
por partes contínua
}
• Cúbica
{
as funções sk(x)
são de grau 3
}
• Interpolante
{
S3(xi) = f(xi)
∀i = 0, . . . ,n
}
Por simplicidade, denotaremos sk(x) da forma
sk(x) = ak(x− xk)
3 + bk(x− xk)
2 + ck(x− xk) + dk
Com 4 constantes {ak,bk, ck e dk} a serem determinadas para cada k, então no total temos 4 n
variáveis a determinar!
Para obter a expressão da Spline Cúbica Interpolante, vamos impor as condições de continuidade para
Os pontos
1◦ derivada
2◦ derivada
i) s3(xi), i = 0, 1, 2, . . . ,n⇒ n+ 1 condições;
Em
El
ab
or
aç
ão
5.12. Spline Cúbica Interpolante 55
ii) sk(xk) = sk+1(xk), i = 1, 2, 3, . . . ,n− 1⇒ n− 1 condições;
iii) s ′k(xk) = s
′
k+1(xk), i = 1, 2, 3, . . . ,n− 1⇒ n− 1 condições;
iv) s ′′k(xk) = s
′′
k+1(xk), k = 1, 2, 3, . . . ,n− 1⇒ condições.
No total, temos 4n− 2 condições e 4n variáveis, ou seja, temos duas condições em aberto que teremos
que definir para fixar a função (pode ser definida a partir de alguma característica do problema).
Obtendo um expressão para sk(x)
Seja sk(x) = ak(x− xk)3 + bk(x− xk)
2 + ck(x− xk) + dk denotando

hk = xk − xk−1
s
′′
k(xk) = gk
f(xk) = gk
• Pela condição i), que S3(xi) = f(xi), i = 0, 1, 2, . . . ,n.
– Para k = 1, 2, 3, . . . ,n (sem o zero)
sk(xk) = f(xk) = dk, então se dk = f(xk)n condições k = 1, 2, . . . ,n
– Para k = 0, temos que S3(x0) = s1(x0) = f(x0)
s1(x0) = a1 (x0 − x1)
3︸ ︷︷ ︸+b1 (x0 − x1)
2︸ ︷︷ ︸+c1 (x0 − x1)︸ ︷︷ ︸+f(x1) = f(x0)
−h1 −h1 −h1
−ah3
1 + b1h
2
1 − c1h1 + f(x1) = f(x0)1 condição
• Pela ii) Como [pela condição i)] já temos que sk(xk) = f(xk), para que sk+1(xk) = sk(xk) basta
que sk+1(xk) = f(xk) tambémo seja. Logo
sk+1(xk) = ak+1(xk − xk+1)
3 + bk+1(xk − xk+1)
2 + ck(xk − xk+1) + f(xk+1) = f(xk)⇒
−ak+1h
3
k+1 + bk+1h
2
k+1 + ckhk + f(xk+1) = f(xk) com k = 1, 2, 3, . . . ,n− 1. n− 1 condições
• Pela iv) s ′′k(xk) = sk+1(xk), k = 1, 2, . . . ,n− 1.
Como s ′′k(x) = 6ak(x− xk) + 2bk, então
s
′′
k(xk) = 2bk ⇒ bk =
s
′′
k(xk)
2
=
gk
2
s
′′
k(xk−1) = 6ak (xk−1 − xk)︸ ︷︷ ︸+2bk = −6akhk + 2bk
−hk
ak =
2bk − s
′′
k(xk−1)
6hk
=
s
′′
k(xk) − s
′′
(xk−1)
6hk
e como queremos a continuidade da segunda derivada
Em
El
ab
or
aç
ão
56 5. Interpolação
sk(xk) = sk+1(xk) ou s
′′
k−1(xk−1) = s
′′
k(xk−1)
ak =
s
′′
k(xk) − s
′′
k−1
6hk
=
gk − gk−1
6hk
, n− 1 condições
Observação 5.12.3 Para k − 1, introduzimos uma variável g0 = s
′′
0 (x0)) que não existia
anteriormente.
Reoganizando os termos, para obter ck em termos de gk.
Como sk+1(xk) = f(xk) [condição ii)] temos
ak+1(xk − xk+1)
3 + bk(xk − xk+1)
2 + ck(xk − xk+1) + dk = f(xk)⇒
−ak+1h
3
k+1 + bk+1h
2
k+1 − ck+1hk+1 + dk+1 = f(xk)
Que é equivalente a
−akh
3
k + bkh
2
k − ckhk + dk = f(xk−1)
ck =
dk − f(xk − 1) − akh3
k + bkh
2
k
hk
Substituindo
ak =
gk − gk−1
6hk
bk =
gk
2
dk = f(xk = yk
ck =
gk − gk−1
hk
+
−
(gk−gk−1h
3
k
6hk +
gkh
2
k
2
hk
ck =
gk − gk−1
hk
+
(−gk + gk−1 + 3gk)
6
hk ⇒
ck =
gk − gk−1
hk
+
2gk + gk−1
6
hk
Assim podemos calcular todos os coeficientes, de sk(x), sabendo o valor de gk, k = 0, 1, 2, . . . ,n.
• Pela condição iii): s ′k(xk) = s ′k+1(xk), k = 1, 2, 3, . . . ,n − 1. A derivada de sk é s ′k(x) = 3ak(x −
xk)
2 + 2bk(x− xk) + ck. Logo s ′k(xk) = ck
sk+1(xk) = 3ak+1(xk − xk+1)
2 + 2bk+1(xk − xk+1) + ck+1
= 3ak+1h
2
k+1 − 2bk+1hk+1 + ck+1 = ck
Logo,
ck+1 = ck − 3ak+1h
2
k+1 + 2bk+1hk+1
Em
El
ab
or
aç
ão
5.12. Spline Cúbica Interpolante 57
Substituindo as expressões para

ak = gk−gk−1
6hk
ab = gk
2
ck =
[
gk−gk−1
hk
+ 2gk+gk−1hk
6
]
gk+1 − yk
hk+1
+
(2gk+1 + gk)hk+1
6
=
gk − gk−1
hk
+
(2gk + gk−1)hk
6
− 3
(gk+1 − gkhk+1
6hk+1
+ 2
gk+1hk+1
2
6
yk+1 − yk
hk+1
+6
−yk + yk+1
hk
= −2yk+1hk+1−ykhk+1−ykhk+1+2ykhk+yk+1hk−3yk+1hk+1+3ykhk+1+6yk+1hk+1
yk+1hk+1 + 2yk(hk + hk+1) + gk−1hk = 6
yk+1 − yk
hk+1
− 6
yk − yk−1
hk
Um sistema com
{
n+ 1 incógnitas {y0,y1,y2, . . . ,yn}
n− 1 equações k = 1, 2, 3, . . . ,n− 1.
Portanto indeterminado
[A]n−1×(n+1)[x]n+1×1 = [b](n−1)×1
Ax = b
A =

h1 2(h1 + h2) h2 0 0 . . . 0
0 h2 2(h2 + h3) h3 0 . . . 0
0 0 h3 2(h3 + h4) h4 . . . 0
...
...
...
...
...
. . .
...
0 0 0 hn−2 2(hn−2 + hn−1) hn−1 0
0 0 0 0 hn−1 (2h3 + hn−1) hn

x =

y0
y1
y2
...
yn
 b =

y2−y2
h2
(y1−y0)
h1
y3−y2
h3
(y2−y1)
h2
...
yn−yn−1
hn
(yn−1−yn−2)
hn−1

Para obter uma solução unívoca, temos que fixar duas condições:
- Se y0 = 0 e yn = 0, a Spline é chamada de natural (equivalente a supor que nos extremos a aproximação
é linear);
- Se y0 = y1 e yn = yn−1;
- Podemos inserir outros valores para y0 e yn, a depender do problema analisado.
Para o caso em que temos o mesmo espacamento entre os pontos,
{
hk = h, ∀ k = 1, 2, 3, . . . ,n
e y0 = yn = 0, temos a seguinte expressão com k = 1, 2, 3, . . . ,n− 1
com y0 = 0 e yn = 0.
Em
El
ab
or
aç
ão
58 5. Interpolação
hyk−1 + 4hgk + hgk+1 =
6
h
(yk+1 − 2yk + y0)

4hy1 + hy0 = 6
h
(y2 − 2y1 + y0)
hy1 + 4hy2 + hy3 = 6
h
(y3 − 2y2 + y1)
hy2 + 4hy2 + hy3 = 6
h
(y3 − 2y2 + y1)
...
hyn−3 + 4hyn−2 + hyn−1 = 6
h
(yn−1 − 2yn−2 + yn−3)
hyn−2 + 4hyn−1 = 6
h
(yn − 2yn−1 + yn−2)
⇓

4h h 0 0 0 . . . 0 0
h 4h h 0 0 . . . 0 0
0 h 4h h 0 . . . 0 0
0 0 h 4h h . . . 0 0
...
...
...
...
...
. . .
...
...
0 0 0 0 0 h 4h h
0 0 0 0 0 . . . h 4h


y1
y2
y3
y4
...
yn−2
yn−1

=
6
h

y2 − 2y1 + y0
y3 − 2y2 + y1
y4 − 2y3 + y2
y5 − 2y4 + y3
...
yn−1 − 2yn−2 + yn−3
yn − 2yn−1 + yn−2

Portanto, dado um conjunto de pontos igualmente espaçados, e com a condição que y0 = 0 e yn = 0, obtemos
um sistema linear com n− 1 incógnitas e n− 1 equações cuja solução permite obter os coeficientes da Spline
Cúbica Interpolante S3(x).
ak =
yk − yk−1
6h
bk =
yk
2
dk = yk ck =
yk − yk−1
h
+
2yk + yk−1
6
h
Exemplo 5.12.1
Obtenha a Spline Cúbica Interpolante do conjunto de dados
x 0 0, 5 1, 0 1, 5 2, 0
f(x) y0 = 3 y1 = 2 y2 = −1 y3 = −4 y4 = −9Como temos S pontos, temos 4 regiões, então
S3(x) =

s1(x), x ∈ [0; 0, 5]
s2(x), x ∈ [0, 5; 1, 0]
s3(x), x ∈ [1, 0; 1, 5]
s4(x), x ∈ [1, 5; 2, 0]
Como o espaçamento é constante
h = 0, 5 para obter os coeficientes de
sk(x) precisamos resolver o sistema.
hyk−1 + 4hyk + hyk+1 =
6
h
(yk+1 − 2yk + y0)
Em
El
ab
or
aç
ão
5.12. Spline Cúbica Interpolante 59
k = 1
k = 2
k = 3

4hy1 + hy2 =
6
h
(y2 − 2y1 + y0)
hy1 + 4hy2 + hy3 =
6
h
(y3 − 2y2 + y1)
hy2 + 4hy3 =
6
h
(y4 − 2y3 + y2)
⇒

2y1 + 0, 5y2 =
6
0, 5
[(−1) − 2 · 2 + 3] = −24
0, 5y1 + 2y2 + 0, 5y3 =
6
0, 5
[(−4) − 2(−1) + 2] = 0
0, 5y2 + 2y3 =
6
0, 5
[(−9) − 2(−4) − 1] = +24

4y1 + y2 = −48
y1 + 4y2 + y3 = 0
y2 + 4y3 = 48
4y1 + y2 = −48
−4y1 − 16y2 − 4y3 = 0
y2 + 4y3 = 48
⇒
y0 = 0 y1 = −12 y2 = 0 y3 = 12
y4 = 0
Logo ak =
yk − yk−1
6h
; bk =
yk
2
; dk = yk ck =
yk − yk−1
h
+
2hyk + yk−1h
6
• s1(x)
a1 =
y1 − y0
6h
=
−12 − 0
6× 0, 5
= −4;
b1 =
y1
2
= −6;
d1 = y1 = 2;
c1 =
y1 − y0
h
+ 2
y1 + y0
6
h =
2 − 3
0, 5
+ 2
(−12)
6
· 0, 5
= −2 − 2 = −4
s1(x) = −4(x− 0, 5)3 − 6(x− 0, 5)2 − 4(x− 0, 5)2
• s2(x)
a2 =
y2 − y1
6h
=
12
6× 0, 5
= 4;
b2 =
y2
2
= 0;
d2 = −1;
c2 =
y2 − y1
h
+ 2
y2 + y1
6
h =
−1 − 2
0, 5
+ 2
(0 − 12)
6
· 0, 5
= −6 − 1 = −7
s2(x) = −4(x− 1)3 − 7(x− 1) − 1
Em
El
ab
or
aç
ão
60 5. Interpolação
• s3(x)
a3 =
y3 − y2
6h
=
12
6× 0, 5
= 4;
b3 =
y3
2
= 6;
d3 = y3 = −4;
c3 =
y3 − y2
h
+ 2
y3 + y2
6
h =
−4 + 1
0, 5
+ 2
(12)
6
· 0, 5
= −6 + 2 = −4
s3(x) = 4(x− 1, 5)3 + 6(x− 1, 5)2 − 4(x− 1, 5) − 4
• s4(x)
a4 =
y4 − y3
6h
=
−12
6× 0, 5
= −4;
b4 =
y4
2
= 0;
d4 = y4 = −9;
c4 =
y4 − y3
h
+ 2
y4 + y3
6
h =
−9 − (−4)
0, 5
+ 2
(12)
6
· 0, 5
= −10 + 1 = −9
s4(x) = −4(x− 2)3 − 9(x− 2) − 9
Portanto a Spline Cúbica Interpolante para esse conjunto de dados é:
S3(x) =

s1(x) = −4(x− 0, 5)3 − 6(x− 0, 5)2 − 4(x− 0, 5) + 2, se x ∈ [0; 0, 5]
s2(x) = 4(x− 1)3 − 7(x− 1) − 1, se x ∈ [0, 5; 1, 0]
s3(x) = 4(x− 1, 5)3 + 6(x− 1, 5)2 − 4(x− 1, 5) − 4, se x ∈ [1, 0; 1, 5]
s4(x) = −4(x− 2)3 − 9(x− 2) − 9, se x ∈ [1, 5; 2, 0]
Em
El
ab
or
aç
ão
6
Ajuste de Curvas pelo Método dos
Mínimos Quadrados
Uma outra forma de se analisar um conjunto de dados (x1, f(x1), . . . , (xn, f(xn)) é tentar ajustar um
conjuntos de curvas que melhor aproximam estes dados.
Portanto, fixando n funções conhecidas y1(x),y2(x), . . . ,yn(x) contínuas em [a,b] com xi ∈ [a,b] ∀
i ∈ {1, 2, . . . ,n}, podemos obter as constantes α1,α2, . . . ,αn tais que
ϕ(x) = α1y1(x) + α2y2(x) + α3y3(x) + · · ·+ αnyn(x)
aproxime a função f(x) dos dados.
Observação 6.0.1 As funções yi(x) são escolhidas a partir de observações dos conjuntos dos dados
ou utilizando uma fundamentação teórica.
Exemplo 6.0.1
FIGURAS
Portanto, nosso objetivo é encontrar os melhores parâmetros α1, . . . ,αn que ajustam as funções co-
nhecidas y1(x), . . . ,yn(x) aos (x1, f(x1)), (x2, f(x2)), . . . , (αn, f(xn)), isto é, os melhores parâmetros que
aproximam ϕ(x) = α1y1(x) + α2y2(x) + · · ·+ αnyn(x) a função f(x).
Existem diversas formas de se obter esses parâmetros, mas estudaremos o Método dos Mínimos Qua-
drados.
Observação 6.0.2 Também podemos propor o mesmo problema para analisar uma função contínua,
mas não abordaremos tal tema neste curso.
6.1 Método dos Mínimos Quadrados
Sejam m pontos (x1, f(x1)), . . . , (xn, f(xn)) e n funções conhecidas y1(x),y2(x), . . . ,yn(x), com m > n,
queremos encontrar os parâmetros α1,α2, . . . ,αn que melhor aproxima a função ϕ(x) = α1y1(x) + · · · +
αnyn(x) a f(x).
Definindo dk = f(xk) − ϕ(xk), o método dos mínimos quadrados consiste em obter os αi de forma
que a soma dos quadrados dos desvios seja mínima.
61
Em
El
ab
or
aç
ão
62 6. Ajuste de Curvas pelo Método dos Mínimos Quadrados
F(α1,α2, . . . ,αn) =
m∑
k=1
a2
k =
m∑
k=1
(f(xk) −ϕ(xk))
2)
=
m∑
k=1
[f(xk) − α1y1(xk) − α2y2(xk) − · · ·− αnyn(xk)]2
Observação 6.1.1 Se ϕ(x) for o polinômio interpolador de grau m, que passa por todos os pontos,
então ϕ(xk) = f(xk) e consequentemente o desvio dk = ϕ(xk) −ϕ(xk) = 0⇒ F(α1,α2, . . . ,αn) = 0.
Logo a interpolação polinomial é um caso especial para o ajuste de funções
Para encontrar o mínimo de F(α1, . . . ,αn) vamos derivar em relação aos parâmetros e igualar a zero.
∂F
∂αj
= 2 ·
m∑
k=1
[f(xk) − α1y1(xk) − · · ·− αnyn(xk)][−yj(xk)]
e como
∂F
∂αj
= 0 temos que
m∑
k=1
[f(xk) − α1y1(xk) − · · ·− αnyn(xk)] · yj(xk) = 0
Fazendo isso para todos os parâmetros j = 1, 2, . . . ,n
m∑
k=1
[f(xk) − α1y1(xk) − · · ·− αnyn(xk)] · y1(xk) = 0
m∑
k=1
[f(xk) − α1y1(xk) − · · ·− αnyn(xk)] · y2(xk) = 0
...
m∑
k=1
[f(xk) − α1y1(xk) − · · ·− αnyn(xk)] · yn(xk) = 0
Efetuando o produto
m∑
k=1
f(xk)y1(xk) − α1
m∑
k=1
y1(xk)y1(xk) − · · ·− αn
m∑
k=1
yn(xk)y1(xk) = 0
m∑
k=1
f(xk)y2(xk) − α1
m∑
k=1
y1(xk)y2(xk) − · · ·− αn
m∑
k=1
yn(xk)y1(xk) = 0
...
...
...
m∑
k=1
f(xk)yn(xk) − α1
m∑
k=1
y1(xk)yn(xk) − · · ·− αn
m∑
k=1
yn(xk)y1(xk) = 0
que irá formar um sistema linear, com incógnitas αi, i = {1,n}
Denotando aij =
∑m
k=1 yi(xk)yj(xk) e bi =
∑m
k=1 f(xk)yi(xk).
a11α1 + a12α2 + a13α3 + · · ·+ a1nαn = b1
a21α1 + a22α2 + a23α3 + · · ·+ a2nαn = b2
...
an1α1 + an2α2 + an3α3 + · · ·+ annαn = bn
Em
El
ab
or
aç
ão
6.1. Método dos Mínimos Quadrados 63
Esse sistema de equações lineares recebe o nome de equações normais e é possível mostrar que a sua
solução é realmente um mínimo e existe se yi(x), . . . ,yn(x) são funções linearmente independentes.
Se analisarmos (yi(x1),yi(x2), . . . ,yi(xn)) como um vetor os elementos aij podem ser vistos como
produtos escalares.
yi = (yi(x1),yi(x2), . . . ,yn(xm))
então
aij = yiyj =< yi;yj >=
m∑
k=1
yi(xk)yj(xk)
Observação 6.1.2 aij = aji, então a matriz [A] é simétrica, ou seja, A = AT .
Definição 6.1.1 Os vetores y1, . . . ,yn com yi = (yi(x1), . . . ,yi(xn)) não ditos ortogonais, e conse-
quentemente, as funções y1(x), . . . ,yn(x) não ortogonais, se
< yi;yj >=
m∑
k=1
yi(xk)y
xk
j =
{
= 0, sei 6= j
6= 0, sei = j
Para polinômios ortogonais, o sistema anterior se torna diagonal, com solução direta
αi =
< f;yi >
< yi;yi >
Exemplo 6.1.1
Seja
x −1 −0, 5 0 0, 5 1
f(x) 2, 05 0, 4 0 0, 55 2
determine o melhor ajuste pelo método dos mínimos quadrados para o conjunto de dados acima,
para as funções
a)
{
y1(x) = 1
y2(x) = x
ou seja, encontrar a melhor reta da forma
ϕ(x) = α1y1(x) + α2y2(x)
ϕ(x) = α1 + α2x
yi = (yi(x1),yi(x2), . . . ,yi(x5)
y1 = (1, 1, 1, 1, 1); y1(x) = 1
y2 = (−1;−0, 5; 0; 0, 5; 1); y2(x) = x
Em
El
ab
or
aç
ão
64 6. Ajuste de Curvas pelo Método dos Mínimos Quadrados

a11 = y1 · y1 = (1, 1, 1, 1, 1)(1, 1, 1, 1, 1) = 5
a12 = y1 · y2 = (1, 1, 1, 1, 1)(−1;−0, 5; 0; 0, 5; 1) = 0
a22 = y2 · y2 = 1 + 0, 52 + 0 + 0, 52 + 1 = 2, 5
f = (2, 05; 0, 4; 0; 0, 51; 2)
então
b1 = f · y1 = 2, 05 + 0, 4 + 0, 51 + 2 = 4, 96
b2 = f · y2 = −2, 05 − 0, 4 · 0, 5 + 0 + 0, 51 · 0, 5 + 2 · 1 = 0, 005
logo
{
a11α1 + a12α2 = b1
a21α1 + a22α2 = b2
{
5α1 = 4, 96
2, 5α2 = 0, 005
⇒
α1 =
4, 96
5
= 0, 992
α2 = 0, 002
portanto,
ϕ(x) = 0, 992 + 0, 002x
b) Se
{
y1(x) = x
y2(x) = x
2 − 1
, então
y1 = (−1;−0, 5; 0; 0, 5; 1)
y2 = (0;−0, 75;−1;−0, 75; 0)
a11 = y1 · y1 = 1 + 0, 52 + 0 + 0, 52 + 1 = 2, 5
a12 = y1 · y2 = 0, 5 · 0, 75 − 0, 5 · 0, 75 = 0
a22 = y2 · y2 = 2, 125
b1 = f · y1 = −2, 05 − 0, 2 + 0, 255 + 2 = 0, 005
b2 = f · y2 = −0, 4 · 0, 75 − 0, 75 · 0, 51 = 0, 6825
portanto,
{
a11α1 + a12α2 = b1
a21α1 + a22α2 = b2
⇒
2, 5α1 = 0, 005⇒ α1 = 0, 002
2, 125α2 = 0, 6825⇒ α2 = 0, 321
Com isso,
Em
El
ab
or
aç
ão
6.1. Método dos Mínimos Quadrados 65
ϕ(x) = α1y1(x) + α2y2(x) = α1x+ α2x
2 − α2 ⇒
ϕ(x) = −0, 321 + 0, 002x+ 0, 321x2
Em
El
ab
or
aç
ão
7
Soluções Numéricas para Equações
Diferenciais Ordinárias
Equações Diferenciais são equações que envolvem derivadas das funções
dy
dx
= y+ x2 ;
d2y
dx2 +
dy
dx
+ y = 2
y2∂y(x, t)
∂x
+
∂y(x, t)
∂t
− 2= 0
Quando a função só possui uma variável, a equação diferencial é chamada de Equação Diferencial
Ordinário (EDO). Se a ED envolve mais de uma variável dependente, temos Equações Diferenciais Parciais
(EDP). Neste curso analisaremos somente as EDO. A ordem de uma equação diferencial é a ordem de
maior derivada.
Exemplo 7.0.1
y ′′(x) + xy ′ = 5 (EDO de 2◦ ordem)
y ′′′′′(x) − y ′′(x) + y = x (EDO de 5◦ ordem)
Uma EDO é dita linear, se a função y(x) e as suas derivadas aparecem linearmente na equação.
Estudaremos as Equações Diferenciais Ordinárias Lineares. Dado uma EDO, dydx = y, queremos en-
contrar a função y(x) que satisfaz a equação.
Neste caso se y = ex, temos que dy(x)dx = ex = y(x). Logo, y(x) = ex é uma solução. Mas também
podemos perceber que y1 = 2ex ⇒ dy
dx = 2ex = y também é solução.
De modo geral, qualquer y(x) = aex, com a ∈ R∗ é uma solução de equação dy
dx = y (temos uma
família de soluções). Para determinar univocamente a solução, precisamos de um valor inicial, ou seja,
sobre alguma informação adicional sobre a função, como

y(0) = 1
y(2) = 5
y(1) = 3
. ⇒

y(x) = ex
y(x) = 5ex−2
y(1) = 3ex−1
Se tivermos uma equação de ordemm seria necessáriom valores iniciais para determinar univocamente
a solução.
Nos ateremos a soluções para EDO de 1◦ ordem, já que sempre podemos reduzir uma EDO de ordem
m para um sistema de EDO de 1◦ ordem.
Exemplo 7.0.2
y(
′′′)(x) + x2y(
′′) +
(cos x)
x3 y
′
(x) = tg(x)
Fazendo
{
y ′(x) = z(x)
y ′′(x) = w(x)
66
Em
El
ab
or
aç
ão
7.1. Método de Euler 67
ẏ
′
(x) =

y(x)
z(x)
w(x)


y ′(x) = z(x)
z ′(x) = w(x)
w ′(x) = −x2w(x) −
cos x
x3 z(x) + y(x) + tg(x)
então resolver a equação diferencial ordinária de 3◦ ordem y(
′′′)(x) + x2y(
′′) +
cos x
x3 y
′
− y(x) = tg x é
equivalente a resolver o sistema
y ′(x)
z ′(x)
w ′(x)
 =

z(x) = z(x)
w(x) = w(x)
−x2w(x) −
cos x
x3 z(x) + y(x) + tg(x)

fazendo
{
z(x) = y ′(x)
w(x) = z ′(x) = y ′′(x)
.
Portanto, focaremos na resolução da EDO de 1◦ ordem, a partir do seguinte problema de valor
inicial (PVI). {
y ′(x) = f(x,y)
y(x0) = y0
7.1 Método de Euler
Seja
{
y ′ = f(x,y)
y(x0) = y0
o PVI, então o Método de Euler consiste em obter os valores de y(xi) = yi a
partir da equação da reta que passa pelo ponto inicial (x0,y0) e tem coeficiente angular f(x,y), ou seja,
r0(x) = y(x0) + (x− x0)y
′(x0).
FIGURA
denotando
{
h = xk+1 − xk e considerando que
y(x1) = y1 ≈ r(x1) = y0 + hf(x0,y0)
. Fazendo isso para k = 0, 1, 2, 3, . . . temos
yk+1 = yk + hf(xk,yk)
Observação 7.1.1 Para obter o próximo ponto, só precisamos saber do anterior. O erro pode ser
estipulado por
E(xn) =
y
′′
(xn)
2
h2 (Erro no cálculo do n-ésimo ponto xn)
Exemplo 7.1.1
Seja
{
y
′
= y
y(0) = 1
, utilizando um passo de h = 0, 2, calcule y(1) pelo Método de Euler.
Seja

y
′
= y︸︷︷︸
f(x,y)
y(0) = 1
, temos que
{
x0 = 0
y0 = 1
, tomando um passo de h = 0, 2:
Em
El
ab
or
aç
ão
68 7. Soluções Numéricas para Equações Diferenciais Ordinárias
x1 = x0 + h = 0, 2
x2 = x0 + 2h = 0, 4
x3 = x0 + 3h = 0, 6
x4 = x0 + 4h = 0, 8
x5 = x0 + 5h = 1, 0
Logo yk+1 = yk + hf(xk,yk),
{
h = 0, 2 e y0 = y(0) = 1
f(xk,yk) = yk
, logo
yk+1 = yk + 0, 2yk = 1, 2yk ⇒
y1 = 1, 2 · 1 = 1, 2
y2 = 1, 2y1 = 1, 2 · 1, 2 = 1, 44
y3 = 1, 2y2 = 1, 2 · 1, 44 = 1, 728
y4 = 1, 2y3 = 1, 2 · 1, 728 = 2, 0736
y5 = 1, 2y4 = 2, 48832
Logo,
y(1) ' y5 = 2, 49
FIGURA
Observação 7.1.2 A solução exta é y(x) = ex onde y(1) = e = 2, 71.
O Método de Euler faz parte de um conjunto de métodos que utilizam a Série de Taylor. Expandindo
a solução y(x) do PVI,
{
y
′
(x) = f(x,y)
y(x0) = y0
em torno do ponto xn temos
y(x) = y(xn) + (x− xn)y
′(xn) +
(x− xn)
2
2!
y
′′
(xn) +
(x− xn)
3
3!
y
′′′
(xn) + · · ·
então, se tomarmos a aproximação até a 1◦ ordem da série de Taylor
y(x) ' y(xn) + (x− xn)y
′
(xn)
e como y ′(xn) = f(xn,yn) e x = xn+1, com xn+1 − xn = h
y(xn+1) = yn+1 ' yn + hf(xn,yn)
Observação 7.1.3 Por isso o erro é E(x) =
h2
2
f
′′
(x)!
Podemos considerar aproximações de ordem superiores, para melhorar resultados, mas ao aumentarmos
a ordem, teríamos que calcular
•
y
′′
(x) =
d
dx
f(x,y(x)) =
∂f
∂x
+
∂f
∂y
dy
dx
= fx + fyy
′
onde estamos denotando ∂f
∂x = fx e ∂f∂y = fy. Como y ′ = f(x,y)
Em
El
ab
or
aç
ão
7.2. Método de Runge-Kutta 69
•
y
′′
(x) = fx + fyf ; y
′′′
(x) = ; y
′′′′
(x) = ; yk(x)
Apesar de termos uma expressão para as derivadas, dependendo da complexidade da função f(x,y)
encontrar as diversas derivadas de f(x,y) pode ser extremamente trabalhoso.
Para contornar o trabalho de se calcular as diversas derivadas de f(x,y), surge o Método de Runge-
Kutta.
7.2 Método de Runge-Kutta
Seja o PVI
{
y ′ = f(x,y)
y(x0) = y0
, um método é dito ser de Runge-Kutta de ordem p se
i) São de passo simples, ou seja, para calcular yn+1 só é necessário saber o ponto anterior yn;
ii) Não exigem o cálculo de qualquer uma das derivadas de f(x,y);
iii) A expressão para a solução coincide com a expressão da série de Taylor até a ordem p.
Observação 7.2.1 O Método de Euler é o Método de Runge-Kutta de ordem 1 (RK1).
Método de Runge-Kutta de Ordem 2 (RK2)
O Método de Runge-Kutta de ordem 2, também conhecido como Método de Euler Aperfeiçoado, pode ser
obtido substituindo no Método de Euler o coeficiente angular da reta, que é f(xn,yn) por
1
2
[
f(xn,yn) + f(xn + h,yn + hy
′
n)
]
o que resulta em:
yn+1 = yn +
h
2
[f(xn,yn) + f(xn + h,yn + hf(xn,yn))]
Exemplo 7.2.1
Seja o PVI
{
y ′ = y
y(0) = 1
, calcule y(1) pelo Método de RK2 com passo h = 0, 2.
Como h = 0, 2, temos que
{
x1 = 0, 2
x2 = 0, 4
;
{
x1 = 0, 6
x2 = 0, 8
; {x5 = 1, 0, logo
yn+1 = yn +
h
2
[f(xn,yn) + f(xn + h,yn + hf(xn,yn))]
Como f(x,y) = y, obtemos
f(xn + h,yn + hyn) = yn + hyn
yn+1 = yn +
h
2
[yn + (yn + hyn)] = yn + hyn +
h2
2
yn
yn+1 = yn
(
1 + h+
h2
2
)
Em
El
ab
or
aç
ão
70 7. Soluções Numéricas para Equações Diferenciais Ordinárias
yn+1 = yn
(
1 + 0, 2 +
0, 04
2
)
= yn(1, 22)
y1 = y0 × 1, 22 = 1× 1, 22 = 1, 22
y2 = y1 × 1, 22 = 1, 22× 1, 22 = 1, 4884
y3 = y2 × 1, 22 = 1, 4884× 1, 22 = 1, 815848
y4 = y3 × 1, 22 = 1, 815484× 1, 22 = 2, 21535456
y5 = y4 × 1, 22 = 2, 702708
Observação 7.2.2 O valor exato é y(x) = ex ⇒ y(1) = e = 2, 71828.
Podemos ter várias expressões do Método Runge-Kutta para uma mesma ordem, além do Método RK
para diversas ordens,
Este método é extessivamente empregado na resolução de diversas equações diferenciais ordinárias
(EDO).
Observação 7.2.3 O método de Euler modificado temos que
yn+1 = yn +
h
2
(k1 + k2) ;
{
k1 = f(xn,yn)
k2 = f(xn+h,yn = hk1)
O método é dito de 2◦ ordem porque coincide com a expanssão em série de Taylor até 2◦ ordem. De
fato
y ′ = f(x,y)⇒ y(x+ h) = y(x) + h · y ′(x) + h2
2!
y ′′(x) +O(3)
= y(x) + hf(x,y) +
h2
2
d
dx
[f(x,y)]
y(x = h) = y(x) + hf(x,y) +
h2
2
[
∂f
∂x
+
∂f
∂y
y ′
]
+O(3)
= y(x) + hf(x,y) +
h2
2
(fx + fyf) +O(3) (Série de Taylor)
y(x+ h) = y(x) +
h
2
f(x,y) +
h
2
[f(x+ h ; y+ hf(x,y)] RK2
Para comparar as expressões, vamos expandir o 3◦ termo (K2) de RK2 em série de Taylor em torno de
(x y) h = 0.
f[x+ h ; y+ hf(x,y)] = f(x,y) + h
[
∂
∂x
f
]
|h=0 + h
[
∂f
∂y
]
|h=0f(x,y) +O(2)
substituindo em RK2
y(x+ h) = y(x) +
h
2
f+
h
2
[f+ hfx + hyfyf+O(2)]
= y(x) + hf+
h2
2
[fx + fyf] +O(3)
Em
El
ab
or
aç
ão
7.3. Método das Diferenças Finitas 71
Coincide até a 2◦ ordem com a expansão em série de Taylor.
Podemos ter outras expressões para RK2, desde que coincidam com a série de Taylor até 2◦ ordem.
De forma similar, podemos definir RK 3◦ ordem.
• Método de Runge-kutta de 3◦ ordem (ou Método de Nystrom)
yn+1 = yn +
h
4
[
k1 +
3
2
k2 +
3
2
k3
]
;

k1 = f(xn;yn)
k2 = f(xn + 2
3h ; yn + 2
3hk1)
k3 = f(xn + 2
3h ; yn + 2
3hk2)
• Método de Runge-kutta de 4◦ ordem
yn+1 = yn +
h
8
[k1 + 3k23k3 + k4]

k1 = f(xn,yn)
k2 = f(xn + h
3 ; yn + h
3 k1
k3 = f(xn + 2
3h; yn − h
3 k1 + hk2)
k4 = f(xn + h; yn + hk1 − hk2 + hk3)
7.3 Métododas Diferenças Finitas
A ideia geral do Método de Diferenças Finitas é a discretização do domínio e a substituição das derivadas
por aproximação envolvendo somente os valores numéricos da função.
Seja x0 ∈ R, definimos a maçha de passo h como os pontos
xi = x0 ± ih; i = 1, 2, 3, . . ., N
Expandindo em série de Taylor y(x+ h), temos que
y(x± h) = y(x)± hy ′(x4) +
h2
2
y ′′(x) +O(3)⇒
y1(x) = ±y(x+ h) − y(x)
h
∓ h
2
y ′′(x) +O(2)︸ ︷︷ ︸ ' y(x+ h)∓ y(x)h
O(1)
• Se considerarmos x+ h, temos:
y ′(x) ' y(x+ h) − y(x)
h
[Diferenças Progressivas]; erro O(1)
• Se considerarmos x− h, temos:
y(x) ' −
y(x− h) − y(x)
h
[Diferenças Regressivas]; erro O(1)
Em
El
ab
or
aç
ão
72 7. Soluções Numéricas para Equações Diferenciais Ordinárias
• Somando as duas diferenças
y ′(x) ' y(x+ h) − y(x− h)
h
[Diferença Central]; erroO(2)
De modo análogo podemos obter uma expressão para a 2◦ derivada de y(x).
y(x+ h) = y(x) + hy ′(x) +
h2
2!
y ′′(x) +
h3
3!
y ′′′(x) +
h4
4
yiv(x) +O(5)
y(x− h) = y(x) − hy ′(x) +
h2
2!
y ′′(x) −
h3
3!
y ′′′(x) +
h4
4
yiv(x) +O(5)
Somando as duas expressões, temos
y ′′(x) =
y(x+ h) − ∂y(x) + y(x− h)
h2 −
h2
12
yiv(x) +O(3)⇒
y ′′(x) ' y(x+ h) − ∂y(x) + y(x− h)
h2 ; com erro O(2)
Desta forma, podemos obter expressões para derivadas de ordem superior. Seguindo o mesmo raciocí-
nio, podemos estender as relações obtidas para gerar aproximações para as derivadas parciais.
Seja u(x, t) um função de duas variáveis
∂tu(x, t) ' u(x, t+ k) −U(x, t)
k
[Diferença Progressiva]
∂xu(x, t) ' u(x+ h, t) − u(x, t)
h
• Diferença regressiva
∂tu(x, t) =
u(x, t) − u(x, t− k)
k
∂xu(x, t) =
u(x, t) − u(x− h, t)
h
• Diferença Central
Ux(x, t) ' U(x+ h, t) −U(x− h, t)
2h
Ut(x, t) ' U(x, t+ k) −U(x, t− k)
2k
Uxx(x, t) ' U(x+ h, t) − ∂U(x, t) +U(x− h, t)
h2
Utt(x, t) ' U(x, t+ h) − ∂U(x, t) +U(x, t− k)
k2
Uxt(x, t) ' U(x+ h, t+ k) −U(x− h, t+ k) −U(x+ h, t− k) +U(x− h, t− k)
4hk
Assim, equações diferenciais parciais de 2◦ ordem podem ser revolvidas discretizendo o domínio, a
partir do Método das Diferenças Finitas.
Exemplo 7.3.1
a(x,y)Uxx + b(x,y)Uxy + C(x,y)Uyy + F(x,y,Ux,Uy,U) = 0
Em
El
ab
or
aç
ão
7.3. Método das Diferenças Finitas 73
Seja 4 = b2 − 4ac, então
→ ∗4 = 0⇒ Equação Parabólica; Ex: Ut − α(x, t)Uxx = r(x, t) (Equação do calor, difusão)
→ ∗4 < 0⇒ Equação Elíptica; Ex: Uxx +U(y,y) = 0 (Equação de Laplace)
→4 > 0⇒ Equação Hiperbólica; Ex: Uxx −U(t, t) = 0 (Equação de Onda)
Observação 7.3.1 (*) Apresentam soluções altamente regulares.
Resolver uma EDP pelo Método das Diferenças Finitas consiste em discretizar a região do espaço.
Em
El
ab
or
aç
ão
Em
El
ab
or
aç
ão
Índice Remissivo
Sistemas, 1
75
	Sistemas Numéricos
	Introdução
	Sistemas Binário e Decimal
	Representação de Inteiros
	Conversão Binário para Decimal
	Conversão Decimal para Binário
	ExercÃcios
	Representação de Frações
	ExercÃcios
	Aritmética de Ponto Flutuante
	Instabilidades Numéricas
	ExercÃcios
	Zeros Reais de Funções Reais
	Fase I: Isolamento das Raízes
	Fase II: Refinamento
	Método da Bissecção
	Estimativa do Número de Iterações
	Método da Posição Falsa
	Método do Ponto Fixo
	Método de Newton-Raphson
	Método da Secante
	Métodos Estudados
	Sistemas de Equações Lineares
	Introdução
	Sistemas Lineares
	Resolução de Sistemas Lineares
	Métodos Diretos
	Introdução à Resolução de Sistemas Não Lineares
	Método de Newton
	Método de Newton Modificado
	Interpolação
	Interpolação Polinomial
	Formas para se obter o polinômio interpolador Pn(x)
	Resolução do Sistema Linear
	Forma de Lagrange
	Forma de Newton
	Forma de Newton para pn(x)
	Estudo do Erro na Interpolação
	Interpolação Inversa
	Fenômeno de Runge
	Spline
	Spline Linear Interpolante
	Spline Cúbica Interpolante
	Ajuste de Curvas pelo Método dos Mínimos Quadrados
	Método dos Mínimos Quadrados
	Soluções Numéricas para Equações Diferenciais Ordinárias
	Método de Euler
	Método de Runge-Kutta
	Método das Diferenças Finitas
	Índice Remissivo