apostila_mned

Prévia do material em texto

Métodos Numéricos
para
Equações Diferenciais
Helio Pedro Amaral Souto
Nova Friburgo, 3 de Abril de 2014
Graduação em Engenharia
Conteúdo
Lista de Figuras v
Lista de Tabelas vii
1 Equações Diferenciais Ordinárias 1
1.1 Métodos “One–Step” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1 Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1.1 Análise do Erro do Método de Euler . . . . . . . . . . . . 3
1.1.1.2 Estabilidade do Método de Euler . . . . . . . . . . . . . . 5
1.1.2 Métodos de Mais Alta Ordem . . . . . . . . . . . . . . . . . . . . . 7
1.1.3 Método de Heun . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.4 Método de Euler Modificado . . . . . . . . . . . . . . . . . . . . . . 9
1.1.5 Métodos de Runge–Kutta . . . . . . . . . . . . . . . . . . . . . . . 10
1.1.5.1 Métodos de Runge–Kutta de Segunda Ordem . . . . . . . 11
1.1.5.2 Métodos de Runge–Kutta de Terceira Ordem . . . . . . . 13
1.1.5.3 Métodos de Runge–Kutta de Quarta Ordem . . . . . . . . 13
1.1.5.4 Sistema de Equações . . . . . . . . . . . . . . . . . . . . . 14
1.2 Formulação de Butcher . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.1 Condições de Ordem . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.2 Alguns Métodos Explícitos . . . . . . . . . . . . . . . . . . . . . . . 19
1.3 Métodos “Multistep” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.3.1 Métodos do Tipo Preditor–Corretor . . . . . . . . . . . . . . . . . . 21
1.3.2 Métodos de Mais Alta Ordem . . . . . . . . . . . . . . . . . . . . . 24
1.3.3 Método de Milne . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
1.3.4 Método de Adams de Quarta Ordem . . . . . . . . . . . . . . . . . 26
1.3.5 Método de Hamming . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.4 Convergência, Consistência e Estabilidade . . . . . . . . . . . . . . . . . . 28
1.4.1 Métodos One-Step . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.4.1.1 Convergência . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.4.1.2 Consistência . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.4.1.3 Estabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . 29
1.4.2 Métodos Multistep . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
1.4.2.1 Convergência . . . . . . . . . . . . . . . . . . . . . . . . . 31
1.4.2.2 Consistência . . . . . . . . . . . . . . . . . . . . . . . . . . 31
i
ii Conteúdo
1.4.2.3 Estabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2 Equações Diferenciais Parciais 35
2.1 Algumas Equações Diferenciais de Interesse . . . . . . . . . . . . . . . . . 35
2.2 Classificação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.2.1 Classificação pelas Características . . . . . . . . . . . . . . . . . . . 38
2.2.2 Equações a N Variáveis Independentes . . . . . . . . . . . . . . . . 41
2.2.3 Transformação de Coordenadas . . . . . . . . . . . . . . . . . . . . 41
2.2.4 Sistema de Equações . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.2.5 Classificação pela Análise de Fourier . . . . . . . . . . . . . . . . . 45
2.3 Problema Bem-Posto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.3.1 Condições de Contorno e Iniciais . . . . . . . . . . . . . . . . . . . 47
2.4 Equações Diferenciais Hiperbólicas . . . . . . . . . . . . . . . . . . . . . . 48
2.4.1 Interpretação em Bases Físicas . . . . . . . . . . . . . . . . . . . . . 49
2.4.2 Condições de Contorno e Inicial Apropriadas . . . . . . . . . . . . . 52
2.5 Equações Diferenciais Parabólicas . . . . . . . . . . . . . . . . . . . . . . . 53
2.5.1 Interpretação em Bases Físicas . . . . . . . . . . . . . . . . . . . . . 55
2.5.2 Condições de Contorno e Inicial Apropriadas . . . . . . . . . . . . . 55
2.6 Equações Diferenciais Elípticas . . . . . . . . . . . . . . . . . . . . . . . . 55
2.6.1 Interpretação em Bases Físicas . . . . . . . . . . . . . . . . . . . . . 56
2.6.2 Condições de Contorno Apropriadas . . . . . . . . . . . . . . . . . 57
3 Equações de Balanço 65
3.1 Equação da Continuidade . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.2 Equações do Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.3 Equação de Energia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4 Discretização 73
4.1 Discretização Espacial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.2 Discretização no Tempo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.3 Expansões em Séries de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.4 Técnica Geral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.5 Acurácia do Processo de Discretização . . . . . . . . . . . . . . . . . . . . 79
4.6 Representação do Tipo Onda . . . . . . . . . . . . . . . . . . . . . . . . . 82
5 Convergência, Consistência e Estabilidade 87
5.1 Convergência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.1.1 Teorema de Lax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.2 Consistência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.2.1 O Esquema FTCS . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.2.2 Esquema Totalmente Implícito . . . . . . . . . . . . . . . . . . . . . 91
5.3 Estabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.3.1 Método da Matriz . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.3.2 Método da Matriz: Condições de Contorno . . . . . . . . . . . . . . 95
5.3.3 Método de Von Neumann . . . . . . . . . . . . . . . . . . . . . . . 96
Conteúdo iii
5.3.4 Método de Von Neumann: Esquema Geral . . . . . . . . . . . . . . 97
5.4 Acurácia da Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
5.4.1 Extrapolação de Richardson . . . . . . . . . . . . . . . . . . . . . . 100
6 Equação de Transporte Linear 103
6.1 Métodos Explícitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
6.2 Métodos Implícitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.3 Equação de Difusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.3.1 Equação de Difusão Unidimensional . . . . . . . . . . . . . . . . . 106
6.3.2 Equação de Difusão em Regime Permanente . . . . . . . . . . . . . 113
6.3.3 Equação de Difusão Bidimensional . . . . . . . . . . . . . . . . . . 114
6.4 Métodos Multidimensionais de Fracionamento . . . . . . . . . . . . . . . . 117
6.5 Equação de Advecção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
6.6 Dissipação e Dispersão Numéricas . . . . . . . . . . . . . . . . . . . . . . . 129
6.6.1 Análise de Fourier . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
6.6.2 Equação Modificada . . . . . . . . . . . . . . . . . . . . . . . . . . 133
6.7 Equação de Transporte Unidimensional . . . . . . . . . . . . . . . . . . . . 135
6.8 Equação de Transporte Bidimensional . . . . . . . . . . . . . . . . . . . . . 147
A O Método de Separação de Variáveis 151
B Algoritmo de Thomas 157
B.1 O Algoritmo de Thomas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
B.2 O Método TDMA Linha a Linha . . . . . . . . . . . . . . . . . . . . . . . 160
B.3 Relaxação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
Bibliografia 163
Índice 165
iv Conteúdo
Lista de Figuras
1.1 Método de Euler. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Método de Heun. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Exemplo de uma árvore. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.1 Curva características no plano t-x. . . . . . . . . . . . . . . . . . . . . . . 39
2.2 Domínio computacional R. . . . . . .. . . . . . . . . . . . . . . . . . . . 48
2.3 Características para a equação da onda. . . . . . . . . . . . . . . . . . . . 49
2.4 Propagação de uma onda senoidal. . . . . . . . . . . . . . . . . . . . . . . 51
2.5 Condições de contorno para o caso hiperbólico. . . . . . . . . . . . . . . . 52
2.6 Condições auxiliares para o caso transiente. . . . . . . . . . . . . . . . . . 53
2.7 Domínio computacional para uma EDP parabólica. . . . . . . . . . . . . . 54
2.8 Decaimento exponencial de uma senóide. . . . . . . . . . . . . . . . . . . 54
2.9 Domínio computacional para uma EDP elíptica. . . . . . . . . . . . . . . 56
4.1 Malha discreta uniforme. . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.1 Exemplo da representação matricial. . . . . . . . . . . . . . . . . . . . . . 118
6.2 Ângulo de fase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
B.1 Método TDMA linha a linha. . . . . . . . . . . . . . . . . . . . . . . . . . 160
v
vi Lista de Figuras
Lista de Tabelas
1.1 Exemplos de leis fundamentais . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Representação da árvore e do operador diferencial elementar . . . . . . . . 17
1.3 Funções da árvore na expansão da série de Taylor . . . . . . . . . . . . . . 18
4.1 Comparação da acurácia dos diferentes esquemas. . . . . . . . . . . . . . . 81
4.2 Relação de amplitudes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
A.1 Separação da equação de Helmohtz. . . . . . . . . . . . . . . . . . . . . . . 154
vii
viii Lista de Tabelas
Capítulo 1
Equações Diferenciais Ordinárias
Vários problemas práticos de engenharia são descritos em termos de taxas de variação
das propriedades físicas do sistema que estamos interessados em estudar. Estas variações
são, normalmente, descritas empregando as leis fundamentais da física, mecânica, eletri-
cidade, termodinâmica, etc (vide a Tabela 1.1). Quando combinadas com as equações
de balanço de massa, momentum e energia elas resultam nas equações diferenciais que
governam diferentes fenômenos físicos. A partir da posterior integração destas equações,
nós obtemos as funções matemáticas que descrevem o estado do sistema, no tempo e no
espaço, em termos da sua massa, energia ou variação da quantidade de movimento.
Tabela 1.1: Exemplos de leis fundamentais que são escritas em termos de taxas de varia-
ção.
Lei Expressão Matemática Variáveis e Parâmetros
Segunda lei de Newton
dv
dt
=
F
m
velocidade (v), força (F )
e massa (m)
Lei de Fourier ~q = −k dT
dx
condutividade térmica
(k) e temperatura (T )
Lei de Fick ~j = −D dc
dx
coeficiente de difusão (D)
e concentração (c)
Como exemplo de uma equação diferencial ordinária, derivada de uma lei fundamental
(lei de Newton), temos a equação que governa a queda de um paraquedista
dv
dt
= g − c
m
v
1
2 Capítulo 1. Equações Diferenciais Ordinárias
onde v é a velocidade do paraquedista, g é a aceleração da gravidade, m é a massa e
c o coeficiente de arrasto. Nesta equação a variável que esta sendo diferenciada, v, é
chamada de variável dependente. A variável em relação a qual v está sendo diferenciada,
t, é denominada de variável independente.
A equação diferencial ordinária, acima, é um exemplo de uma equação de primeira
ordem, porque a sua derivada de mais alta ordem é uma derivada primeira. Em contrapar-
tida, a equação que descreve a posição x de um sistema massa–mola com amortecimento
é um exemplo de uma equação diferencial ordinária de segunda ordem
m
d2x
dt2
+ c
dx
dt
+ k x = 0
onde m é a massa, c o coeficiente de amortecimento e k a constante elástica da mola.
Vale ressaltar que equações diferenciais ordinárias de alta ordem podem ser reduzidas
a um sistema de equações diferenciais ordinárias de primeira ordem. No exemplo acima,
podemos introduzir uma nova variável y definida por
y =
dx
dt
que pode ser derivada com relação ao tempo, de modo que obtemos o sistema
y =
dx
dt
dy
dt
= −c y + k x
m
que é equivalente a equação original de segunda ordem. Uma vez que uma equação
diferencial ordinária de ordem n pode, de modo similar, ser reduzido a um sistema de
equações diferenciais de primeira ordem, nós focaremos nossa atenção na resolução das
equações de primeira ordem.
Na prática, muitas destas equações diferenciais ordinárias, que representam um grande
número de problemas em engenharia, não podem ser resolvidas empregando os métodos
analíticos do cálculo diferencial. Assim, os métodos numéricos de resolução destas equa-
ções diferenciais adquirem uma grande importância em vários campos da engenharia.
A fim de que possamos especificar completamente a solução de uma equação diferencial
ordinária, devemos fornecer condições auxiliares. Em se tratando de equações diferenciais
de primeira ordem, uma condição auxiliar chamada de condição inicial é necessária para
garantirmos a obtenção de uma única solução.
Quando estivermos trabalhando com equações diferenciais de ordem n, n condições
serão requeridas para que possamos obter uma única solução. Caso todas as condições
sejam especificadas para um mesmo valor de uma variável independente (por exemplo para
t = 0), então o problema será chamado de um problema de valor inical. Por outro lado,
um problema de valor de contorno será definido quando estas condições forem fornecidas
para diferentes valores da variável independente.
1.1. Métodos “One–Step” 3
1.1 Métodos “One–Step”
No que se segue, estaremos interessados em resolvermos equações diferenciais ordiná-
rias do tipo
dy
dx
= f(x, y) (1.1)
para uma condição inicial dada por y(x0) = y0, que é conhecido como um Problema de
Valor Inicial (PVI).
Uma primeira possibilidade, para que possamos resolver esta equação, pode ser ob-
tida a partir da discretização do termo da derivada de primeira ordem, empregando um
esquema de diferenças avançadas:
dy
dx
∣∣∣∣
i
=
yi+1 − yi
∆x
ou ainda
yi+1 = yi +
dy
dx
∣∣∣∣
i
∆x
De acordo com esta equação, a partir da estimativa da derivada (que fornece a inclinação
da tangente à curva no ponto xi) nós determinamos, por extrapolação, um novo valor de
y no ponto xi+1 a partir de um valor conhecido de yi. Esta fórmula pode ser aplicada
passo a passo a fim de determinarmos a evolução futura da solução.
Todos os métodos do tipo “one–step” podem ser expressos nesta mesma forma geral,
sendo que a única diferença entre eles será a maneira de avaliar a derivada primeira.
1.1.1 Método de Euler
Uma primeira possibilidade, de estimativa da derivada, é fornecida pela própria ex-
pressão da equação diferencial ordinária e consiste em empregarmos a própria função a fim
de avaliarmos a inclinação da tangente no ponto xi, conforme esquematizado na Fig. 1.1.
Portanto, neste método a estimativa é dada por
yi+1 = yi + f (xi, yi) ∆x (1.2)
onde esta fórmula é conhecida como o método de Euler ou de Euler–Cauchy. Nela, o
novo valor de y é determinado a partir de uma extrapolação linear empregando o valor
da derivada em xi (Fig. 1.1).
1.1.1.1 Análise do Erro do Método de Euler
A resolução numérica da equação diferencial ordinária, de modo semelhante ao caso
de uma equação diferencial parcial, introduz dois tipos de erros:
• Um erro de truncamento ou discretização devido à natureza da técnica empregada
na aproximação de y;
4 Capítulo 1. Equações Diferenciais Ordinárias
y
xx xi i+1
yi
y i+1
verdadeiro
erro
Figura 1.1: Método de Euler.
• Um erro de arredondamento devido ao fato do número limitado de digitos que são
utilizados pelos computadores.
O erro de truncamento é composto de duas partes. Uma primeira, chamada de erro
local de truncamento, que resulta da aplicação local do método numérico para um único
incremento. A segunda é a parcela devida à propogação do erro de truncamento proveni-ente das aproximações decorrentes dos incrementos anteriores. A soma destes dois erros
fornece o erro global de truncamento.
Expandindo em série de Taylor a função y, que possui as suas derivadas contínuas, em
torno do ponto xi,
yi+1 = yi + y
′
i ∆x+
y′′i
2
∆x2 + . . .+
y
(n)
i
n!
∆xn +Rn
onde y′ = dy/dx e Rn representa o resto da série e é definido por
Rn =
y(n+1)(ξ)
(n+ 1)!
∆xn+1
onde ξ está contido no intervalo de xi a xi+1. Esta expansão ainda pode ser escrita numa
forma alternativa, se considerarmos que estamos tratando da equação diferencial ordinária
y′ = f(x, y)
yi+1 = yi + f(xi, yi) ∆x+
f ′(xi, yi)
2
∆x2 + . . .+
f (n−1)(xi, yi)
n!
∆xn +O (∆xn+1) (1.3)
Comparando esta expansão com o método de Euler, Eq. (1.2), vemos que este último
pode ser escrito como
yi+1 = yi + f(xi, yi) ∆x+ Et
1.1. Métodos “One–Step” 5
onde Et é o verdadeiro erro local devido ao truncamento, sendo dado por
Et =
f ′(xi, yi)
2
∆x2 + . . .+O (∆xn+1)
Para valores suficientemente pequenos do incremento ∆x, os diferentes termos do erro
decrescem à medida que a sua ordem aumenta. Portanto, usualmente representamos o
erro como
Ea =
f ′(xi, yi)
2
∆x2
ou
Ea = O
(
∆x2
)
onde Ea representa o erro local de truncamento aproximado.
Conforme pudemos ver, a expansão em série de Taylor nos fornece um meio para quan-
tificarmos o erro no método de Euler. Entretanto, existem algumas limitações associadas
ao seu uso:
• A série de Taylor nos fornece somente uma estimativa do erro de truncamento local.
Ela não nos fornece o erro propagado, ou seja, o erro global de truncamento;
• Na prática, podemos trabalhar com funções que são mais complexas do que simples
polinômios. Como consequência, o cálculo das derivadas que aparecem na série de
Taylor podem não ser facilmente obtidas.
Apesar dessas limitações impedirem uma análise exata do erro, na maior parte dos
problemas de interesse a série de Taylor ainda pode nos fornecer informações valiosas
sobre o comportamento do método de Euler. Do desenvolvimento, em série de Taylor,
pudemos verificar que o erro local é proporcional ao quadrado do incremento de espaço e
à derivada primeira da equação diferencial. É possível mostrar, também, que o erro global
de truncamento é de O(∆x) [1], ou seja, ele é proporcional ao incremento de espaço. Estas
observações nos levam a algumas conclusões importantes:
• O erro pode ser reduzido mediante um decréscimo do incremento;
• O método fornecerá resultados livres de erro caso o solução da equação diferencial
ordinária seja linear, uma vez que as derivadas de ordem superiores a um serão
nulas.
Pelas razões mencionadas, acima, o método de Euler é denominado um método de
primeira ordem.
1.1.1.2 Estabilidade do Método de Euler
O conceito matemático de estabilidade do método numérico está associado ao cresci-
mento ilimitado da solução numérica ou, conforme veremos mais adiante, do erro intro-
duzido nos dados iniciais ou no processo de discretização. Neste caso, a solução numérica
6 Capítulo 1. Equações Diferenciais Ordinárias
não irá tender para a solução da equação diferencial quando o incremento tender a zero e
o esquema numérico será dito ser instável [1, 2].
Este conceito pode ser entendido tomando como exemplo o seguinte problema de valor
inicial:
dy
dx
= −α y
apresentando como condição inicial y(0) = y0. Do Cálculo Diferencial sabemos que esta
equação diferencial ordinária possui a seguinte solução:
y = y0e
−αx
Assim, essa solução tem como valor inicial y0 e aproxima–se assintoticamente de zero à
medida que o valor de x aumenta.
Empregando o método de Euler na resolução numérica deste problema, temos que:
yi+1 = yi +
dy
dx
∣∣∣
i
∆x = yi + f(xi, yi)∆x = yi − αyi∆x
ou ainda,
yi+1 = (1− α∆x) yi
Deste resultado vemos claramente que a estabilidade do método numérico depende do
incremento ∆x, isto é, o valor de |1− α∆x| deve ser menor do que 1. Esta restrição
implica no fato de que o incremento deve verificar a condição ∆x < 2/α a fim de que
a solução numérica acompanhe a solução analítica. Neste caso, dizemos que o método é
condicionalmente estável.
Uma maneira de contornarmos esta restrição seria a de empregarmos uma abordagem
implícita para o método de Euler. Esta representação é chamada de implícita devido ao
fato de que a variável que buscamos determinar, yi+1, aparece em ambos os lados do sinal
de igualdade
yi+1 = yi +
dy
dx
∣∣∣
i+1
∆x = yi + f(xi+1, yi+1)∆x = yi − αyi+1∆x
ou ainda,
yi+1 =
1
1 + α∆x
yi
sendo que agora a condição de estabilidade,
∣∣∣∣ 11 + α∆x
∣∣∣∣ < 1, é sempre verificada para
α > 0. Assim, dizemos que o método de Euler implícito é incondicionalmente estável.
Este conceito de estabilidade pode ser estendido a outros métodos numéricos destina-
dos à resolução de um problema de valor inicial do tipo:
dy
dx
= f(x, y), x > a
para uma condição inicial y(a) = y0. Os leitores interessados devem se reportar às refe-
rências [2, 3] para maiores detalhes. Este assunto será retomado ao final deste capítulo.
1.1. Métodos “One–Step” 7
1.1.2 Métodos de Mais Alta Ordem
Uma maneira fácil de reduzirmos o erro do método de Euler seria a de incluirmos
termos de mais alta ordem da série de Taylor (1.3) na solução. Por exemplo, retendo o
termo de segunda ordem
yi+1 = yi + f(xi, yi)∆x+
f ′(xi, yi)
2
∆x2 + Ea
onde o erro de truncamento aproximado seria dado por
Ea =
f ′′(xi, yi)
6
∆x3
Embora este método seja facilmente implementado no caso de funções polinomiais, a
inclusão de termos de mais alta ordem pode não ser trivial no caso de estarmos consi-
derando equações diferenciais ordinárias mais complicadas. Em particular, este é o caso
em se tratando de equações diferenciais ordinárias que são funções de ambas as variáveis
dependentes e independentes, sendo necessário o emprego da regra da cadeia na determi-
nação destas derivadas. Por exemplo, a derivada primeira de f(x, y) será
f ′(x, y) =
∂f
∂x
+
∂f
∂y
dy
dx
enquanto que a derivada segunda seria dada por
f ′′(x, y) =
∂f ′
∂x
+
∂f ′
∂y
dy
dx
com o grau de complexidade aumentando para as derivadas de mais alta ordem.
Assim sendo, alternativas foram pensadas para que possamos aumentar a performance
dos métodos one–step, sem que tenhamos que avaliar derivadas superiores a primeira
ordem. Conforme veremos a seguir, estes métodos são comparáveis aos métodos que
empregam termos de mais alta ordem na série de Taylor.
1.1.3 Método de Heun
O método de Euler apresenta uma fonte de erro fundamental que é devida ao fato de
que a derivada, avaliada no início do intervalo, é suposta prevalecer sobre todo o intervalo.
Em seguida, veremos que algumas modificações simples podem ser introduzidas a fim de
superarmos esta dificuldade.
O método de Heun utiliza a determinação de duas derivadas no intervalo, como medida
para melhorar a estimativa da inclinação da curva no intervalo considerado. As derivadas
são avaliadas uma no ponto inicial e a outra no ponto final. A estimativa da inclinação
é obtida, então, a partir da média dos valores das duas derivadas, conforme mostrado
esquematicamente na Fig. 1.2.
No método de Euler, a inclinação é avaliada no ponto inicial
y′i = f(xi, yi)
8 Capítulo 1. Equações Diferenciais Ordinárias
y
xx xi i+1
yi
y i+1
0
Figura 1.2: Método de Heun.
e é empregada na extrapolação linear que determina yi+1
y0i+1 = yi + f(xi, yi)∆x
Enquanto que para o método de Euler esta seria a estimativa final, no método de Heun
trata-se de uma estimativa intermediária e ela é conhecida como a equação preditora. Ela
fornece uma estimativa de yi+1 com a qual determinamos uma estimativa da inclinação à
curva no ponto final do intervalo:
y′i+1 = f(xi+1, y
0
i+1)
Assim, as duas expressões para as inclinações podem sercombinadas a fim de obtermos
um valor médio para a inclinação no intervalo considerado:
y¯′ =
y′i + y
′
i+1
2
=
f(xi, yi) + f(xi+1, y
0
i+1)
2
Este valor é, então, empregado na extrapolação linear de yi a yi+1 utilizando o método de
Euler
yi+1 = yi +
f(xi, yi) + f(xi+1, y
0
i+1)
2
∆x
onde esta é a equação corretora.
Portanto, o método de Heun é um método do tipo preditor–corretor. Neste método,
devemos observar que yi+1 aparece em ambos os lados do sinal de igualdade e, portanto,
esta equação pode ser empregada de modo iterativo a fim de corrigirmos este valor. Isto
quer dizer que a estimativa anterior y0i+1 pode ser usada repetidamente a fim de prover
uma melhor estimativa de yi+1.
1.1. Métodos “One–Step” 9
Na introdução do método de Heun, consideramos que a derivada era uma função
tanto da variável dependente quanto da independente. Entretanto, para casos como os
dos polinômios, onde a equação diferencial ordinária é função unicamente da variável
independente, não necessitamos empregar a equação preditora e basta utilizarmos uma
única vez a cada iteração a equação corretora
yi+1 = yi +
f(xi) + f(xi+1)
2
∆x
onde podemos notar uma similaridade entre o lado direito do sinal de igualdade e a regra
do trapézio para a integração. De fato, partamos da equação diferencial ordinária
dy
dx
= f(x)
Esta equação pode ser resolvida por integração na forma:∫ yi+1
yi
dy =
∫ xi+1
xi
f(x) dx
que resulta em
yi+1 = yi +
∫ xi+1
xi
f(x) dx
Agora, aplicando a regra do trapézio no intuito de avaliarmos a integral acima
yi+1 = yi +
f(xi) + f(xi+1)
2
∆x− f
′′(ξ)
12
∆x3
onde ξ encontra–se no intervalo entre xi e xi+1. Assim, trata–se de um método de segunda
ordem, uma vez que a derivada segunda da equação diferencial ordinária é zero quando a
solução for quadrática. Além do mais, os erros local e global são de O (∆x3) e de O (∆x2)
respectivamente [1].
1.1.4 Método de Euler Modificado
Uma outra modificação simples do método de Euler resulta no método de Euler modi-
ficado. Esta técnica usa o método de Euler na predição de um valor de y no ponto médio
do intervalo
yi+1/2 = yi + f(xi, yi)
∆x
2
Em seguida, este valor é empregado na estimativa do valor da inclinação à curva no ponto
médio:
y′i+1/2 = f(xi+1/2, yi+1/2)
onde assumimos que este valor é uma aproximação da inclinação média para todo o
intervalo. Este valor da inclinação é então usado na extrapolação linear entre xi e xi+1
empregando o método de Euler
yi+1 = yi + f(xi+1/2, yi+1/2)∆x
10 Capítulo 1. Equações Diferenciais Ordinárias
Devemos notar que devido ao fato de yi+1 não aparecer em ambos os lados desta equação,
esta equação corretora não pode ser empregada iterativamente a fim de melhorarmos a
solução.
O método de Euler modificado é superior ao método de Euler visto que ele avalia a
inclinação à curva no ponto médio do intervalo de predição. Da teoria sobre as técnicas
de discretização, sabemos que a técnica de diferenças finitas conhecida como diferenças
centradas aproxima melhor a derivada do que as técnicas de diferenças avançadas ou
recuadas. No mesmo sentido, uma aproximação centrada como a utilizada aqui possui um
erro de truncamento de O (∆x2) em comparação com a aproximação avançada do método
de Euler que possui um erro de O (∆x). Conseqüentemente, este método apresenta um
erro local e global que são de O (∆x3) e de O (∆x2) respectivamente [1].
1.1.5 Métodos de Runge–Kutta
Os métodos do tipo Runge–Kutta podem alcançar a mesma acurácia dos métodos de
mais alta ordem que empregam o desenvolvimento em série de Taylor sem, no entanto,
requerer a determinação das derivadas de mais alta ordem. Muitas variações destes méto-
dos existem, mas todos podem ser escritos na seguinte forma generalizada em se tratando
dos métodos explícitos:
yi+1 = yi + φ(xi, yi,∆x)∆x (1.4)
onde φ(xi, yi,∆x) é chamada de função incremento a qual pode ser interpretada como a
inclinação representativa de todo o intervalo. A função incremento pode ser escrita numa
forma geral como
φ = b1k1 + b2k2 + · · ·+ bnkn
onde os diferentes coeficientes bi, i = 1, 2, . . . são constantes e as funções ki, i = 1, 2, . . .
são dadas por
k1 = f(xi, yi)
k2 = f(xi + c2∆x, yi + a21k1∆x)
k3 = f(xi + c3∆x, yi + a31k1∆x+ a32k2∆x)
...
kn = f(xi + cn∆x, yi + an1k1∆x+ an2k2∆x+ · · ·+ ann−1kn−1∆x)
Devemos notar que nestas expressões as funções ki são relações de recorrência. Isto quer
dizer que k1 aparece na equação para k2, que aparece na equação para k3 e assim por
diante. Tal fato torna os métodos de Runge–Kutta eficientes ferramentas do cálculo
numérico.
Vários tipos de métodos de Runge–Kutta podem ser imaginados a partir do emprego
de diferentes números de termos da função incremento. Devemos notar que o método de
Runge–Kutta de primeira ordem com n = 1 corresponde, de fato, ao método de Euler.
Uma vez escolhido o número n de termos, os valores de a, p e q são avaliados igualando a
equação generalizada do método de Runge–Kutta aos termos de uma expansão em série
de Taylor. Deste modo, ao menos para as versões de baixa ordem, o número de termos n
normalmente representa o ordem do método.
1.1. Métodos “One–Step” 11
1.1.5.1 Métodos de Runge–Kutta de Segunda Ordem
A forma generalizada para um método de segunda ordem é dada por
yi+1 = yi + (b1k1 + b2k2)∆x
onde
k1 = f(xi, yi)
k2 = f(xi + c2∆x, yi + a21k1∆x)
A fim de determinarmos os valores de b1, b2, c2 e a21 nós iremos expandir em série de
Taylor yi+1 e empregaremos f(xi, yi) como sendo a derivada dy/dx
yi+1 = yi + f(xi, yi)∆x+ f
′(xi, yi)
∆x2
2
+O (∆x3)
onde f ′(xi, yi) deve ser determinada a partir do emprego de:
df(xi, yi) =
∂f
∂x
∣∣∣
i
dx+
∂f
∂y
∣∣∣
i
dy
ou ainda
f ′(xi, yi) =
df
dx
=
∂f
∂x
∣∣∣
i
+
∂f
∂y
∣∣∣
i
dy
dx
Substituindo este resultado na expansão em série de Taylor
yi+1 = yi + f(xi, yi)∆x+
(
∂f
∂x
+
∂f
∂y
dy
dx
)
i
∆x2
2
+O (∆x3) (1.5)
Em seguida, nós vamos expandir a equação geral (Eq. (1.4)) para o método de segunda
ordem, lembrando que
g(x+ r, y + s) = g(x, y) + r
∂g
∂x
+ s
∂g
∂y
+
1
2!
(
∂2g
∂x2
r2 + 2
∂2g
∂x∂y
r s+
∂2g
∂y2
s2
)
+ · · ·
Portanto,
k2 = f(xi, yi) + c2 ∆x
∂f
∂x
+ a21 k1 ∆x
∂f
∂y
+O (∆x2)
Substituindo, agora, k1 e k2 na expressão geral (1.4) nós obtemos
yi+1 = yi+ b1∆x f(xi, yi) + b2∆x f(xi, yi) + b2c2∆x
2 ∂f
∂x
+ b2a21∆x
2 f(xi, yi)
∂f
∂y
+O (∆x3)
ou ainda, agrupando os termos de mesma ordem em ∆x
yi+1 = yi + (b1 + b2)f(xi, yi)∆x+
[
b2c2
∂f
∂x
+ b2a21 f(xi, yi)
∂f
∂y
]
∆x2 +O (∆x3) (1.6)
12 Capítulo 1. Equações Diferenciais Ordinárias
Finalmente, comparando termo a termo as expansões (1.5) e (1.6), vemos que as
seguintes condições devem ser verificadas para que tenhamos a equivalência destas duas
equações
b1 + b2 = 1
b2c2 =
1
2
b2a21 =
1
2
Como existem quatro incógnitas e somente três equações, não há uma única solução
que satisfaça este sistema de equações. Necessitamos, portanto, especificar o valor de uma
destas incógnitas a fim de determinarmos as outras três. Em consequência, teremos uma
família de métodos de segunda ordem. Admitamos que nós especifiquemos o valor de b2,
das equações acima obtemos
b1 = 1− b2
c2 = a21 =
1
2b2
Devido ao fato de podermos escolher uma infinidade de valores para a2, existirão uma
infinidade de métodos de Runge–Kutta de segunda ordem, que fornecerão o mesmo resul-
tado para a solução da equação diferencial ordinária em se tratando de uma função f(x, y)
constante, linear ou quadrática. Entretanto, obteremos soluções diferentes quando esta
função for mais complicada. Em seguida, apresentaremos três dos métodos de segunda
ordem mais comumente utilizados:
Método de Heun (b2 = 1/2).Para b2 = 1/2 nós podemos resolver o sistema de
equações e obtemos b1 = 1/2 e c2 = a21 = 1. Substituindo esses resultados na equação
geral obtemos:
yi+1 = yi +
1
2
(k1 + k2) ∆x
onde
k1 = f(xi, yi)
k2 = f(xi + ∆x, yi + k1∆x)
Notemos que k1 representa a inclinação à curva no ponto inicial do intervalo enquanto
que k2 é a inclinação no final do intervalo. Portanto, este método de Runge–Kutta de
segunda ordem corresponde ao método de Heun com uma única iteração corretora.
Método do Polígono (b2 = 1). Caso assumamos que b2 = 1, então b1 = 0, c2 = a21 =
1/2 e a equação geral torna–se
yi+1 = yi + k2∆x
onde
k1 = f(xi, yi)
k2 = f (xi + ∆x/2, yi + k1∆x/2)
o que corresponde ao método do polígono melhorado.
1.1. Métodos “One–Step” 13
Método de Ralston (b2 = 2/3). Ralston (1962) e Ralston e Rabinowitz (1978) deter-
minaram que para b2 = 2/3 obtemos o valor mínimo para o erro de truncamento para os
algoritmos de Runge–Kutta de segunda ordem. Para esta versão temos que b1 = 1/3 e
c2 = a21 = 3/4
yi+1 = yi +
1
3
(k1 + 2k2) ∆x
onde
k1 = f(xi, yi)
k2 = f (xi + 3∆x/4, yi + 3k1∆x/4)
1.1.5.2 Métodos de Runge–Kutta de Terceira Ordem
Podemos empregar o mesmo desenvolvimento aplicado na obtenção do método de
segunda ordem ao caso de n = 3. O resultado desta derivação produzirá um sistema
de seis equações a oito incógnitas. Portanto, devemos fornecer valores para duas destas
incógnitas a fim de que possamos determinar os parâmetros restantes. Apresentaremos
agora uma versão comumente empregada
yi+1 = yi +
1
6
(k1 + 4k2 + k3) ∆x
onde
k1 = f(xi, yi)
k2 = f (xi + ∆x/2, yi + k1∆x/2)
k3 = f (xi + ∆x, yi −∆x k1 + 2∆x k2)
Para uma função f dependente unicamente de x, o método de terceira ordem reduz-se à
regra de Simpson (1/3) quando empregamos um polinômio interpolador de Lagrange de
segunda ordem
yi+1 = yi +
∫ b
a
f(x) dx
onde
I =
∫ b
a
f(x) dx ≈ ∆x
6
[f(x0) + 4f(x1) + f(x2)]
onde x0 e x2 representam os pontos inicial e final do intervalo ∆x.
Os métodos de Runge-Kutta de terceira ordem apresentam um erro local de O (∆x4)
e um erro global de O (∆x3) e os seus resultados são exatos quando a solução for cúbica.
1.1.5.3 Métodos de Runge–Kutta de Quarta Ordem
Os métodos de Runge–Kutta mais populares são os de quarta ordem. Assim como
no caso dos métodos de segunda ordem, existe uma infinidade de versões destes métodos.
Apresentaremos, aqui, o método que é conhecido como o método clássico de Runge–Kutta
de quarta ordem:
yi+1 = yi +
1
6
(k1 + 2k2 + 2k3 + k4) ∆x
14 Capítulo 1. Equações Diferenciais Ordinárias
onde
k1 = f(xi, yi)
k2 = f (xi + ∆x/2, yi + k1∆x/2)
k3 = f (xi + ∆x/2, yi + k2∆x/2)
k4 = f (xi + ∆x, yi + k3∆x)
De modo análogo ao caso do método de terceira ordem, em se tratando de uma função de
x somente (k2 = k3), o método de Runge–Kutta de quarta ordem é equivalente à regra
de Simpson (1/3) para a integração conforme mostrado anteriormente.
1.1.5.4 Sistema de Equações
Muitos problemas de engenharia e das ciências exatas necessitam da solução de um
sistema de equações diferenciais simultâneas, ao invés da resolução de uma única equação.
Tais sistemas podem ser representados como
dy1
dx
= f1 (x, y1, y2, . . . , yn)
dy2
dx
= f2 (x, y1, y2, . . . , yn)
...
...
dyn
dx
= fn (x, y1, y2, . . . , yn)
A fim de que possamos obter a solução de tal sistema, necessitamos fornecer n condições
iniciais para o valor inicial de x.
Todos os métodos estudados neste capítulo podem ser aplicados na resolução de sis-
temas de equações. Alguns sistemas de aplicações em engenharia podem ser composto
de centenas de equações. Em cada caso, o procedimento adotado na resolução destes
sistemas consiste na aplicação dos métodos “one–step” a todas as equações a cada passo
antes de prosseguirmos para o passo posterior.
1.2 Formulação de Butcher
A família dos métodos de Runge-Kutta pode ser representada seguindo o princípio
descrito por Butcher, nas quais os coeficientes são estrategicamente colocados em uma
tabela para facilitar o cálculo de suas propriedades. A forma geral para determinarmos a
solução aproximada é dada por
yn+1 = yn + ∆x
m∑
i=1
biki
1.2. Formulação de Butcher 15
para 1 ≤ i ≤ s, sendo que nesta equação ki = f(xn + ci∆x, Yi), i = 1, 2, . . . , s para a
forma non-autonomous, ou ki = f(Yi) para os sistemas autonomous, e Yi é dado por
Yi = yn + ∆x
s∑
j=1
aijkj
sendo s o número de estágios do método.
Os coeficientes característicos (aij, bi e ci), que vão determinar a acurácia e a estabili-
dade desses métodos, são determinados convenientemente através da tabela de Butcher [4]
c A
bT
c1 a11 a12 · · · a1s
c2 a21 a22 · · · a2s
...
...
... . . .
...
cs as1 as2 . . . ass
b1 b2 · · · bs
Os coeficientes do vetor c estão relacionados com os elementos da matriz A segundo
a seguinte forma
ci =
s∑
j=1
aij i = 1, 2, ..., s.
e o vetor b deve ser de tal que
s∑
i=1
bi = 1
As componentes do vetor c são os valores intermediários dos estágios, ou seja,
ci ∈ [xn, xn+1] para i = 1, 2, ..., s
As formulações para as quais aij = 0 se j ≥ i são chamadas de formulações explícitas
ou ERK (Explicit Runge-Kutta)
c A
bT
c1 0 0 · · · 0
c2 a21 0 · · · 0
...
...
... . . .
...
cs as1 as2 . . . 0
b1 b2 · · · bs
Neste curso iremos considerar somente as formulações explícitas dos métodos de Runge-
Kutta.
1.2.1 Condições de Ordem
A ordem p dos métodos multi-passos lineares são construídas usando aproximações
para y(xn+1) em termos de y e ∆x y′ avaliados nos respectivos pontos xi das soluções.
16 Capítulo 1. Equações Diferenciais Ordinárias
Dessa forma, a aproximação deve satisfazer exatamente a um polinômio para a expansão
na série de Taylor para um grau menor que p, dado que os erro são estimados em termos
de O(∆xp+1). Para o método de Runge-Kutta torna-se um pouco mais complicado visto
que a aproximação descrita pela Eq. (1.2) é baseada na integral
y(xn+1) = y(xn) + ∆x
∫
y′(xn + ξ∆x)dξ ≈ y(xn) + ∆x
s∑
i=1
biy
′(xn + ci∆x)
Para lidar com este problema, devemos proceder da seguinte forma: encontrar a ex-
pansão da série de Taylor da solução exata, depois encontrar a expansão da série de Taylor
para o cálculo da solução aproximada utilizando o método de Runge-Kutta e, por fim,
comparar as duas expansões termo a termo de modo a que elas sejam equivalentes até o
termo de O(∆xp+1).
Primeiro devemos encontrar a expansão para y que satisfaça o PVI (1.1)
dy
dx
= f(y)
com condição inicial dada por y(xn) = yn e, então, uma formulação para as demais
derivadas dessa função
y′(x) = f(y(x)),
y′′(x) = f ′(y(x))y′(x)
= f ′(y(x))f(y(x)),
y′′′(x) = f ′′(y(x))f(y(x))y′(x) + f ′(y(x))f ′(y(x))y′(x)
= f ′′(y(x))(f(y(x)), f(y(x))) + f ′(y(x))f ′(y(x))f(y(x)).
considerando somente os termos até a derivada terceira. Esses termos são, então, reescritos
na forma
y′(x) = f
y′′(x) = f′f
y′′′(x) = f′′(f, f) + f′f′f
sendo f = f(y(x)), f′ = f ′(y(x)), f′′ = f ′′(y(x)). A partir do resultado destas derivações,
podemos encontrar um padrão capaz de ser utilizado em uma estrutura do tipo árvore.
Nesta estrutura relacionamos f(y(x)) a uma folha da árvore, f ′(y(x)) a um vértice com
uma única extremidade exterior e f ′′(y(x)) a um vértice com duas extremidades exter-
nas. Em se tratando de f ′ e f ′′ as extremidades externas estão conectadas a sub-árvores,
representando estes operadores linear e bilinear respectivamente.
A Tabela 1.2 mostra a representação da estrutura em árvores e as suas respectivas
representações para o operador elementar diferencial F(ψ), considerando cada árvore
ψ ∈ Ψ, sendo Ψ o conjunto de todas as árvores com raízes.
1.2. Formulação de Butcher 17
Tabela 1.2: Representação da árvore e do operador diferencial elementar.
Diagrama F(ψ) Termos
• f f f(y(x))•
•
f′
f
f′f f ′(y(x))f(y(x))
•
• •
��@@
ff
f′′ f
′′(f, f) f ′′(y(x))(f(y(x)), f(y(x)))
•
•
•
f′
f′
f
f′f′f f ′(y(x))f ′(y(x))f(y(x))
A solução exata do problema de valor inicial y(x0 + ∆x) pode ser expandida em série
de Taylor em torno de x0 na forma
y(x0 + ∆x) = y0 +
∑
ψ∈Ψ
α(ψ)∆xr(ψ)
r(ψ)!
F(ψ)(y0)
que pode ser reescrita como
y(x0 + ∆x) = y0 +
∑
ψ∈Ψ
∆xr(ψ)
σ(ψ)γ(ψ)
F(ψ)(y0) (1.7)
onde r(ψ) representa a ordem da árvore ψ (número de vértices ou nós da árvore), σ(ψ)
a simetria de ψ (a ordem do grupo de automorfismo), γ(ψ) a densidade da árvore ψ,
α(ψ) =
r(ψ)!
σ(ψ)γ(ψ)
o número de combinações ordenadas em ψ, β(ψ) =
r(ψ)!
σ(ψ)
o número de
combinações não ordenadas em ψ, F(ψ)(y0) a diferencial elementar. O valor de γ(ψ) para
cada árvore pode ser obtido da seguinte forma, associamos um fator para cada vértice da
árvore, as folhas possuem fator unitário e todos os demais a soma dos fatores associados
aos vizinhos para onde crescem mais um, o produto de todos os fatores é o valor de γ(ψ).
O polinômio Φ(ψ), peso elementar, é construído a partir dos coeficientes de cada método,
baseado em sua respectiva árvore, onde b representa a raiz, c as folhas e A os vértices
internos. Em seguida é fornecido um exemplo de uma árvore na Fig. 1.3.
Por outro lado, também é possível expandir em série de Taylor a solução aproximada
y1 = y0 +
∑
ψ∈Ψ
β(ψ)∆xr(ψ)
r(ψ)!
Φ(ψ)F(ψ)(y0)
ou ainda
y1 = y0 +
∑
ψ∈Ψ
∆xr(ψ)
σ(ψ)
Φ(ψ)F(ψ)(y0) (1.8)
18 Capítulo 1. Equações Diferenciais Ordinárias
Figura 1.3: Exemplo de uma árvore.
•
• ••
• •
�
�
��
@
@
@@
@
@
@@
�
�
��
i
j −→
•
• ••
• •
�
�
��
@
@
@@
@
@
@@
�
�
��
bi
aij
cici
cj cj
•
• ••
• •
�
�
��
@
@
@@
@
@
@@
�
�
��
6
3
11
1 1
Φ(ψ) =
s∑
i,j=1
bic
2
i aijc
2
j γ(ψ) = 18
Como a obtenção das famílias de métodos de Runge-Kutta baseia-se no fato de que
os termos das séries de Taylor (Eqs. (1.7) e (1.8)) devam coincidir até o termo de ordem
∆xp, portanto, é necessário garantir que as condições de ordem [5]
Φ(ψ) =
1
γ(ψ)
, ∀ r(ψ) ≤ p
sejam verificadas.
A Tabela 1.3 apresenta as funções F(ψ), Φ(ψ) e os valores de r(ψ), σ(ψ), γ(ψ), α(ψ)
e β(ψ) para uma expansão considerando somente os termos até quarta ordem.
Tabela 1.3: Funções da árvore na expansão da série de Taylor.
ψ F(ψ) Φ(ψ) r(ψ) σ(ψ) γ(ψ) α(ψ) β(ψ)
• f ∑ bi 1 1 1 1 1
•
•
f′f
∑
bici 2 1 2 1 2
•
• •
��@@ f′′(f, f)
∑
bic
2
i 3 2 3 1 3
•
•
•
f′f′f
∑
biaijcj 3 1 6 1 6
•
•• •
@@�� f(3)(f, f, f)
∑
bic
3
i 4 6 4 1 4
•
• •
•
��@@ f′′(f, f′f)
∑
biciaijcj 4 1 8 3 24
•
•
• •
��@@
f′f′′(f, f)
∑
biaijc
2
j 4 2 12 1 12
•
•
•
•
f′f′f′f
∑
biaijajkck 4 1 24 1 24
1.2. Formulação de Butcher 19
1.2.2 Alguns Métodos Explícitos
Na abordagem para a determinação dos coeficientes característicos, o usual é primeiro
determinarmos c para, em seguida, calcularmos b. Feito isto, podemos encontrar os
coeficientes da matriz A mediante a resolução de um sistema de equações lineares.
Para um método explícito de segunda ordem, da Tabela 1.3 vemos que as condições
que devem ser verificadas são dadas por
• b1 + b2 = 1
•
•
b2c2 =
1
2
chegamos à solução arbitrária em termos do parâmetro não nulo c2, uma vez que temos
duas equações e três incógnitas,
0
c2 c2
1− 1
2c2
1
2c2
Escolhendo c2 =
1
2
ou c2 = 1, obtemos dois casos especiais que são respectivamente as
formulações apresentadas nas regras do ponto médio e do trapézio desenvolvidas por [6]
0
1
2
1
2
0 1
0
1 1
1
2
1
2
Para um método de terceira ordem, as condições necessárias agora são dadas por
• b1 + b2 + b3 = 1
•
•
b2c2 + b3c3 =
1
2
•
• •
��@@ b2c22 + b3c
2
3 =
1
3
•
•
•
b3a32c2 =
1
6
sendo que alguns casos especiais podem ser recuperados em função dos resultados das
tabelas
0
1
2
1
2
1 −1 2
1
6
2
3
1
6
0
2
3
2
3
2
3
0 2
3
1
4
3
8
3
8
0
2
3
2
3
0 −1 1
0 3
4
1
4
20 Capítulo 1. Equações Diferenciais Ordinárias
bem como a formulação conhecida de Heun
0
1
3
1
3
2
3
0 2
3
1
4
0 3
4
No último exemplo, consideramos um método de quarta ordem
• b1 + b2 + b3 + b4 = 1 (1.9)
•
•
b2c2 + b3c3 + b4c4 =
1
2
(1.10)
•
• •
��@@ b2c22 + b3c
2
3 + b4c
2
4 =
1
3
(1.11)
•
•
•
b3a32c3 + b4a42c2 + b4a43c3 =
1
6
(1.12)
•
•• •
@@�� b2c32 + b3c
3
3 + b4c
3
4 =
1
4
(1.13)
•
• •
•
��@@ b3c3a32c3 + b4c4a42c2 + b4c4a43c3 =
1
8
(1.14)
•
•
• •
��@@
b3a32c
2
2 + b4a42c
2
2 + b4a43c
2
3 =
1
12
(1.15)
•
•
•
•
b4a43a32c2 =
1
24
(1.16)
Essas equações devem ser resolvidas em termos de c como parâmetros e no caso de b a
partir das Equações (1.9), (1.10), (1.11) e (1.13). Posteriormente, devemos determinar os
coeficientes da matriz A a partir das Eqs. (1.12), (1.14) e (1.15). E, finalmente, utilizar a
Eq. (1.16) para obtermos a condição de consistência em c. Esta condição de consistência
é encontrada para c4 = 1.
Teorema 1 (Butcher [5]) Se um método de Runge-Kutta explícito com s = 4 possui
ordem 4, então
s∑
i=j+1
biaij = bj(1− cj), j = 1, 2, 3, 4
e, em particular, c4 = 1.
Kutta [7] classificou todas as soluções possíveis para as condições de um método de
quarta ordem e, em particular, o famoso método clássico de quarta ordem já apresentado
1.3. Métodos “Multistep” 21
anteriormente nesse capítulo
0
1
2
1
2
1
2
0 1
2
1 0 0 1
1
6
1
3
1
3
1
6
1.3 Métodos “Multistep”
Os métodos one–step vistos na Seção 1.1 utilizam somente informações do ponto xi
a fim de predizer o valor da variável dependente yi+1 no ponto futuro xi+1. Os métodos
multistep, que veremos a seguir, baseiam-se no fato de que uma vez iniciado o processo
numérico nós temos ao nosso dispor informações importantes provenientes dos pontos
anteriores. A curva obtida a partir da conexão destes valores prévios nos fornece infor-
mações a respeito da trajetória da solução. Tal fato é explorado pelos métodos que serão
apresentados nesta seção.
1.3.1 Métodos do Tipo Preditor–Corretor
Iremos mostrar, agora, como as fórmulas do tipo Preditor–Corretor podem ser deriva-
das matematicamente. Esta derivação é interessante no sentido de que serão empregados
conceitos como a integração numérica e a resolução de equações diferenciais ordinárias.
Além do mais, este procedimento pode ser empregado na obtenção de métodos multistep
de alta ordem e na estimativa dos erros associados a estes métodos.
O nosso ponto de partida está baseado na resolução da seguinte equação diferencial
ordinária
dy
dx
= f(x, y)
Sabemos que a solução desta equação pode ser obtida mediante a sua integração entre os
pontos i e i+ 1 ∫ yi+1
yi
dy =
∫ xi+1
xi
f(x, y) dx
onde a primeira integral pode ser facilmente avaliada:
yi+1 = yi +
∫ xi+1
xi
f(x, y) dx (1.17)
Esta equação nos fornece a solução da equação diferencial ordinária caso a integral, do lado
direito do sinal de igualdade, possa ser avaliada e o valor futuro da variável dependente
yi+1 pode ser determinado a partir do seu valor precedente yi e da resolução da integral
de f(x, y).
22 Capítulo 1. Equações Diferenciais Ordinárias
Uma primeira possibilidade de avaliarmos a integral em (1.17) pode ser dada pela
aplicação de uma integração numérica como, por exemplo, a regra do trapézio, ou seja:∫ xi+1
xi
f(x, y) dx =
f(xi, yi) + f(xi+1, yi+1)
2
∆x
Após substituição deste resultado em (1.17) obtemos:
yi+1 = yi +
f(xi, yi) + f(xi+1, yi+1)
2
∆x
que corresponde à equação corretora para o método de Heun. Do Cálculo Numérico[1]
sabemos que o erro de truncamento associado à regra do trapézio é dado por:
Ec = − 1
12
∆x3f ′′(ξc) = − 1
12
∆x3y′′′(ξc)
onde o subscrito c indica o erro do corretor e ξc encontra–se no intervalo determinado por
xi e xi+1.
Um procedimento similar pode ser empregado na obtenção do passo preditor, a partir
de uma integração entre os limites i− 1 e i+ 1:∫ yi+1
yi−1
dy =
∫ xi+1
xi−1
f(x, y) dx
e como resultado desta integração nós obtemos a forma da equação preditora:
yi+1 = yi−1 +
∫ xi+1
xi−1
f(x, y) dx (1.18)
Agora, ao invés de empregarmos uma fórmula fechada1, na avaliação da integral em (1.18),
vamos empregar a primeira fórmula aberta de Newton–Cotes, que é conhecida como o
método do ponto médio ∫ xi+1
xi−1
f(x, y) dx = 2∆xf(xi, yi)
cujo erro de truncamento é [1]:
Ep =
1
3
∆x3f ′′(ξp) =
1
3
∆x3y′′′(ξp)
sendo que o subscrito p indica o erro do preditor e ξp está contido no intervalo xi−1 a xi+1.
Após substituição da integral na Eq. (1.18) pela sua aproximação chegamos à forma
final da equação preditora
yi+1 = yi−1 + 2∆xf(xi, yi)
Vamos resumir, abaixo, as equações preditora e corretora associadas a este método,
cujos os erros local de truncamento são de O (∆x3):
1Uma fórmula fechada é aquela onde os pontos correspondentes aos limites de integração são co-
nhecidos, em oposição ao caso da fórmula aberta, onde os limites de integração vão além dos pontos
conhecidos.
1.3. Métodos “Multistep” 23
Preditor
y0i+1 = y
m
i−1 + 2∆xf(xi, y
m
i )
Corretor
yji+1 = y
m
i +
f(xi, y
m
i ) + f(xi+1, y
j−1
i+1 )
2
∆x
onde o sobrescrito j denota que a equação corretora é aplicada iterativamente de j = 1
até j = m para que possamos obter soluções refinadas. Portanto, os valores ymi e ymi−1
correspondem aos valores finais do processo corretivo.
Conforme podemos observar a equação preditora depende do conhecimento prévio do
valor da variável dependente yi−1 a fim de que ela possa ser empregada. Esta informação
não está normalmente disponível em um típico problema de valor inicial. Por este motivo
este método é chamado de non–self–starting Heun [1].
Caso os passos preditor e corretor de um método multistep apresentem a mesma ordem
do erro de truncamento, uma estimativa do erro local de truncamento pode ser obtida
durante o processo computacional de obtenção da solução numérica. Esta informação
pode, então, ser empregada como um critério de ajuste do incremento.
Devido à existência do erro de truncamento sabemos que a solução numérica é uma
aproximação da solução exata (yex = yi+1 +Et) da equação diferencial ordinária, ou seja:
yex = y
0
i+1 +
1
3
∆x3y′′′(ξp) (1.19)
e
yex = y
m
i+1 −
1
12
∆x3y′′′(ξc) (1.20)
Subtraindo a Eq. (1.19) da Eq. (1.20) chegamos ao seguinte resultado
0 = ymi+1 − y0i+1 −
5
12
∆x3y′′′(ξ)
onde ξ encontra-se agora no intervalo compreendido entre xi−1 e xi+1. Esta equação pode
ainda ser reescrita na forma:
y0i+1 − ymi+1
5
= − 1
12
∆x3y′′′(ξ)
Deste resultado podemos notar que ele é idêntico à expressão do erro de truncamento
do corretor, Ec, a menos do argumento da derivada terceira. Assumindo que a derivada
terceira de y não varie de forma apreciável sobre o intervalo entre xi−1 e xi+1, podemos
considerar como sendo válida a seguinte igualdade:
Ec ∼= −y
m
i+1 − y0i+1
5
24 Capítulo 1. Equações Diferenciais Ordinárias
e a partir desta relação nós podemos estimar o erro de truncamento com base em duas
quantidades que já são calculadas pelo método numérico. Assim sendo, esta estimativa
pode ser empregada como um critério para ajustarmos o incremento ∆x de modo que o
erro fique dentro de uma faixa de valores aceitáveis.
As derivações dos erros Ep e Ec podem ser generalizadas uma vez conhecidos os erros
de truncamento das fórmulas de integração empregadas na avaliação das integrais do passo
preditor:
yex = y
0
i+1 +
ηp
δp
∆xn+1y(n+1)(ξp)
e do passo corretor:
yex = y
m
i+1 −
ηc
δc
∆xn+1y(n+1)(ξc)
A partir do conhecimento destes valores os erros de truncamento do passo preditor
pode ser estimado como sendo dado por [1]:
Ep =
ηpδc
ηcδp + ηpδc
(
ymi − y0i
)
(1.21)
enquanto que no caso do passo corretor:
Ec ∼= − ηcδp
ηcδp + ηpδc
(
ymi+1 − y0i+1
)
(1.22)
1.3.2 Métodos de Mais Alta Ordem
Do desenvolvimento precedente fica claro que uma estratégia óbvia para que possamos
desenvolver métodos multistep de mais alta ordem passa pelo emprego de fórmulas de
integração de mais alta ordem para os passos preditor e corretor. Os métodos deste tipo
são compostos de um método explícito (preditor) e de outro implícito (corretor). Antes de
prosseguirmos com a apresentação destes métodos de mais alta ordem, vamos relembrar
as fórmulas de integração de Newton–Cotes e de Adams. As fórmulas de integração de
Newton–Cotes, para n pontos igualmente espaçados (n+1 segmentos), podem ser do tipo
aberta:
yi+1 = yi−n +
∫ xi+1
xi−n
f(x, y) dx
ou fechada, para n segmentos e n+ 1 pontos,
yi+1 = yi−n+1 +
∫ xi+1
xi−n+1
f(x, y) dx
onde estas integrais são aproximadas por um desenvolvimento do tipo:
In(f) = ∆x
n−1∑
k=0
wkf(xi−k, yi−k) +O
(
∆xp+2
)
1.3. Métodos “Multistep” 25
para a forma aberta e
In(f) = ∆x
n−1∑
k=−1
wkf(xi−k, yi−k) +O
(
∆xp+2
)
para a forma fechada, sendo n o número de segmentos empregados e os valores dos coe-
ficientes wk (k = −1, 0, 1, . . . , n) podem ser encontrados, por exemplo, na referência [1].
Nestas equações temos que p = n + 1 se n for par e p = n caso n seja ímpar. Portanto,
as fórmulas gerais de Newton-Cotes são dadas por
yi+1 = yi−n + ∆x
n−1∑
k=0
wkf(xi−k, yi−k) +O
(
∆xp+2
)
(1.23)
para a forma aberta e
yi+1 = yi−n+1 + ∆x
n−1∑
k=−1
wkf(xi−k, yi−k) +O
(
∆xp+2
)
(1.24)
para a forma fechada.
O outro tipo de fórmula de integração que é com frequência empregado na obtenção
de algoritmos multistep voltados para a resolução de equações diferenciais ordinárias é a
fórmula aberta de Adams–Bashforth. Uma forma de derivarmos estas fórmulas tem como
ponto de partida um desenvolvimento em série de Taylor em torno do ponto xi:
yi+1 = yi + fi∆x+ f
′
i
∆x2
2
+ f ′′i
∆x3
6
+ . . . = yi + ∆x
(
fi + f
′
i
∆x
2
+ f ′′i
∆x2
6
+ . . .
)
Empregando uma aproximação do tipo diferença recuada na aproximação de f ′i , este
resultado pode ser reescrito como:
yi+1 = yi + ∆x
{
fi +
∆x
2
[
fi − fi−1
∆x
+ f ′′i
∆x
2
+O (∆x2)]+ f ′′i ∆x26 + . . .
}
ou ainda,
yi+1 = yi + ∆x
(
3
2
fi − 1
2
fi−1
)
+
5
12
∆x3f ′′i +O
(
∆x4
)
Outras expressões de mais alta ordem podem ser obtidas se empregarmos aproximações
de mais alta ordem para f ′i . A fórmula aberta geral de Adams de ordem n pode ser
expressa como
yi+1 = yi + ∆x
n−1∑
k=0
βkfi−k +O
(
∆xn+1
)
(1.25)
Os valores dos coeficientes βk podem ser encontrados, para diferentes valores de n, na
referência [1].
26 Capítulo 1. Equações Diferenciais Ordinárias
Um procedimento similar pode ser feito para obtermos a expressão geral para a fórmula
fechada, conhecida como a fórmula de Adams–Moulton. Neste caso, nosso ponto de
partida será um desenvolvimento em série de Taylor em torno do ponto xi+1
yi = yi+1 − fi+1∆x+ f ′i+1
∆x2
2
− f ′′i+1
∆x3
6
+ . . .
onde devemos aproximar agora o valor da derivada f ′i+1. Deste modo obtemos a fórmula
fechada geral de Adams de ordem n
yi+1 = yi + ∆x
n−1∑
k=0
βkfi+1−k +O
(
∆xn+1
)
(1.26)
com os valores de βk também fornecidos em [1].
1.3.3 Método de Milne
Um dos mais comuns métodos do tipo multistep é conhecido como o método de Milne
que emprega fórmulas de integração de Newton–Cotes. Neste método uma formulação a
três pontos (4 segmentos) do tipo aberta de Newton–Cotes (1.23), n=3, é usadapara o
passo preditor:
y0i+1 = y
m
i−3 +
4
3
(
2fmi − fmi−1 + 2fmi−2
)
∆x
enquanto que uma fórmula fechada com dois segmentos (três pontos) de Newton–Cotes (1.24)
(regra de Simpson 1/3), n=2, é empregada para o passo corretor:
yji+1 = y
m
i−1 +
1
3
(
fmi−1 + 4f
m
i + f
j−1
i+1
)
∆x
Para este método temos que os erros de truncamento preditor e corretor são de O(∆x5) e
ηp = 14, δp = 45, ηc = 1 e δc = 90 e, a partir das fórmulas gerais, os erros de truncamento
associados a estes dois passos podem ser determinados pelas Eqs. (1.21) e (1.22)
Ep =
28
29
(
ymi − y0i
)
Ec ∼= − 1
29
(
ymi+1 − y0i+1
)
1.3.4 Método de Adams de Quarta Ordem
Outro método multistep popular pode ser implementado a partir do emprego de fór-
mulas do tipo Adams–Bashforth (1.25) de quarta ordem (n = 4) para o passo preditor
y0i+1 = y
m
i +
(
55
24
fmi −
59
24
fmi−1 +
37
24
fmi−2 −
9
24
fmi−3
)
∆x
1.3. Métodos “Multistep” 27
e, para o passo corretor, uma fórmula de quarta ordem (n=4) do tipo Adams–Moulton (1.26)
yji+1 = y
m
i +
(
9
24
f j−1i+1 +
19
24
fmi −
5
24
fmi−1 +
1
24
fmi−2
)
∆x
Os erros de truncamento associados a estes dois passos são de O(∆x5) e dados por [1]
Ep =
251
270
(
ymi − y0i
)
Ec = − 19
270
(
ymi+1 − y0i+1
)
uma vez que ηp = 251, δp = 720, ηc = 19 e δc = 720 para o método em questão.
À primeira vista pode parecer que o método de Milne é melhor que o método de quarta
ordem de Adams, visto que o primeiro necessita de menos avaliações da função f(x, y) e
apresenta um erro do passo corretor bem inferior ao do método de Adams. Embora essa
conclusão aplique–se a muitos casos, existem problemas para os quais o método de Milne
apresenta resultados numericamente inaceitáveis. Esta fraca performance se deve a uma
instabilidade do método de Milne, sendo oriunda especificamente do passo corretor [1].
1.3.5 Método de Hamming
O método de Hamming vem a ser uma alternativa estável ao método de Milne. Ele
emprega o mesmo passo preditor que o do método de Milne:
y0i+1 = y
m
i−3 +
4
3
(
2fmi − fmi−1 + 2fmi−2
)
∆x+
14
45
∆x5y(5)
enquanto que o passo corretor do método de Milne é substituído por uma versão estável
que apresenta um erro de truncamento de mesma ordem de grandeza que a do método de
Milne [1]:
yji+1 =
1
8
[
9ymi − ymi−2 + 3
(
f j−1i+1 + 2f
m
i − fmi−1
)
∆x
]− 1
40
∆x5y(5)
Destes valores dos erros de truncamento e das fórmulas gerais (1.21) e (1.22) chegamos
aos seguintes resultados para as estimativas do erro de truncamento
Ep =
112
121
(
ymi − y0i
)
Ec ∼= − 9
121
(
ymi+1 − y0i+1
)
28 Capítulo 1. Equações Diferenciais Ordinárias
1.4 Convergência, Consistência e Estabilidade
Até o presente momento nós apresentamos alguns métodos numéricos do tipo one-step
emultistep e destacamos o fato de que quanto maior a ordem do erro de truncamento maior
será a acurácia destes métodos. Entretanto, conforme nós vimos nas Seções 1.1.1.2 e 1.3.4,
estes métodos numéricos podem apresentar resultados que não são acurados. Nesta seção
nós vamos introduzir três conceitos que são fundamentais para que possamos entender
porquê alguns destes métodos podem falhar. Estas propriedades que estão associadas aos
diferentes métodos do tipo diferenças são a convergência, a consistência e a estabilidade
e elas serão introduzidas separadamente para os métodos do tipo one-step e multistep.
1.4.1 Métodos One-Step
A noção de convergência será introduzida para o problema de valor inicial
dy
dx
= f(x, y)
para uma condição inicial dada por y(x0) = y0 e vamos assumir que a função f(x, y)
satisfaz a condição de Lipschitz em y. Uma função é dita satisfazer a condição de Lipschitz
na variável y em um conjunto D ⊂ R2 se
|f(x, y)− f(x, z)| ≤ L |y − z|
para alguma constante L > 0 e x, y e z ∈ D [3].
1.4.1.1 Convergência
Um método numérico é dito convergir se [3]
lim
i→∞
|y(x)− yi| = 0
para todo x0 ≤ x ≤ X, onde x = x0 + i∆x e yi é a aproximação numérica correspondente
à solução exata y(x).
Entretanto, esta definição não pode ser aplicada diretamente a fim de estabelecermos
a convergência de um método numérico, uma vez que nós não conhecemos, em geral, a
solução exata do problema y(x). A convergência será comprovada indiretamente mediante
a introdução dos conceitos de consistência e de estabilidade.
1.4.1.2 Consistência
Um método numérico é dito ser consistente se
lim
∆x→0
Et
∆x
= 0
onde Et representa o erro de truncamento local do método do tipo diferenças que nós
queremos estudar. Esta condição implica no fato de que o erro de truncamento local de
ummétodo do tipo diferenças deve ser pelo menos deO (∆x2) para que ele seja consistente.
1.4. Convergência, Consistência e Estabilidade 29
Conforme nós já vimos, todos os métodos one-step podem ser postos na forma geral
yi+1 = yi + φ(xi, yi,∆x)∆x
onde φ(x, y,∆x) é conhecida com a função incremento. A partir do desenvolvimento em
série de Taylor, em torno do ponto (xi, yi), do lado esquerdo desta igualdade nós obtemos
yi + f(xi, yi)∆x+ f
′
(xi, yi)
∆x2
2
+ . . . = yi + φ(xi, yi,∆x)∆x
Este resultado pode ainda ser rearrumado de modo a obtermos
φ(xi, yi,∆x) = f(xi, yi) +O(∆x)
e um método do tipo one-step será consistente caso
lim
∆x→0
φ(x, y,∆x) = f(x, y)
Uma outra maneira de vermos este resultado seria a de reescrevermos o desenvolvi-
mento em série de Taylor na forma:
dy
dx
∣∣∣
i
∆x = φ(x, y,∆x)∆x− d
2y
dx2
∣∣∣
i
∆x2
2
+ . . .
ou ainda,
dy
dx
∣∣∣
i
= φ(x, y,∆x)− d
2y
dx2
∣∣∣
i
∆x
2
+ . . .
e observarmos que para que esta equação seja consistente com a equação diferencial ordi-
nária modelo, em cada ponto da malha, a seguinte condição deve ser verificada à medida
que o incremento ∆x tende a zero:
lim
∆x→0
φ(x, y,∆x) = f(x, y)
1.4.1.3 Estabilidade
Um método numérico do tipo diferenças é dito ser estável para um problema de valor
inicial, do tipo apresentado no início desta seção, se existir um ∆x0 > 0 tal que uma
variação na condição inicial, por uma quantidade fixa, produz uma variação limitada na
solução numérica para todo 0 < ∆x ≤ ∆x0.
Antes de prosseguirmos com o nosso estudo da estabilidade, vamos enunciar o seguinte
teorema que será fundamental para que possamos garantir que um método numérico do
tipo diferenças seja convergente.
Teorema 2 (Asaithambi [3]) Um método do tipo diferenças consistente é convergente
se e somente se ele for estável.
30 Capítulo 1. Equações Diferenciais Ordinárias
Portanto, a convergência estará assegurada se o método numérico for consistente e
estável.
Uma outra maneira de enunciarmos a estabilidade de um método numérico do tipo
one-step é fornecida se para duas soluções numéricas diferentes yi e zi, correspondentes à
aproximação numérica de y(xi), nós verificarmos a seguinte desigualdade
|yi − zi| ≤ K |y0 − z0|
para uma constante K > 0 e para todo ∆x tal que 0 < ∆x ≤ ∆x0 para algum ∆x0 > 0,
e para todo x0 ≤ xi ≤ X. Nesta desigualdade y0 e z0 correspondem aos valores iniciais.
A noção de estabilidade também pode ser relacionada com a regularidade da função
incremento φ(x, y,∆x) mediante o teorema
Teorema 3 (Asaithambi [3]) Um método numérico do tipo one-step é estável se ele
for regular.
Em seguida, veremos em que condições o método numérico pode ser considerado re-
gular.
Da forma geral dos métodos one-step podemos dizer que ele será regular se a função
incremento satisfaz a condição de Lipschitz em y no domínio D = {(x,w)|x0 ≤ x ≤
X,−∞ < w < ∞, 0 < ∆x < ∆x0} para algum ∆x0 > 0 de modo que a desigualdade
abaixo é verificada
|φ(x, y,∆x)− φ(x, z,∆x)| ≤ M|y − z|
para y e z ∈ D e M é chamada de constante de Lipschitz para a função incrementoφ(x, y,∆x).
A prova deste teorema é relativamente simples e podemos demonstrá-la se nós consi-
derarmos, mais uma vez, duas aproximações numéricas yi e zi geradas a partir da forma
geral, de modo que
|yi+1 − zi+1| ≤ |yi − zi|+ ∆x|φ(xi, yi,∆x)− φ(xi, zi,∆x)|
Como estamos considerando que o método numérico é regular
|yi+1 − zi+1| ≤ |yi − zi|+ ∆xM|yi − zi|
ou ainda
|yi+1 − zi+1| ≤ (1 + ∆xM)|yi − zi|
Por indução podemos reescrever este resultado como
|yi+1 − zi+1| ≤ (1 + ∆xM)i+1|y0 − z0|
(1 + ∆xM)|yi − zi| ≤ (1 + ∆xM)i+1 |y0 − z0|
|yi − zi| ≤ (1 + ∆xM)i |y0 − z0|
1.4. Convergência, Consistência e Estabilidade 31
Agora, do desenvolvimento em série de Taylor da função exponencial ex
ex = 1 + x+
x2
2!
+
x3
3!
+ . . .
nós podemos substituir o termo entre parênteses na desigualdade acima
|yi − zi| ≤
(
e∆xM
)i |y0 − z0|
|yi − zi| ≤ eM(i∆x) |y0 − z0|
|yi − zi| ≤ eM(xi−x0) |y0 − z0|
|yi − zi| ≤ eM(X−x0) |y0 − z0|
que para K = eM(X−x0) representa a condição para que o método numérico seja estável.
Finalmente, podemos enunciar o seguinte teorema
Teorema 4 (Asaithambi [3]) Um método numérico one-step regular é convergente se
e somente se ele for consistente.
1.4.2 Métodos Multistep
Visto os conceitos de convergência, consistência e estabilidade para os métodos do tipo
one-step, vamos aplicar estas mesmas propriedades ao caso dos métodos multistep.
1.4.2.1 Convergência
A mesma definição de convergência, aplicada aos métodos one-step, também se aplica
ao caso dos métodos multistep. Assim como no caso anterior, a fim de mostrarmos que
um método numérico do tipo diferenças é convergente, devemos estabelecer que ele é
consistente e estável [3].
1.4.2.2 Consistência
Os métodos multistep, a exemplo dos métodos one-step, também podem ser postos
numa forma geral. A forma mais geral de um método (k + 1) multistep é dada por [3]
yi+1 =
k∑
j=0
αjyi−j + ∆x
k∑
j=−1
βjfi−j para i ≥ k
onde para β−1 6= 0 nós obtemos um método implícito e para β−1 = 0 um método explícito.
Nesta equação os coeficientes α e β são escolhidos de modo a obtermos um método com
a mais alta ordem de acurácia possível. Entretanto, a escolha destes coeficientes baseada
somente na ordem de acurácia do método pode nos levar a métodos que não são estáveis.
32 Capítulo 1. Equações Diferenciais Ordinárias
Como nós sabemos que a solução exata menos a solução numérica é igual ao erro de
truncamento, Et = yex(xi+1)− yi+1,
Et = yex(xi+1)−
[
k∑
j=0
αjyi−j + ∆x
k∑
j=−1
βjfi−j
]
= yex(xi+1)−
[
k∑
j=0
αjyi−j + ∆x
k∑
j=−1
βjy
′
i−j
]
Substituindo a solução aproximada pela solução exata e expandindo em série de Taylor,
em torno do ponto xi, o lado direito desta igualdade [3]
Et ≈
m∑
r=0
y(r)ex
∆xr
r!
−
[
k∑
j=0
αj
m∑
r=0
y(r)ex
(−j∆x)r
r!
+ ∆x
k∑
j=−1
βj
m−1∑
r=0
y(r+1)ex
(−j∆x)r
r!
]
ou ainda, agrupando os termos semelhantes:
Et ≈
[
1−
k∑
j=0
αj
]
yex +
[
1−
(
−
k∑
j=0
jαj +
k∑
j=−1
βj
)]
∆x
1!
y
′
ex
+
m∑
r=2
[
1−
k∑
j=0
αj(−j)r − r
k∑
j=−1
βj(−j)r−1
]
∆xr
r!
y(r)ex
Deste resultado e da condição para que um método numérico seja consistente,
lim
∆x→0
Et
∆x
= 0
vemos claramente que as seguintes igualdades devem ser verificadas para que um método
multistep seja consistente:
k∑
j=0
αj = 1
−
k∑
j=0
jαj +
k∑
j=−1
βj = 1
1.4.2.3 Estabilidade
Finalizando esta seção iremos abordar o problema da estabilidade dos métodos mul-
tistep. As condições de estabilidade que devem ser verificadas pelos métodos multistep
serão introduzidas a partir do comportamento da solução numérica do seguinte problema
de valor de contorno trivial
dy
dx
= 0
1.4. Convergência, Consistência e Estabilidade 33
tendo por condição inicial y(0) = c e é evidente, para este problema, que f(x, y) = 0.
Caso a condição inicial seja introduzida sem nenhum erro de arredondamento, qualquer
método multistep forneceria a solução exata deste problema para todo i ≥ k. Por outro
lado, se a condição inicial estiver contaminada pelos erros de arredondamento, poderão
ocorrer desvios nas soluções numéricas. São estes desvios que nós queremos estudar aqui.
Da equação geral dos métodos multistep aplicada a este problema específico nós obtemos:
yi+1 =
k∑
j=0
αjyi−j para i ≥ k
ou ainda, numa forma equivalente (i = r + k),
yr+(k+1) =
k∑
j=0
αjyr+(k−j) para r ≥ 0
Usando o operador deslocamento E [3] a equação acima pode ser reescrita como
Ek+1yr −
k∑
j=0
αjE
k−jyr = p(E)yr = 0
onde o polinômio característico, p(λ), associado a esta equação de diferenças é dado por
p(λ)yr =
[
λk+1 −
k∑
j=0
αjλ
k−j
]
yr = 0
e cujas raízes (que podem ser complexas) são λ1, λ2, . . . , λk+1. As condições de estabi-
lidade dos métodos multistep serão estabelecidas em função destas raízes do polinômio
característico, de modo que os erros de arredondamento não cresçam exponencialmente [3].
Para uma equação linear de diferenças de ordem k do tipo
αkyi+k + αk−1yi+k−1 + . . .+ α1yi+1 + α0yi = 0
a solução geral desta equação é dada por:
yi =
k∑
j=1
pj(i)λ
i
j
onde λj é uma raíz do polinômio característico e o polinômio pj(.) possui um grau a menos
que a multiplicidade das raízes do polinômio característico.
Portanto, caso todas as raízes fossem distintas nós poderíamos escrever a solução do
problema de valor inicial, para λ1 ∼= c, como [3]
yi = c+
k∑
j=2
γjλ
i
j
34 Capítulo 1. Equações Diferenciais Ordinárias
onde o somatório pode ser visto como os erros de arredondamento associados à solução
aproximada. Deste resultado vemos que o erro de arredondamento irá crescer exponenci-
almente a menos que |λj| ≤ 1 para j = 2, 3, . . . , k.
O próximo teorema generaliza este resultado e impõe as condições que devem ser
verificadas para que o método numérico seja estável
Teorema 5 (Asaithambi [3]) O método numérico multistep dado pela sua forma geral
é estável se e somente se ele satisfizer a condição da raíz:
|λj| ≤ 1, j = 1, 2, . . . , k
|λj| = 1 ⇒ dp(λj)
dλ
6= 0
A primeira condição requer que todas as raízes do polinômio característico estejam
contidas no círculo {λ | |λ| ≤ 1} no plano complexo. A segunda condição estabelece que
todas as raízes que se encontram na fronteira do círculo unitário devem ser simples.
Um método multistep é dito ser fortemente estável se ele verifica a condição da raíz
e λ = 1 é a única raíz característica de módulo igual a um. Ele é dito ser fracamente
estável se satisfaz a condição da raíz e tem mais de uma raíz característica distinta com
magnitude igual a um. Por outro lado, ele será dito instável se ele não satisfizer a
condição da raíz.
O último teorema a ser enunciado estabelece as condições para que tenhamos um
método multistep convergente
Teorema 6 (Asaithambi [3]) Um método multistep é convergente se e somente se ele
for consistente e satisfizer as condição da raíz.
Capítulo 2
Equações Diferenciais Parciais
Como veremos mais adiante, as equações que governam o escoamento de fluidos e o
transporte de massa e energia são equações diferenciais parciais, contendo termos em de-
rivadas primeira e segunda nas coordenadas espaciais e, geralmente, em derivada primeira
no tempo. Nestas equações as derivadas no tempo aparecem como termos lineares, mas
frequentemente as derivadas espaciais estão associadas a termos não-lineares.
Neste capítulo, veremos como classificar as equações diferenciais parciais em elípticas,
parabólicas e hiperbólicas. Em seguida elas serão examinadas segundo as suas caracterís-
ticas matemáticas e físicas.
2.1 Algumas Equações Diferenciais de Interesse
No nosso curso iremos nos concentrar no estudo da equação de transporte linear.
Esta equação inclui os mecanismos de transportedo tipo convectivo e difusivo e vários
problemas de interesse nas engenharias são governados por este tipo de equação.
A título de exemplo vamos listar abaixo algumas equações diferenciais parciais impor-
tantes que governam vários problemas físicos que aparecem com frequência. Em alguns
casos específicos podemos obter uma solução analítica exata para estas equações [8].
1. A equação linear da onda de primeira ordem
∂ϕ
∂t
+ c
∂ϕ
∂x
= 0
governa a propagação de uma onda que se movendo da esquerda para a direita com
uma velocidade constante c.
2. A equação de Burger sem viscosidade
∂ϕ
∂t
+ ϕ
∂ϕ
∂x
= 0
que governa a propagação não–linear de uma onda simples em uma dimensão.
35
36 Capítulo 2. Equações Diferenciais Parciais
3. A equação de Burger
∂ϕ
∂t
+ ϕ
∂ϕ
∂x
= ν
∂2ϕ
∂x2
governa a propagação não–linear de uma onda simples em uma dimensão com dis-
sipação viscosa.
4. A equação de Poisson
∂2ϕ
∂x2
+
∂2ϕ
∂y2
= f(x, y)
que governa a distribuição de temperatura em um sólido com fontes de calor prescri-
tas pela função f(x, y). A equação de Poisson também determina o campo elétrico
em uma região contendo uma densidade de carga dada por f(x, y).
5. A equação de convecção–difusão
∂ϕ
∂t
+ u
∂ϕ
∂x
= α
∂2ϕ
∂x2
que representa o transporte, unidimensional, por convecção e difusão da quantidade
ϕ em uma região que encontra-se em movimento com velocidade u. A quantidade
α representa o coeficiente de difusão.
6. A equação de Korteweg–de Vries
∂ϕ
∂t
+ ϕ
∂ϕ
∂x
+ β
∂3ϕ
∂x3
= 0
que governa o movimento não–linear de ondas dispersivas.
7. A equação de Helmholtz
∂2ϕ
∂x2
+
∂2ϕ
∂y2
+ k2ϕ = 0
que governa o movimento de ondas harmônicas, sendo k o parâmetro de frequência.
Como exemplo temos a propagação de ondas acústicas.
8. A equação biharmônica
∂4ϕ
∂x4
+
∂4ϕ
∂y4
= 0
que determina as linhas de corrente para escoamentos com muito baixo número de
Reynolds, governado pela equação de Stokes.
2.2. Classificação 37
9. Equação do telégrafo
∂2ϕ
∂t2
+ a
∂ϕ
∂t
+ b ϕ = c2
∂2ϕ
∂x2
que governa a transmissão de impulsos elétricos num fio longo com uma certa dis-
tribuição de capacitância, indutância e resistência. As aplicações desta equação
também incluem o movimento de uma corda com uma força de amortecimento pro-
porcional à velocidade e a condução de calor com uma velocidade de propagação
térmica finita.
2.2 Classificação
Uma simples classificação das equações diferenciais parciais é possível no caso das
equações lineares de segunda ordem,
A
∂2u
∂x2
+B
∂2u
∂x∂y
+ C
∂2u
∂y2
+D
∂u
∂x
+ E
∂u
∂y
+ Fu+G = 0
onde A, B, . . ., G são constantes. Três categorias podem ser distinguidas:
i) elíptica: B2 − 4AC < 0
ii) parabólica: B2 − 4AC = 0
iii) hiperbólica: B2 − 4AC > 0
Como podemos observar, a classificação depende somente dos coeficientes dos termos de
derivadas de maior ordem nas variáveis independentes.
Caso os coeficientes A, B, . . ., G sejam funções de x, y, ∂u/∂x ou ∂u/∂y, a classificação
pode ser aplicada, contanto que A, B e C tenham uma interpretação local. Isto implica
no fato de que a classificação das equações pode mudar em diferentes partes do domínio
de resolução.
As diferentes categorias de EDPs (Equações Diferenciais Parciais) podem ser, gros-
seiramente, associadas aos diferentes tipos de escoamentos. Em geral, problemas de-
pendentes do tempo estão associados a EDPs do tipo parabólico ou hiperbólico. EDPs
parabólicas governam escoamentos contendo mecanismos de dissipação, como a dissipação
viscosa e térmica. Neste caso, os gradientes decrescerão para valores crescentes do tempo,
se as condições de contorno não forem dependentes do tempo. Não havendo mecanismos
de dissipação, a solução terá uma amplitude constante no caso de uma EDP linear e
poderá aumentar caso ela seja não-linear. Esta solução é governada por EDPs do tipo
hiperbólicas. As EDPs elípticas estão normalmente associadas a problemas em regime
permanente e de equilíbrio. Entretanto, alguns problemas em regime permanente são
descritos por EDPs parabólicas (escoamento do tipo camada limite) e EDPs hiperbólicas
(escoamento supersônico não-viscoso).
38 Capítulo 2. Equações Diferenciais Parciais
Em se tratando de sistemas de equações diferenciais, a classificação não pode ser feita
mediante uma simples inspeção das equações. Usualmente é necessário examinarmos
as características associadas às equações, a fim de determinarmos de maneira correta a
classificação.
2.2.1 Classificação pelas Características
A resolução numérica de uma equação diferencial parcial requer que sejam fornecidas
condições auxiliares (inicial e de contorno) adequadas. A especificação destas condições
auxiliares pode ser facilitada com a introdução do conceito das características. Em se
tratando de um problema bidimensional, caso existam características reais, elas serão
representadas por curvas no interior do domínio de resolução. Ao longo destas curvas
características são propagadas, para o interior do domínio, as informações provenientes
das condições auxiliares. No caso de problemas hiperbólicos e para condições auxiliares
contendo descontinuidades, estas condições acarretam na existência de soluções que serão
igualmente descontínuas.
Tomemos como exemplo o caso de uma equação diferencial parcial de primeira ordem
na forma:
A
∂u
∂t
+B
∂u
∂x
= C (2.1)
na qual A, B e C podem ser funções de t, x e u, mas não podem ser funções das de-
rivadas da variável dependente u, pois elas podem ser descontínuas nas características.
Suponhamos, agora, que u seja especificado ao longo de uma curva L no plano t − x,
Fig. 2.1. Vamos também introduzir, neste plano, uma curva característica C que inter-
cepta a curva L no ponto p de coordenadas (xp,tp), conforme representado na Fig. 2.1.
Ao longo da curva característica C temos que [9]
du =
∂u
∂t
dt+
∂u
∂x
dx (2.2)
para dt e dx infinitesimais.
Combinando as Eqs. (2.1) e (2.2) de modo a eliminarmos o termo ∂u/∂t(
dx
dt
− B
A
)
∂u
∂x
=
du
dt
− C
A
Para que ∂u/∂x possa ser indeterminado ao longo da curva C, devemos fazer
dx
dt
=
B
A
de modo que esta equação define a curva característica C e vemos que
du
dt
=
C
A
ao longo da característica. A resolução numérica desta equação fornece o valor de u.
2.2. Classificação 39
x
t
C
L
p
p px
t
Figura 2.1: Curva características no plano t-x.
Da integração da equação que define a curva característica vemos que existe uma
infinidade de curvas características no plano t− x. Em particular, para B/A =constante,
temos que as características são retas e as suas equações são dadas por:
x(t)− x0 =
(
B
A
)
(t− t0)
onde as constantes x0 e t0 representam as coordenadas de um ponto qualquer sobre a
característica.
Finalizando, o conhecimento do valor de u no ponto p sobre a curva L, que intercepta
a curva característica C, permite que determinemos a variável dependente u mediante a
resolução de dois problemas de valor inicial (PVI):
dx
dt
=
B
A
x(tp) = xp
que fornece os valores de x e t de todos os pontos que se encontram sobre a característica
C que passa pelo ponto (tp,xp) e 
du
dt
=
C
A
u(tp) = up
que provê os valores de u ao longo da característica. Os valores de u que não se encontram
sobre esta característica devem ser determinados pela resolução dos mesmos problemas
de valor inicial, mas para novas condições iniciais.
40 Capítulo 2. Equações Diferenciais Parciais
Problemas para os quais mais de uma característica passam por um mesmo ponto
apresentam a formação de choques nestes pontos. O choque é resultante da propagação
de informações diferentes ou conflitantes do valor da variável dependente