apostila mned

Prévia do material em texto

Métodos Numéricos
para
Equações Diferenciais
Helio Pedro Amaral Souto
Nova Friburgo, 3 de Agosto de 2017
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 . . . . . . . . . . . . . . . . . 6
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 . . . . . . . . . . . 14
1.1.5.4 Sistema de Equações . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2 Formulação de Butcher . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.1 Condições de Ordem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.2.2 Alguns Métodos Explícitos . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.3 Métodos “Multistep” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.3.1 Métodos do Tipo Preditor–Corretor . . . . . . . . . . . . . . . . . . . . . 24
1.3.2 Métodos de Mais Alta Ordem . . . . . . . . . . . . . . . . . . . . . . . . 27
1.3.3 Método de Milne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
1.3.4 Método de Adams de Quarta Ordem . . . . . . . . . . . . . . . . . . . . 29
1.3.5 Método de Hamming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
1.4 Convergência, Consistência e Estabilidade . . . . . . . . . . . . . . . . . . . . . 30
i
ii CONTEÚDO
1.4.1 Métodos One-Step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
1.4.1.1 Convergência . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
1.4.1.2 Consistência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
1.4.1.3 Estabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
1.4.2 Métodos Multistep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.4.2.1 Convergência . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.4.2.2 Consistência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.4.2.3 Estabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2 Equações Diferenciais Parciais 39
2.1 Algumas Equações Diferenciais de Interesse . . . . . . . . . . . . . . . . . . . . 39
2.2 Classificação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.2.1 Classificação pelas Características . . . . . . . . . . . . . . . . . . . . . . 42
2.2.2 Equações a N Variáveis Independentes . . . . . . . . . . . . . . . . . . . 45
2.2.3 Transformação de Coordenadas . . . . . . . . . . . . . . . . . . . . . . . 46
2.2.4 Sistema de Equações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.2.5 Classificação pela Análise de Fourier . . . . . . . . . . . . . . . . . . . . 49
2.3 Problema Bem-Posto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.3.1 Condições de Contorno e Inicial . . . . . . . . . . . . . . . . . . . . . . . 51
2.4 Equações Diferenciais Hiperbólicas . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.4.1 Interpretação em Bases Físicas . . . . . . . . . . . . . . . . . . . . . . . . 53
2.4.2 Condições de Contorno e Inicial Apropriadas . . . . . . . . . . . . . . . . 56
2.5 Equações Diferenciais Parabólicas . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.5.1 Interpretação em Bases Físicas . . . . . . . . . . . . . . . . . . . . . . . . 59
2.5.2 Condições de Contorno e Inicial Apropriadas . . . . . . . . . . . . . . . . 59
2.6 Equações Diferenciais Elípticas . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
2.6.1 Interpretação em Bases Físicas . . . . . . . . . . . . . . . . . . . . . . . . 61
2.6.2 Condições de Contorno Apropriadas . . . . . . . . . . . . . . . . . . . . 61
3 Equações de Balanço 63
3.1 Equação da Continuidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.2 Equações do Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.3 Equação de Energia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4 Discretização 69
4.1 Discretização Espacial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.2 Discretização no Tempo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.3 Expansões em Séries de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.4 Técnica Geral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.5 Acurácia do Processo de Discretização . . . . . . . . . . . . . . . . . . . . . . . 74
4.6 Representação do Tipo Onda . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
CONTEÚDO iii
5 Convergência, Consistência e Estabilidade 81
5.1 Convergência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.1.1 Teorema de Lax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.2 Consistência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.2.1 O Esquema FTCS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.2.2 Esquema Totalmente Implícito . . . . . . . . . . . . . . . . . . . . . . . . 84
5.3 Estabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.3.1 Método da Matriz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.3.2 Método da Matriz: Condições de Contorno . . . . . . . . . . . . . . . . . 88
5.3.3 Método de Von Neumann . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.3.4 Método de Von Neumann: Esquema Geral . . . . . . . . . . . . . . . . . 90
5.4 Acurácia da Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.4.1 Extrapolação de Richardson . . . . . . . . . . . . . . . . . . . . . . . . . 92
6 Equação de Transporte Linear 95
6.1 Métodos Explícitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
6.2 Métodos Implícitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.3 Equação de Difusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.3.1 Equação de Difusão Unidimensional . . . . . . . . . . . . . . . . . . . . 97
6.3.2 Equação de Difusão em Regime Permanente . . . . . . . . . . . . . . . . 103
6.3.3 Equação de Difusão Bidimensional . . . . . . . . . . . . . . . . . . . . . 104
6.4 Métodos Multidimensionais de Fracionamento . . . . . . . . . . . . . . . . . . . 106
6.5 Equação de Advecção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.6 Dissipação e Dispersão Numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . 117
6.6.1 Análise de Fourier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.6.2 Equação Modificada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.7 Equação de Transporte Unidimensional . . . . . . . . . . . . . . . . . . . . . . . 123
6.8 Equação de Transporte Bidimensional . . . . . . . . . . . . . . . . . . . . . . . . 133
A O Método de Separação de Variáveis 137
B Algoritmo de Thomas 143
B.1 O Algoritmo de Thomas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
B.2 O Método TDMA Linha a Linha . . . . . . . . . . . . . . . . . . . . . . . . . . 145
B.3 Relaxação . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . 147
Bibliografia 149
Índice 151
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. . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.2 Domínio computacional R. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.3 Características para a equação da onda. . . . . . . . . . . . . . . . . . . . . . . 54
2.4 Propagação de uma onda senoidal. . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.5 Condições de contorno para o caso hiperbólico. . . . . . . . . . . . . . . . . . . 57
2.6 Condições auxiliares para o caso transiente. . . . . . . . . . . . . . . . . . . . . 57
2.7 Domínio computacional para uma EDP parabólica. . . . . . . . . . . . . . . . . 58
2.8 Decaimento exponencial de uma senóide. . . . . . . . . . . . . . . . . . . . . . 59
2.9 Solução da equação de Laplace. . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
2.10 Domínio computacional para uma EDP elíptica. . . . . . . . . . . . . . . . . . 61
4.1 Malha discreta uniforme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.1 Exemplo da representação matricial. . . . . . . . . . . . . . . . . . . . . . . . . 107
6.2 Ângulo de fase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
B.1 Método TDMA linha a linha. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
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 . . . . . . . . . . . . . . . . . 19
4.1 Comparação da acurácia dos diferentes esquemas. . . . . . . . . . . . . . . . . . 76
4.2 Relação entre as derivadas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
A.1 Separação da equação de Helmohtz. . . . . . . . . . . . . . . . . . . . . . . . . . 140
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, eletricidade, ter-
modinâ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
1
22
de Newton), temos a equação que governa a queda de um paraquedista
dv
dt
= g − c
m
v
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 contrapartida, 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 garantir-
mos 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 especi-
ficadas 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.
Capítulo 1. Equações Diferenciais Ordinárias 3Capítulo 1. Equações Diferenciais Ordinárias 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 obtida 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 expressã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çãodiferencial 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 1.1. Métodos “One–Step”4 1.1. Métodos “One–Step”
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 utili-
zados 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 proveniente 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
Capítulo 1. Equações Diferenciais Ordinárias 5Capítulo 1. Equações Diferenciais Ordinárias 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 quantifi-
carmos 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 comporta-
mento 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 função
f(x, y). É 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 ordi-
ná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.
6 1.1. Métodos “One–Step”6 1.1. Métodos “One–Step”
1.1.1.2 Estabilidade do Método de Euler
O conceito matemático de estabilidade do método numérico está associado ao crescimento
ilimitado da solução numérica ou, conforme veremos mais adiante, do erro introduzido nos
dados iniciais ou no processo de discretização. Neste caso, a solução numérica 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 incre-
mento ∆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 im-
plí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 destinados à
resolução de um problema de valor inicial do tipo:
dy
dx
= f(x, y), x > a
Capítulo 1. Equações Diferenciais Ordinárias 7Capítulo 1. Equações Diferenciais Ordinárias 7
para uma condição inicial y(a) = y0. Os leitores interessados devem se reportar às referências
[2, 3] para maiores detalhes. Este assunto será retomado ao final deste capítulo.
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 numérica. 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 in-
clusão de termos de mais alta ordem pode não ser trivial no caso de estarmos considerando
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 inde-
pendentes, sendo necessário o emprego da regra da cadeia na determinaçã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. Con-
forme 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ão8 1.1. Métodos “One–Step”8 1.1. Métodos “One–Step”
y
xx xi i+1
yi
y i+1
0
Figura 1.2: Método de Heun.
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)
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 ser combinadas 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
Capítulo 1. Equações Diferenciais Ordinárias 9Capítulo 1. Equações Diferenciais Ordinárias 9
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.
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 ne-
cessitamos 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 modificado.
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
10 1.1. Métodos “One–Step”10 1.1. Métodos “One–Step”
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
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 deO (∆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étodos 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 incli-
naçã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)
Capítulo 1. Equações Diferenciais Ordinárias 11Capítulo 1. Equações Diferenciais Ordinárias 11
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.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
)
+ · · ·
12 1.1. Métodos “One–Step”12 1.1. Métodos “One–Step”
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)
Finalmente, comparando termo a termo as expansões (1.5) e (1.6), vemos que as seguintes
condições devem serverificadas 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 infini-
dade de métodos de Runge–Kutta de segunda ordem, que fornecerão o mesmo resultado 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)
Capítulo 1. Equações Diferenciais Ordinárias 13Capítulo 1. Equações Diferenciais Ordinárias 13
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.
Método de Ralston (b2 = 2/3). Ralston (1962) e Ralston e Rabinowitz (1978) determina-
ram 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)
14 1.1. Métodos “One–Step”14 1.1. Métodos “One–Step”
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. Apresen-
taremos, 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
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)
...
...
Capítulo 1. Equações Diferenciais Ordinárias 15Capítulo 1. Equações Diferenciais Ordinárias 15
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 sistemas
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 faci-
litar 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
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 estabilidade
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
16 1.2. Formulação de Butcher16 1.2. Formulação de Butcher
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. 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 expansã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)),
yIV (x) = f ′′′fff + 2f ′′ff ′f + 2f ′f ′′ff + f ′f ′f ′f
= f ′′′(f, f, f) + 2f ′′ff ′f + 2f ′f ′′(f, f) + f ′f ′f ′f.
Capítulo 1. Equações Diferenciais Ordinárias 17Capítulo 1. Equações Diferenciais Ordinárias 17
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)) e 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 externas. Em se tra-
tando 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 repre-
sentações para o operador elementar diferencial F(ψ), considerando cada árvore ψ ∈ Ψ, sendo
Ψ o conjunto de todas as árvores com raízes.
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)
18 1.2. Formulação de Butcher18 1.2. Formulação de Butcher
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
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)
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.
Capítulo 1. Equações Diferenciais Ordinárias 19Capítulo 1. Equações Diferenciais Ordinárias 19
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.2 Alguns Métodos Explícitos
Na abordagem para a determinação dos coeficientes característicos, o usual é primeiro de-
terminarmos 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
20 1.2. Formulação de Butcher20 1.2. Formulação de Butcher
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
Das expansões (1.7), (1.8) e da Tabela 1.3 podemos verificar que os resultados determinados,
nos exemplos mostrados, correspondem de fato às séries serem coincidentes. Para um método
de segunda ordem teremos que
y(x0 + ∆x) = y0 +
∆xr(1)
σ(1)γ(1)
F(1)(y0) + ∆x
r(2)
σ(2)γ(2)
F(2)(y0)
= y0 +
∆x1
1 · 1f(y0) +
∆x2
1 · 2f
′f(y0)
e
y1 = y0 +
∆xr(1)
σ(1)
Φ(1)F(1)(y0) + ∆x
r(2)
σ(2)
Φ(2)F(2)(y0)
= y0 +
∆x1
1
(∑
bi
)
f(y0) +
∆x2
1
(∑
bici
)
f ′f(y0)
verificando a exatidão do sistema obtido anteriormente para a determinações dos coeficientes.
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
Capítulo 1. Equações Diferenciais Ordinárias 21Capítulo 1. Equações Diferenciais Ordinárias 21
bem como a formulação conhecida de Heun
0
1
3
1
3
2
3
0 2
3
1
4
0 3
4
De forma análoga ao caso do método anterior, obtemos das expansões (1.7), (1.8) e da
Tabela 1.3 os desenvolvimentos em série de Taylor da solução exata e da solução aproximada,
fornecida por um método de Runge-Kutta de terceira ordem. Em se tratando da solução exata,
temos que
y(x0 + ∆x) = y0 +
∆xr(1)
σ(1)γ(1)
F(1)(y0) + ∆x
r(2)
σ(2)γ(2)
F(2)(y0)
+
∆xr(3)
σ(3)γ(3)
F(3)(y0) + ∆x
r(4)
σ(4)γ(4)
F(4)(y0)
= y0 +
∆x1
1 · 1f(y0) +
∆x2
1 · 2f
′f(y0)
+
∆x3
2 · 3f
′′(f, f)(y0) +
∆x3
1 · 6f
′f ′f(y0)
e, para a solução numérica,
y1 = y0 +
∆xr(1)
σ(1)
Φ(1)F(1)(y0) + ∆x
r(2)
σ(2)
Φ(2)F(2)(y0)
+
∆xr(3)
σ(3)
Φ(3)F(3)(y0) + ∆x
r(4)
σ(4)
Φ(4)F(4)(y0)
= y0 +
∆x1
1
(∑
bi
)
f(y0) +
∆x2
1
(∑
bici
)
f ′f(y0)
+
∆x3
2
(∑
bic
2
i
)
f ′′(f, f)(y0) +
∆x3
1
(∑
biaijcj
)
f ′f ′f(y0)
e, mais uma vez, recuperamos o sistema de equações apresentado anteriormente. É importante
lembrar-nos de que estes resultados são válidos para os métodos explícitos, de modo que aij = 0
para j ≥ i.
22 1.2. Formulação de Butcher22 1.2. Formulação de Butcher
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)
•
•
•
b3a32c2 + b4a42c2 + b4a43c3=
1
6
(1.12)
•
•• •
@@�� b2c32 + b3c
3
3 + b4c
3
4 =
1
4
(1.13)
•
• •
•
��@@ b3c3a32c2 + 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. Escrevendo em termos dos desenvolvimentos em série de Taylor
y(x0 + ∆x) = y0 +
∆xr(1)
σ(1)γ(1)
F(1)(y0) + ∆x
r(2)
σ(2)γ(2)
F(2)(y0) + ∆x
r(3)
σ(3)γ(3)
F(3)(y0) + ∆x
r(4)
σ(4)γ(4)
F(4)(y0)
+
∆xr(5)
σ(5)γ(5)
F(5)(y0) + ∆x
r(6)
σ(6)γ(6)
F(6)(y0) + ∆x
r(7)
σ(7)γ(7)
F(7)(y0) + ∆x
r(8)
σ(8)γ(8)
F(8)(y0)
= y0 +
∆x1
1 · 1 f(y0) +
∆x2
1 · 2 f
′f(y0) +
∆x3
2 · 3 f
′′(f, f)(y0) +
∆x3
1 · 6 f
′f ′f(y0)
+
∆x4
6 · 4 f
′′′(f, f, f)(y0) +
∆x4
1 · 8 f
′′(f, f ′f)(y0) +
∆x4
2 · 12f
′f ′′(f, f)(y0) +
∆x4
1 · 24f
′f ′f ′f(y0)
Capítulo 1. Equações Diferenciais Ordinárias 23Capítulo 1. Equações Diferenciais Ordinárias 23
e, para a solução numérica,
y1 = y0 +
∆xr(1)
σ(1)
Φ(1)F(1)(y0) + ∆x
r(2)
σ(2)
Φ(2)F(2)(y0) + ∆x
r(3)
σ(3)
Φ(3)F(3)(y0)
+
∆xr(4)
σ(4)
Φ(4)F(4)(y0) + ∆x
r(5)
σ(5)
Φ(5)F(5)(y0)
+
∆xr(6)
σ(6)
Φ(6)F(6)(y0) + ∆x
r(7)
σ(7)
Φ(7)F(7)(y0)
+
∆xr(8)
σ(8)
Φ(8)F(8)(y0)
= y0 +
∆x1
1
(∑
bi
)
f(y0) +
∆x2
1
(∑
bici
)
f ′f(y0) +
∆x3
2
(∑
bic
2
i
)
f ′′(f, f)(y0)
+
∆x3
1
(∑
biaijcj
)
f ′f ′f(y0) +
∆x4
6
(∑
bic
3
i
)
f ′′′(f, f, f)(y0)
+
∆x4
1
(∑
biciaijcj
)
f ′′(f, f ′f)(y0) +
∆x4
2
(∑
biaijc
2
j
)
f ′f ′′(f, f)(y0)
+
∆x4
1
(∑
biaijajkck
)
f ′f ′f ′f(y0)
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 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
24 1.3. Métodos “Multistep”24 1.3. Métodos “Multistep”
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 informaçõ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 derivadas
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).
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
Capítulo 1. Equações Diferenciais Ordinárias 25Capítulo 1. Equações Diferenciais Ordinárias 25
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):
Preditor
y0i+1 = y
m
i−1 + 2∆xf(xi, y
m
i )
1Uma fórmula fechada é aquela onde os pontos correspondentes aos limites de integração são conhecidos, em
oposição ao caso da fórmula aberta, onde os limites de integração vão além dos pontos conhecidos.
26 1.3. Métodos “Multistep”26 1.3. Métodos “Multistep”
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 aproxi-
maçã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 podemosnotar 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
e a partir desta relação nós podemos estimar o erro de truncamento com base em duas quan-
tidades 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.
Capítulo 1. Equações Diferenciais Ordinárias 27Capítulo 1. Equações Diferenciais Ordinárias 27
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
)
para a forma aberta e
In(f) = ∆x
n−1∑
k=−1
wkf(xi−k, yi−k) +O
(
∆xp+2
)
28 1.3. Métodos “Multistep”28 1.3. Métodos “Multistep”
para a forma fechada, sendo n o número de segmentos empregados e os valores dos coeficientes
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ên-
cia [1].
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
+ . . .
Capítulo 1. Equações Diferenciais Ordinárias 29Capítulo 1. Equações Diferenciais Ordinárias 29
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, é usada para 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órmulas
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
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
30 1.4. Convergência, Consistência e Estabilidade30 1.4. Convergência, Consistência e Estabilidade
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 =
9
8
ymi −
1
8
ymi−2 +
3
8
(
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
)
1.4 Convergência, Consistência e Estabilidade
Até o presente momento nós apresentamos alguns métodos numéricosdo tipo one-step e
multistep 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.
Capítulo 1. Equações Diferenciais Ordinárias 31Capítulo 1. Equações Diferenciais Ordinárias 31
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(xi)− yi| = 0
para todo x0 ≤ x ≤ X, onde xi = 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 um método do
tipo diferenças deve ser pelo menos de O (∆x2) para que ele seja consistente.
Conforme nós já vimos, todos os métodos one-step podem ser postos na forma geral
yi+1 = yi + φ(x, y,∆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 + φ(x, y,∆x)∆x
32 1.4. Convergência, Consistência e Estabilidade32 1.4. Convergência, Consistência e Estabilidade
Este resultado pode ainda ser rearrumado de modo a obtermos
φ(x, y,∆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 desenvolvimento
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 ordiná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.
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 incre-
mento φ(x, y,∆x) mediante o teorema
Capítulo 1. Equações Diferenciais Ordinárias 33Capítulo 1. Equações Diferenciais Ordinárias 33
Teorema 3 (Asaithambi [3]) Um método numérico do tipo one-step é estável se ele for re-
gular.
Em seguida, veremos em que condições o método numérico pode ser considerado regular.
Da forma geral dos métodos one-step podemos dizer que ele será regular se a função incre-
mento 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 eM é 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 considerar-
mos, 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|φ(x, y,∆x)− φ(x, z,∆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|
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
34 1.4. Convergência, Consistência e Estabilidade34 1.4. Convergência, Consistência e Estabilidade
Teorema 4 (Asaithambi [3]) Um método numérico one-step regular é convergente se e so-
mente 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.
Como nós sabemos que a solução exata menos a solução numérica é igual ao erro de trun-
camento, 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, yi),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!
]
Capítulo 1. Equações Diferenciais Ordinárias 35Capítulo 1. Equações Diferenciais Ordinárias 35
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 mul-
tistep 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étodosmultistep. 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 trivial
dy
dx
= 0
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
36 1.4. Convergência, Consistência e Estabilidade36 1.4. Convergência, Consistência e Estabilidade
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 estabilidade
dos métodos multistep serão estabelecidas em função dessas 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 raiz 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 pro-
blema de valor inicial, para λ1 ∼= c, como [3]
yi = c+
k∑
j=2
γjλ
i
j
onde o somatório pode ser visto como os erros de arredondamento associados à solução apro-
ximada. Deste resultado vemos que o erro de arredondamento irá crescer exponencialmente 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 raiz:
|λj| ≤ 1, j = 1, 2, . . . , k
|λj| = 1 ⇒ dp(λj)
dλ
6= 0
Capítulo 1. Equações Diferenciais Ordinárias 37Capítulo 1. Equações Diferenciais Ordinárias 37
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 raiz e
λ = 1 é a única raiz característica de módulo igual a um. Ele é dito ser fracamente estável se
satisfaz a condição da raiz e tem mais de uma raiz 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 raiz.
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 a condição da raiz.
38 1.4. Convergência, Consistência e Estabilidade38 1.4. Convergência, Consistência e Estabilidade
CAPÍTULO 2
EQUAÇÕES DIFERENCIAIS PARCIAIS
Como veremos mais adiante, as equações que governam o escoamento de fluidos e o trans-
porte de massa e energia são equações diferenciais parciais, contendo termos em derivadas
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, pa-
rabólicas e hiperbólicas. Em seguida elas serão examinadas segundo as suas características
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 transporte do 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 importan-
tes 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.
39
40 2.1. Algumas Equações Diferenciais de Interesse40 2.1. Algumas Equações Diferenciais de Interesse
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.
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 dissipaçã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 prescritas
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 advecçã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.
Capítulo 2. Equações Diferenciais Parciais 41Capítulo 2. Equações Diferenciais Parciais 41
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 Rey-
nolds, governado pela equação de Stokes.
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 distribuição
de capacitância, indutância e resistência. As aplicações desta equação também incluem
o movimento