Embed Size (px)
Citation preview
Universidade Federal do Paraná � UFPR
Programa de Pós-Graduação em Engenharia Química� PPGEQ
Métodos Numéricos em EngenhariaQuímica
Prof. Éliton Fontana
2018/1
Conteúdo
1. Introdução 3
1.1. Classi�cação das Equações Diferenciais . . . . . . . . . . . . . . . . . . . . . 4
1.1.1. Equações Diferenciais Ordinárias (EDO) e Parciais (EDP) . . . . . . 4
1.1.2. Ordem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.3. Linearidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1.4. Homogeneidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1.5. Coe�cientes constantes e variáveis . . . . . . . . . . . . . . . . . . . . 6
I. Métodos Numéricos para Análise de Equações Diferenciais 7
2. Métodos Numéricos Para PVI's 8
2.1. Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.1.1. Erro Associado ao Método de Euler . . . . . . . . . . . . . . . . . . . 11
2.1.2. Método de Euler para Sistemas de EDO's . . . . . . . . . . . . . . . 13
2.2. Derivações do Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.1. Método de Euler Aprimorado (Fórmula de Heun) . . . . . . . . . . . 16
2.2.2. Método de Euler Implícito (Inverso) . . . . . . . . . . . . . . . . . . . 18
2.2.3. Regra do Ponto Médio . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3. Métodos de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1. Método de Runge-Kutta de 2a Ordem . . . . . . . . . . . . . . . . . 22
2.3.2. Método de Runge-Kutta de 4a Ordem . . . . . . . . . . . . . . . . . 24
2.3.3. Método de Runge-Kutta para Sistemas . . . . . . . . . . . . . . . . . 25
2.3.4. Métodos de Runge-Kutta Adaptativos . . . . . . . . . . . . . . . . . 28
2.4. Métodos de Passos Múltiplos . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.4.1. Método de Adams-Bashforth . . . . . . . . . . . . . . . . . . . . . . . 33
2.4.2. Método de Adams-Moulton . . . . . . . . . . . . . . . . . . . . . . . 34
2
2.5. Estabilidade dos Métodos de Resolução de EDO's . . . . . . . . . . . . . . . 37
2.6. Análise de Estabilidade de Métodos de Passo Único . . . . . . . . . . . . . . 38
2.6.1. Equação Modelo: Decaimento de Primeira Ordem . . . . . . . . . . . 39
2.7. Problemas Rígidos (Sti�) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.7.1. Sistemas de Equações Diferenciais Rígidas Lineares . . . . . . . . . . 43
2.7.2. Sistemas Não-Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3. Métodos Numéricos para Problemas de Valor de Contorno 46
3.1. Estratégias de Solução de PVC's . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2. Aproximações por Diferenças Finitas . . . . . . . . . . . . . . . . . . . . . . 48
3.2.1. Aproximação da Derivada Primeira . . . . . . . . . . . . . . . . . . . 49
3.2.2. Aproximação da Derivada Segunda . . . . . . . . . . . . . . . . . . . 52
3.2.3. Discretização das Condições de Contorno . . . . . . . . . . . . . . . . 53
3.2.4. Discretização de PVC's utilizando Diferenças Finitas . . . . . . . . . 54
4. Método das Linhas 60
3
1. Introdução
A maioria dos sistemas de interesse na engenharia envolvem mais de uma variável de-
pendente que são função da mesma variável independente, como por exemplo a variação na
concentração de diversas espécies químicas em um reator ao longo do tempo ou a variação
nas três componentes do vetor velocidade ao longo de uma dada direção espacial. De forma
geral, não é possível obter a solução para uma única variável independentemente, pois usu-
almente existe uma dependência de uma sobre a outra. Isto faz com que os modelos gerem
sistemas de equações que devem ser resolvidas ao mesmo tempo, podendo estas equações
serem algébricas, diferencias, integrais ou uma mistura delas.
Além disso, uma equação diferencia de ordem n pode sempre ser transforma em um sis-
tema de n equações diferenciais de primeira ordem. Esta abordagem é particularmente útil
para obter soluções numéricas de Problemas de Valor Inicial (PVI's), já que os métodos de
resolução (Euler, Runge-Kutta, etc.) são baseados na aproximação da derivada primeira e
portanto só podem ser utilizados para resolver equações de primeira ordem.
No momento, o objetivo será a análise de sistemas de equações diferenciais ordinárias
(EDO's), ou seja, equações que possuem somente uma variável independente. Esta análise
pode ser conduzida de forma quantitativa ou qualitativa, sendo que cada forma possui suas
vantagens e desvantagens. A análise quantitativa se refere à resolução direta das equações,
seja por via analítica ou numérica, de forma a se obter como as variáveis dependentes va-
riam em função da independente. Esta forma de análise é útil quando se deseja conhecer esta
dependência para uma condição especí�ca, porém trás pouca informação sobre o comporta-
mento global do sistema e como uma variável in�uencia as demais. A análise qualitativa,
por sua vez, consiste em examinar o comportamento geral de todas as famílias de solução
do sistema, o que permite obter informações valiosas sobre a estabilidade das soluções.
O material a seguir será dividido em duas partes. Na primeira parte, serão apresentados
métodos numéricos que podem ser aplicados na análise de equações diferenciais ordinárias,
além do Método das Linhas, que pode ser utilizado para transformar equações diferenciais
4
parciais em um sistema de EDO's. Na segunda parte serão apresentados métodos qualitativos
de análise e uma introdução à análise de estabilidade de equações não-lineares. Ao longo
do texto serão feitas referências às principais características das equações diferenciais, pois a
forma de análise pode depender de determinadas classi�cações. Para facilitar, a seguir será
apresentada uma breve revisão da classi�cação das equações diferenciais .
1.1. Classi�cação das Equações Diferenciais
As equações diferenciais são classi�cadas de acordo com suas características. Esta classi-
�cação é essencial para a determinação dos métodos de análise mais adequados para cada
equação.
1.1.1. Equações Diferenciais Ordinárias (EDO) e Parciais (EDP)
Quando a função desconhecida depende somente de uma variável independente a equação
é chamada de Equação Diferencial Ordinária (EDO). Por exemplo:
4d2y
dt2+ sin(t)
dy
dt= et
De forma genérica, uma EDO pode ser expressa como:
f(t, y, y′, y′′, . . . , y(n)) = 0
Caso a variável dependente dependa de duas ou mais variáveis, a equação é classi�cada
como Equação Diferencial Parcial (EDP). Por exemplo:
∂T
∂t= α2∂
2T
∂x2
1.1.2. Ordem
A ordem de uma equação diferencial é a ordem da derivada de maior ordem que aparece
na equação. Por exemplo:d2y
dt2= 0 segunda ordem
dy
dt+ y3 = 0 primeira ordem
A ordem da equação não está relacionada com a maior potência na qual a variável está
elevada.
5
1.1.3. Linearidade
Uma equação algébrica é dita linear quando pode ser escrita da forma:
a1x1 + a2x2 + a3x3 + . . . anxn + b = 0
ou seja, quando as variáveis x1, x2, . . . xn não aparecem em termos não-lineares.
De forma semelhante, a equação diferencial ordinária
F (t, y, y′, y′′, . . . , y(n)) = 0
é dita linear se F é uma função linear em relação à variável dependente e suas derivadas
(y, y′, y′′, . . . , y(n)).
Isto implica que uma EDO linear pode ser escrita da forma:
an(t)dny
dtn+ an−1(t)
d(n−1)y
dt(n−1)+ . . .+ a0(t)y = g(t)
ou seja, a variável y ou suas derivadas não estão presentes em termos não-lineares, como
por exemplo produto entre a variável e as derivadas, yn com n 6= 1 e funções não lineares
contendo a variável y (funções trigonométricas, exponencial, etc.). A variável independente
(t) pode estar presente em termos não-lineares. Exemplos de equações lineares:
d2y
dt2+ t
dy
dt+ y = 0
dy
dt+ t2y = sin(t)
dy
dt− yet = 0
Quando as equações não podem ser expressas da forma apresentada anteriormente são
consideradas não-lineares. As equações diferenciais não-lineares possuem em geral uma re-
solução muito mais complexa e por isso são mais difíceis de serem avaliadas analiticamente.
Exemplos de equações não-lineares:
dy
dt+ y2 = 0 y
dy
dt= t
d2y
dt2= sin(y)
dy
dt+ tey = 0
De forma similar, uma EDP é dita não-linear quando possui algum termo não linear
envolvendo qualquer variável dependente.
As EDO's lineares ainda podem ser divididas em outras categorias que são utilizadas para
facilitar a escolha dos métodos de solução:
1.1.4. Homogeneidade
Uma equação diferencial do tipo
F (t, y, y′, y′′, . . . , y(n)) = 0
6
é homogênea quando o termo associado somente à variável independente é nulo. Na forma
apresentada anteriormente para as EDO's lineares, a equação é homogênea quando g(t) = 0.
Caso algum termo diferente de zero não estiver multiplicado pela variável dependente (y) ou
alguma de suas derivadas, a equação é dita não-homogênea;
1.1.5. Coe�cientes constantes e variáveis
Esta classi�cação costuma ser empregada somente para EDO's. Na expressão para uma
EDO genérica vista anteriormente, os coe�cientes an(t), an−1(t)...a0(t) podem ser funções
conhecidas da variável independente t. No entanto, caso estes coe�cientes sejam constantes
(não dependam de t) a equação é conhecida como EDO com coe�cientes constantes, caso
contrário é chamada de EDO com coe�cientes variáveis.
Uma classe de equações muito importante nas ciências exatas são as equações autônomas.
Estas equações surgem quando a variável independente não aparece de forma explícita na
equação.
7
Parte I.
Métodos Numéricos para Análise de
Equações Diferenciais
8
2. Métodos Numéricos Para PVI's
2.1. Método de Euler
A maior parte das equações diferenciais de interesse na engenharia não possuem solução
analítica conhecida. Além disso, em muitos casos as soluções analíticas co-nhecidas são
difíceis de serem utilizadas, por exemplo quando envolvem séries in�nitas ou integrais sem
resolução analítica. Nestes casos é mais conveniente a utilização de ferramentas numéricas
para a resolução das equações.
Existem basicamente três formas de analisar equações diferenciais: através de métodos
analíticos, métodos qualitativos e métodos numéricos. Os métodos analíticos permite esta-
belecer uma relação direta entre as variáveis dependentes e independentes através de funções
válidas em um determinado intervalo. Os métodos qualitativos envolvem a análise do com-
portamento geral da equação sem necessariamente buscar uma solução, como por exemplo
através da construção do campo de direções.
Os métodos de resolução numérica permitem que a solução seja estimada para uma dada
condição inicial e um conjunto de parâmetros especí�cos. A ideia geral é de que é possível
obter aproximações da solução de uma determinada EDO avançando em pontos que estão
a uma distância ∆t da condição inicial. Diferentemente dos métodos analíticos, a solução
obtida através de métodos numéricos é válida somente para um conjunto de parâmetros
especí�cos e para uma dada condição inicial. Caso esta condição inicial ou algum parâmetro
sejam alterados, deve-se resolver o problema novamente.
A determinação de qual forma de análise deve ser utilizada irá depender tanto das caracte-
rísticas das equações diferenciais que estão sendo avaliadas quanto de que tipo de informação
está se buscando.
Assim como para o caso da resolução de sistemas de equações algébricas, existem diversos
métodos que podem ser utilizados para a resolução de sistemas de equações diferenciais.
Inicialmente será apresentado o método de Euler. Este método costuma ser pouco e�ciente
9
para a resolução de sistemas complexos, porém diversos conceitos fundamentais dos métodos
de resolução podem ser mais facilmente apresentados através deste método.
O método de Euler é o método mais simples para a aproximação da solução de EDO's.
Considere o seguinte problema de valor inicial de primeira ordem:
y′ =dy
dt= f(t, y) y(t0) = y0
A solução deste PVI passa obrigatoriamente pelo ponto (t0, y0). O objetivo do método de
Euler é estimar o valor da solução y(t) em um ponto com uma distância ∆t do ponto inicial
t0. Este ponto passa a ser chamado de t1 = t0 + ∆t, de modo que o do método de Euler é
determina o valor da solução em t1, y1 = y(t1). A partir deste valor, pode-se estimar o valor
da função em um ponto t2 distante ∆t do ponto t1 e continuar com este procedimento até
se obter uma aproximação para a solução por quantos pontos forem necessários. Esta classe
de métodos são chamados de métodos passo-a-passo. Para avaliar o termo y1 = y(t0 + ∆t),
pode-se aplicar uma expansão em séries de Taylor em torno do ponto t0:
y(t1) = y(t0 + ∆t) = y0 + ∆ty′(t0) +∆t2
2!y′′(t0) +
∆t3
3!y′′′(t0) + . . .
Considerando que o espaçamento ∆t seja muito pequeno, os termos ∆t2, ∆t3, etc. podem
ser desprezados. Com isso:
y(t0 + ∆t) = y1 ≈ y0 + ∆ty′(t0)
Considerando também que y′ = f(t, y), a relação anterior pode ser reescrita como:
y1 = y0 + ∆tf(t0, y0)
De forma geral, a relação anterior pode ser utilizada para estimar a solução para qualquer
valor yn+1 com base no valor para yn:
yn+1 = yn + ∆tf(tn, yn)
Como o valor no novo ponto yn+1 é calculado com base somente em valores no ponto
anterior (já conhecidos), este método é dito explícito.
Estrutura de um código para resolver uma EDO dy/dt = f(t, y) utilizando o método de
Euler:
Passo 1: De�nir f(t, y);
Passo 2: De�nir os valores iniciais t[1] = t0, y[1] = y0;
10
Passo 3: De�nir o tamanho do passo ∆t e o número de passos n;
Passo 4: Para i = 1 até i = n calcular:
k1 = f(t[i], y[i])
y[i+ 1] = y[i] + ∆t× k1
t[i+ 1] = t[i] + ∆t
(2.1)
Passo 5: Plotar os resultados.
Exemplo 01): Utilizando o método de Euler, obtenha a solução aproximada até t = 0.5
para o seguinte PVI, utilizando ∆t = 0.1:
dy
dt= −2t− y y(0) = −1
Neste caso f(t, y) = −2t− y. No ponto inicial temos:
t0 = 0 y0 = −1
Como o passo ∆t é constante, os valores de tn são automaticamente de�nidos, t1 = t0 +
∆t = 0.1, t2 = t1 + ∆t = 0.2, . . .
Avaliando y1 = y(t1) = y(0.1):
y1 = y0 + ∆tf(t0, y0) = −1 + 0.1(−2× 0− (−1)) = −0.9
Repetindo para os próximos passos:
y2 = y1 + ∆tf(t1, y1) = −0.9 + 0.1(−2× 0.1− (−0.9)) = −0.83
y3 = y2 + ∆tf(t2, y2) = −0.83 + 0.1(−2× 0.2− (−0.83)) = −0.787
y4 = y3 + ∆tf(t3, y3) = −0.787 + 0.1(−2× 0.3− (−0.787)) = −0.7683
y5 = y4 + ∆tf(t4, y4) = −0.7683 + 0.1(−2× 0.4− (−0.7683)) = −0.77147
Para comparação, a solução exata da equação diferencial é:
y(t) = −3e−t − 2t+ 2
sendo que para os pontos avaliados a solução exata é:
y0 = −1 y1 = −0.9145 y2 = −0.85619 y3 = −0.8224
y4 = −0.8109 y5 = −0.81959
11
Exemplo 02): Utilizando o método de Euler, obtenha a solução aproximada para t = 1
para o seguinte PVI, utilizando ∆t = 0.2:
dy
dt= 2t2y2 y(0) = 1
Com base nos dados fornecidos, t0 = 0 e y0 = 1. Como ∆t = 0.2, temos que t1 = 0.2,
t2 = 0.4, t3 = 0.6, t4 = 0.8 e t5 = 1. Avaliando os termos intermediários:
y1 = y0 + 0.2 · (2 · t20y20) = 1 + 0.2 · (2 · 02 · 12) = 1
y2 = y1 + 0.2 · (2 · t21y21) = 1 + 0.2 · (2 · 0.22 · 12) = 1.016
y3 = y2 + 0.2 · (2 · t22y22) = 1 + 0.2 · (2 · 0.42 · 1.0162) = 1.0821
y4 = y3 + 0.2 · (2 · t23y23) = 1 + 0.2 · (2 · 0.62 · 1.08212) = 1.2507
y5 = y4 + 0.2 · (2 · t24y24) = 1 + 0.2 · (2 · 0.82 · 1.25072) = 1.6511
A solução exata em t = 1 é y(1) = 3, portanto a solução obtida está com um erro
signi�cativo. Por este exemplo, �ca claro que é fundamental avaliar os erros associados à
utilização do método de Euler.
2.1.1. Erro Associado ao Método de Euler
A utilização de métodos numéricos sempre leva à obtenção de soluções aproximadas. Uma
importante propriedade dos métodos numéricos é a convergência. Um método é dito
convergente quando a solução obtida tende à solução exata quando o espaçamento ∆t→ 0.
Na resolução de EDO's pelo método de Euler (ou por qualquer outro método), existem
duas fontes de erros:
- Erros de truncamento: erros associados com o truncamento da expansão em série de
Taylor no segundo termo. Como a determinação do ponto yn+1 depende do valor de yn, o
erro de truncamento pode aumentar rapidamente ao longo da resolução. Para reduzir os
erros de truncamento, pode-se reduzir o passo de tempo ∆t.
- Erros de arredondamento: erro associado à precisão dos valores numéricos utilizados
(número de dígitos signi�cativos: 6-9 para precisão única e 15-17 para precisão dupla).
Quanto maior for o número de operações necessárias, maior será o erro de arredondamento.
Assim, uma maneira de reduzir este erro é aumentar o passo de tempo, de modo que menos
passos precisam ser resolvidos para atingir o tempo desejado.
12
Dessa forma, conforme o passo de tempo aumenta, o erro de arredondamento diminui e o
de truncamento aumenta. Por isso, existe um ponto ótimo onde o erro total é mínimo, como
ilustrado na �gura a seguir. Usualmente os erros de arredondamento são menos signi�cativos,
por isso a tendência é que passos pequenos gerem melhores resultados.
Como visto anteriormente, pode-se obter o valor de uma função y(t) em um ponto tn+1
com base no valor em tn através de uma expansão em séries de Taylor:
yn+1 = yn + ∆t y′(tn) +∆t2
2!y′′(tn) +
∆t2
3!y′′′(tn) + . . .
A diferença entre o valor exato de yn+1 e o valor aproximado obtido com o uso do método
de Euler corresponde ao erro local de truncamento (para um determinado valor de n):
en+1 =∆t2
2!y′′(tn) +
∆t2
3!y′′′(tn) + . . .
Nesta forma, o erro ainda representa uma série in�nita. Para evitar este problema, pode-se
utilizar a fórmula de Taylor com resto de Lagrange:
yn+1 = yn + ∆t y′(tn) +∆t2
2!y′′(c)
onde c ∈ [tn, tn+1]. Este resultado é uma extensão do teorema do valor médio, que diz que
dada uma função contínua f de�nida num intervalo fechado [a,b] e diferenciável em (a,b),
existe algum ponto c em (a,b) tal que :
f ′(c) =f(b)− f(a)
b− a
Utilizando a relação anterior, o erro en+1 pode ser avaliado como:
en+1 =y′′(c)
2∆t2
13
Como y′′(c) é uma constante (valor �nito), temos que:
en+1 = O(∆t2)
ou seja, o erro local é da ordem de ∆t2.
Como cada passo gera um erro local, quanto mais passos forem utilizados, maior será o
acúmulo de erros. Por isso, quando o valor yn+1 estiver sendo determinado, o erro associado
à determinação do valor de yn também deve ser considerado. Para n passos, o erro global
En pode ser avaliado como:
En+1 = en+1 · n = O(h2) · n = O(h2)n · hh
= O(h)tn
onde h = ∆t. Novamente, como tn é um escalar:
En+1 = O(h)
Assim, o erro global no método de Euler é da ordem do tamanho do passo utilizado. Esta
relação mostra que o erro tende a zero conforme h→ 0 e portanto o método é convergente.
De forma geral, um dado método é convergente quando:
En+1 = O(hp) p > 0
Neste caso, o método possui ordem p, portanto, o método de Euler é um método de primeira
ordem. Métodos de ordem superior convergem mais rapidamente para a solução exata (não
exigem passos tão pequenos), sendo por isso mais indicados que os métodos de primeira
ordem.
A utilização de um passo muito pequeno (tendendo a zero) é impraticável, pois exigiria
um tempo computacional demasiadamente alto, além de fazer com que os erros de arredon-
damento aumentassem rapidamente. Um procedimento simples que pode na maioria dos
casos ser utilizado para determinar se o passo utilizado é adequado é resolver o problema
com valores de ∆t gradativamente menores. A partir de um determinado ponto as soluções
obtidas serão muito parecidas, sendo que uma maior redução em ∆t a partir deste ponto não
irá reduzir o erro global de forma signi�cativa e irá aumentar o erro de arredondamento.
2.1.2. Método de Euler para Sistemas de EDO's
O método de Euler explícito pode ser estendido para aplicação em sistemas de EDO's de
primeira ordem. Considere o seguinte sistema de PVI's:
dx
dt= f(t, x, y) x(0) = x0
14
dy
dt= g(t, x, y) y(0) = y0
De forma semelhante ao realizado para uma única equação, pode-se partir das condições
iniciais e ir avançando ao longo de t com base nos valores já conhecidos:
xn+1 = xn + ∆t · f(tn, xn, yn)
yn+1 = yn + ∆t · g(tn, xn, yn)
Generalizando para um sistema de m equações da forma:
dy1
dt= f1(t, y1, y2, y3, . . . , ym)
dy2
dt= f2(t, y1, y2, y3, . . . , ym)
dy3
dt= f3(t, y1, y2, y3, . . . , ym)
...dym
dt= fm(t, y1, y2, y3, . . . , ym)
y1(t0) = y10
y2(t0) = y20
y3(t0) = y30
...
ym(t0) = ym0
(2.2)
A utilização do método de Euler explícito implica em:
y1n+1
y2n+1
y3n+1
...
ymn+1
=
y1n
y2n
y3n
...
ymn
+ ∆t
f1(tn, y1n, y
2n, . . . , y
mn )
f2(tn, y1n, y
2n, . . . , y
mn )
f3(tn, y1n, y
2n, . . . , y
mn )
...
fm(tn, y1n, y
2n, . . . , y
mn )
(2.3)
A aplicação do método para resolução de sistemas também permite que ele seja utilizado
para a resolução de equações de ordem maior que 1. Por exemplo, considere o seguinte PVI:
d2y
dx2= f
(x, y,
dy
dx
)y(0) = y0
dy
dx
∣∣∣x=0
= y1
Esta equação pode ser transformada em um sistema de primeira ordem através da de�nição
de uma variável adicional u:
u =dy
dx→ d2y
dx2=du
dx
Com isso, o PVI pode ser reescrito como:
dy
dx= u y(0) = y0
15
du
dx= f (x, y, u) u(0) = y1
Este procedimento pode ser aplicado para escrever uma EDO de qualquer ordem como
um sistema de EDO's de primeira ordem.
Exemplo 03: Utilizando o método de Euler, obtenha a solução aproximada para os três
primeiros passos de tempo para o seguinte sistema, considerando ∆t = 0.1:
dx
dt= y x(0) = 0
dy
dt= x√y + 3y y(0) = 2
Desse modo, f(t, x, y) = y e g(t, x, y) = x√y + 3y.
Neste caso, temos que t0 = 0, x0 = 0 e y0 = 2. Para avaliar os próximos passos, deve-se
determinar todas as variáveis em t1, na sequência todas as variáveis em t2 e assim sucessi-
vamente:
x1 = x0 + ∆t · f(t0, x0, y0) = 0 + 0.1 · (2) = 0.2
y1 = y0 + ∆t · g(t0, x0, y0) = 2 + 0.1 · (0 ·√
2 + 3 · 2) = 2.6
x2 = x1 + ∆t · f(t1, x1, y1) = 0.2 + 0.1 · (2.6) = 0.46
y2 = y1 + ∆t · g(t1, x1, y1) = 2.6 + 0.1 · (0.2 ·√
2.6 + 3 · 2.6) = 3.4122
x3 = x2 + ∆t · f(t2, x2, y2) = 0.46 + 0.1 · (3.4122) = 0.80122
y3 = y2 + ∆t · g(t2, x2, y2) = 3.4122 + 0.1 · (0.46 ·√
3.4122 + 3 · 3.4122) = 4.520896
16
2.2. Derivações do Método de Euler
A utilização do método de Euler para a resolução de problemas típicos da engenharia
usualmente implica no uso de passos de tempo muito pequenos para garantir a convergência.
Normalmente, melhores resultados são obtidos com o uso de métodos de maior ordem ou
que ao menos possuam características de convergência mais atrativas.
Como visto na aula passada, o método de Euler possui diversas limitações que restringem
o seu uso a problemas muito especí�cos. Dentre as principais causas destas limitações pode-
se citar os termos desprezados na expansão em série de Taylor, a necessidade de utilizar
passos de tempo muito pequenos e o fato de que o método considera que o valor da derivada
y′(ti) é constante em todo o intervalo (ti, ti+1). O método de Euler serviu como base para o
desenvolvimento de métodos gradativamente mais complexos e com melhor desempenho.
O método de Euler e outros métodos explícitos de maior ordem fazem parte de uma
família de métodos numéricos mais abrangentes chamados de métodos de Runge-Kutta (RK).
O método de Euler corresponde ao método RK de primeira ordem. Além dos métodos
RK, outras abordagens podem ser empregadas, como por exemplo a utilização de métodos
implícitos e os métodos de passos múltiplos. A seguir serão apresentadas três modi�cações
do método de Euler que pertencem a estas categorias.
Assim como para o método de Euler, o objetivo continua sendo buscar uma solução apro-
ximada para o PVI:dy
dt= f(t, y) y(t0) = y0
2.2.1. Método de Euler Aprimorado (Fórmula de Heun)
Este método faz parte da categoria de passos múltiplos, onde os valores da derivada em
mais de um ponto são utilizados para determinar o próximo ponto.
A fórmula de Heun é uma modi�cação simples do método de Euler e consiste em uma
abordagem preditiva-corretiva. O passo inicial consiste determinar um valor preditor com
base no método de Euler:
y∗n+1 = yn + ∆t · f(tn, yn)
A partir deste valor, é calculada uma etapa corretora com base no valor médio da função
avaliada em yn e no valor preditor y∗n+1:
yn+1 = yn +∆t
2· (f(tn, yn) + f(tn+1, y
∗n+1))
17
Esta fórmula representa uma correção do valor estimado anteriormente y∗n+1.
Para implementar computacionalmente este método, pode-se utilizar a mesma estrutura
geral apresentada para o método de Euler, sendo necessário somente modi�car o passo 4:
Passo 4: Para i = 1 até i = n calcular:
k1 = f(t[i], y[i])
yp[i+ 1] = y[i] + ∆t× k1
t[i+ 1] = t[i] + ∆t
k2 = f(t[i+ 1], yp[i+ 1])
y[i+ 1] = y[i] + (∆t/2)(k1 + k2)
Exemplo 01) Repita o exemplo visto anteriormente utilizando a fórmula de Euler apri-
morada.
Lembrando do PVI (∆t = 0.2):
dy
dt= 2t2y2 y(0) = 1
A resolução com o método de Euler modi�cado envolve mais etapas:
- Passo 1:
k1 = f(t0, y0) = 0 yp1 = y0 + ∆t · k1 = 1 t1 = t0 + ∆t = 0.2
k2 = f(t1, yp1) = 0.08 y1 = y0 + (∆t/2) · (k1 + k2) = 1.008
- Passo 2:
k1 = f(t1, y1) = 0.08128 yp2 = y1 + ∆t · k1 = 1.0243 t2 = t1 + ∆t = 0.4
k2 = f(t2, yp2) = 0.3357 y2 = y1 + (∆t/2) · (k1 + k2) = 1.0497
- Passo 3:
k1 = f(t2, y2) = 0.35259 yp3 = y2 + ∆t · k1 = 1.1202 t3 = t2 + ∆t = 0.6
k2 = f(t3, yp3) = 0.90349 y3 = y2 + (∆t/2) · (k1 + k2) = 1.1753
- Passo 4:
k1 = f(t3, y3) = 0.99456 yp4 = y3 + ∆t · k1 = 1.3742 t4 = t3 + ∆t = 0.8
18
k2 = f(t4, yp4) = 2.4172 y4 = y3 + (∆t/2) · (k1 + k2) = 1.5165
- Passo 5:
k1 = f(t4, y4) = 2.94371 yp5 = y4 + ∆t · k1 = 2.1052 t5 = t4 + ∆t = 1
k2 = f(t5, yp5) = 8.8637 y5 = y4 + (∆t/2) · (k1 + k2) = 2.6973
Para comparação, a solução exata para t = 0.8 é y = 1.5182 e para t = 1 é y = 3.
2.2.2. Método de Euler Implícito (Inverso)
O método de Euler explícito utiliza a inclinação da reta tangente no ponto n para estimar
o valor da solução no ponto n+ 1:
yn+1 = yn + ∆t · (f(tn, yn))
De forma alternativa, poderia-se utilizara a inclinação da reta tangente no próprio ponto
n+ 1. Esta formulação dá origem ao método de Eulet implícito ou inverso ou regressivo:
yn+1 = yn + ∆t · (f(tn+1, yn+1))
Como a variável yn+1 aparece nos dois lados da igualdade, pode ser necessário utilizar
algum método iterativo para a resolução do problema, caso não seja possível explicitar o
valor yn+1. Para EDO's lineares sempre é possível isolar yn+1, no entanto para não-lineares
isso usualmente não é possível.
Esta formulação representa um método implícito, pois a variável yn+1 depende impli-
citamente dela mesmo. Assim como a formulação explícita, este método representa uma
aproximação de primeira ordem, porém, como será visto nas próximas aulas, os métodos
implícitos possuem uma grande vantagem sobre os explícitos em relação à estabilidade do
método.
Exemplo 02) Resolva o seguinte PVI utilizando o método de Euler implícito, com ∆t =
0.1 até t = 0.3:dy
dt= t2 + ty y(0) = 1
19
2.2.3. Regra do Ponto Médio
Lembrando da Fórmula de Taylor com o resto de Lagrange:
yn+1 = yn + ∆t y′(tn) +∆t2
2!y′′(c)
Quando os termos de segunda ordem são desconsiderados, obtém-se o método de Euler
explícito. Uma maneira de se conseguir uma aproximação mais precisa (com maior ordem) é
considerar um maior número de termos na expansão. Por exemplo, caso o termo de terceira
ordem também for considerado, temos que:
yn+1 = yn + ∆t y′(tn) +∆t2
2!y′′(tn) +
∆t3
3!y′′′(c)
onde novamente c ∈ (tn, tn+1).
De forma semelhante, pode-se substituir ∆t por −∆t para se obter uma aproximação para
yn−1:
yn−1 = yn −∆t y′(tn) +∆t2
2!y′′(tn)− ∆t3
3!y′′′(d)
onde d ∈ (tn−1, tn).
Fazendo a expansão para yn+1 menos a expansão para yn−1:
yn+1 − yn−1 = 2∆ty′(tn) +y′′′(c) + y′′′(d)
3!∆t3
Novamente, como as derivadas y′′′(c) e y′′′(d) são dois escalares, a relação anterior pode
ser expressa como:
yn+1 = yn−1 + 2∆ty′(tn) +O(∆t3)
Desconsiderando os termos de terceira ordem:
yn+1 = yn−1 + 2∆ty′(tn)
Esta relação é conhecida como regra do ponto médio e representa uma aproximação de
segunda ordem para yn+1. Esta fórmula depende do conhecimento do valor da variável em
três pontos, por isso não pode ser utilizara para determinar y1, já que y−1 não existe. Outra
relação pode ser utilizada para determinar o ponto y1, e a relação anterior pode ser usada
para os demais pontos.
Comparando com o método de Euler explícito:
yn+1 = yn + ∆ty′(tn) → y′(tn) =yn+1 − yn
∆t
20
yn+1 = yn−1 + 2∆ty′(tn) → y′(tn) =yn+1 − yn−1
2∆t
Exemplo 07:) Utilize a regra do ponto médio para resolver o seguinte PVI até t = 0.3,
utilizando ∆t = 0.1:dy
dt= y + 2t+ t2 y(0) = 1
Para avaliar o primeiro ponto, pode-se utilizar alguma formulação distinta. Como a regra
do ponto médio é um método de segunda ordem, é recomendável a utilização de outro método
de segunda ordem para não haver uma perda muito grande de precisão.
Utilizando o método de Euler aprimorado:
k1 = f(t0, y0) = 1
yp1 = y0 + ∆t · k1 = 1.1
t1 = t0 + ∆t = 0.1
k2 = f(t1, y1) = 1.31
y1 = y0 + (∆t/2)(k1 + k2) = 1.1155
A partir dos valores conhecidos para y0 e y1, pode-se agora determinar os demais pontos:
y2 = y0 + 2∆tf(t1, y1) = 1 + 2 · 0.1 · (1.1155 + 2 · 0.1 + 0.12) = 1.2651
y3 = y1 + 2∆tf(t2, y2) = 1.1155 + 2 · 0.1 · (1.2651 + 2 · 0.2 + 0.22) = 1.4565
21
2.3. Métodos de Runge-Kutta
Os métodos de Runge-Kutta (RK) utilizam a aproximação em série de Taylor para obter
formulações de alta ordem sem necessitar o cálculo das derivadas de alta ordem. De forma
geral, estes métodos podem ser expressos de forma geral como:
yi+1 = yi + ∆tφ(ti, yi,∆t)
onde a função φ(ti, yi,∆t) é chamada de incremento e pode ser interpretada como uma
representação da inclinação no intervalo entre ti e ti+1. Como esta função só depende de
pontos já conhecidos, os métodos de Runge-Kutta clássicos são explícitos.
Para um método de ordem n, a função incremento costuma ser expressa na seguinte forma:
φ(ti, yi,∆t) = a1k1 + a2k2 + a3k3 + . . .+ ankn
onde a1, a2, a3, . . . , an são constantes e os termos ki são dados por:
k1 = f(ti, yi)
k2 = f(ti + p1∆t, yi + q11k1∆t)
k3 = f(ti + p2∆t, yi + q21k1∆t+ q22k2∆t)
...
kn = f(ti + pn−1∆t, yi + qn−1,1k1∆t+ qn−1,2k2∆t+ . . .+ qn−1,n−1kn−1∆t)
onde os termos pi e qij são constantes. Como pode ser observado, a determinação do parâ-
metros kn depende os valores ki onde i < n.
Por exemplo, para um método de primeira ordem a função incremento é avaliada como:
φ(ti, yi,∆t) = a1k1
onde k1 = f(ti, yi). Substituindo na expressão para a forma geral dos métodos de Runge-
Kutta:
yi+1 = yi + ∆t(a1f(ti, yi))
Os métodos RK estão baseados diretamente na expansão em série de Taylor em torno do
ponto yi. Como visto anteriormente, a expansão de primeira ordem resulta em:
yi+1 = yi + ∆tf(ti, yi)
Comparando os termos semelhantes com a equação anterior, observa-se que a1 = 1. O
método de Runge-Kutta de primeira ordem corresponde, de fato, ao método de Euler.
22
2.3.1. Método de Runge-Kutta de 2a Ordem
Considere agora que se deseja obter uma aproximação de segunda ordem. Neste caso,
temos que:
yi+1 = yi + (a1k1 + a2k2)∆t
k1 = f(ti, yi)
k2 = f(ti + p1∆t, yi + q11k1∆t)
Para determinar os valores de a1, a2, p1 e q11 pode-se primeiramente avaliar a expansão
em séries de Taylor de segunda ordem em torno de yi é dada por:
yi+1 = yi + ∆tf(ti, yi) +df
dy
∆t2
2
Considerando que f = f(t, y), a derivada total de f em relação a t é dada por:
df
dt=∂f
∂t+∂f
∂y
dy
dt
Substituindo na expressão anterior:
yi+1 = yi + ∆tf(ti, yi) +
(∂f
∂t+∂f
∂y
dy
dt
)∆t2
2
Utilizando um procedimento semelhante, pode-se utilizar a expansão em séries de Taylor
para expandir a própria função f(t, y). Lembrando, para uma função de duas variáveis, a
expansão em série de Taylor é dada por:
f(x, y) = f(x0, y0) +∂f
∂x(x− x0) +
∂f
∂y(y − y0)+
1
2!
(∂2f
∂x2(x− x0)2 +
∂2f
∂x∂y(x− x0)(y − y0) +
∂2f
∂y2(y − y0)2
)+ . . .
De forma equivalente, pode-se avaliar a função em um ponto qualquer f(x+ a, y+ b) com
base no valor de f(x, y). Considerando uma aproximação de primeira ordem:
f(x+ a, y + b) = f(x, y) +∂f
∂xa+
∂f
∂yb
Assim, considerando a de�nição de k2, o termo f(ti+p1∆t, yi+q11k1∆t) pode ser expresso
como:
k2 = f(ti + p1∆t, yi + q11k1∆t) = f(ti, yi) + p1∆t∂f
∂t+ q11k1∆t
∂f
∂y
Substituindo os valores de k1 e k2 na expressão para yi+1 = yi+(a1k1+a2k2)∆t e colocando
os termos de mesma ordem em evidência:
yi+1 = yi + (a1f(ti, yi) + a2f(ti, yi))∆t+
(a2p1
∂f
∂t+ a2q11f(ti, yi)
∂f
∂y
)∆t2
23
Comparando com a equação obtida anteriormente:
yi+1 = yi + ∆tf(ti, yi) +
(∂f
∂t+∂f
∂y
dy
dt
)∆t2
2
Igualando os coe�cientes com os termos de mesma ordem, observa-se que as seguintes
relações devem ser satisfeitas:
a1 + a1 = 1
a2p1 = 1/2
a2q11 = 1/2
Este sistema possui três equações e quatro variáveis, portanto não existe solução única. No
entanto, a partir do momento em que um dos valores é especi�cado, os demais podem ser
calculados. Isto implica que existem in�nitas formulações para o método de Runge-Kutta
de segunda ordem.
Euler Modi�cado (a2 = 0.5)
Assumindo que a2 = 0.5, obtém-se que a1 = 0.5 e p1 = q11 = 1. Assim, o método de
Runge-Kutta de segunda ordem resulta em:
k1 = f(ti, yi)
k2 = f(ti + ∆t, yi + k1∆t)
yi+1 = yi +(k1 + k2)
2∆t
Este método corresponde ao método de Euler modi�cado (fórmula de Heun).
Regra do Ponto Médio (a2 = 1)
Fazendo a2 = 1, temos que a1 = 0 e p1 = q11 = 1/2:
k1 = f(ti, yi)
k2 = f(ti + ∆t/2, yi + k1∆t/2)
yi+1 = yi + k2∆t
Esta formulação corresponde ao método de Euler modi�cado com base na regra do ponto
médio.
24
Método de Ralston (a2 = 2/3)
Fazendo a2 = 2/3, os demais valores são avaliados como a1 = 1/3 e p1 = q11 = 3/4:
k1 = f(ti, yi)
k2 = f(ti + 3∆t/4, yi + 3k1∆t/4)
yi+1 = yi + (k1/3 + 2k2/3)∆t
Esta formulação é chamada de Método de Ralston
2.3.2. Método de Runge-Kutta de 4a Ordem
Seguindo o mesmo procedimento mostrado anteriormente para uma aproximação de quarta
ordem, obtém-se o método RK4. Este é possivelmente o método mais empregado para a
resolução de PVI's, devido a sua alta precisão e por ser um método explícito.
Neste caso, os parâmetros k1, k2, k3 e k4 são avaliados como:
k1 = f(ti, yi)
k2 = f(ti + ∆t/2, yi + k1∆t/2)
k3 = f(ti + ∆t/2, yi + k2∆t/2)
k4 = f(ti+1, yi + k3∆t)
Com base nestes valores, o ponto yn+1 é avaliado como:
yn+1 = yn + (1/6)(k1 + 2 · k2 + 2 · k3 + k4)
A utilização de passos de tempo ∆t constantes para todos os casos possui uma série de
desvantagens, especialmente porque em muitos casos a equação pode precisar de valores
muito pequenos para manter a precisão em alguns pontos, enquanto que em outros pontos
valores maiores poderiam ser utilizados para reduzir o número de cálculos realizados.
Muitos softwares utilizam uma abordagem que determina localmente o tamanho necessário
para os passos, através da comparação dos resultados obtidos com RK4 e uma formulação
de Runge-Kutta de quinta ordem (ex. ode45 do Matlab).
Exemplo 08:) Utilize o método RK4 para resolver o seguinte PVI, adotando ∆t = 0.1
até t = 0.3:dy
dt= −y + 2t y(0) = 2
25
2.3.3. Método de Runge-Kutta para Sistemas
O método de Runge-Kutta pode ser utilizado para a resolução de sistemas de EDO's.
Considere, por exemplo, um sistema com duas equações:
dy
dt= f(y, x, t) y(t0) = y0
dx
dt= g(y, x, t) x(t0) = x0
Neste caso os parâmetros k1, . . . , k4 irão depender de qual equação está sendo avaliada:
ky1 = f(tn, yn, xn)
kx1 = g(tn, yn, xn)
ky2 = f(tn + ∆t/2, yn + ky1∆t/2, xn + kx1 ∆t/2)
kx2 = g(tn + ∆t/2, yn + ky1∆t/2, xn + kx1 ∆t/2)
ky3 = f(tn + ∆t/2, yn + ky2∆t/2, xn + kx2 ∆t/2)
kx3 = g(tn + ∆t/2, yn + ky2∆t/2, xn + kx2 ∆t/2)
ky4 = f(tn + ∆t, yn + ky3∆t, xn + kx3 ∆t)
kx4 = g(tn + ∆t, yn + ky3∆t, xn + kx3 ∆t)
Com base nestes valores, a solução no ponto n+ 1 pode ser aproximada:
yn+1 = yn +∆t
6(ky1 + 2ky2 + 2ky3 + ky4)
xn+1 = xn +∆t
6(kx1 + 2kx2 + 2kx3 + kx4 )
Considere agora um sistema de m equações:
dy1
dt= f1(t, y1, y2, y3, . . . , ym)
dy2
dt= f2(t, y1, y2, y3, . . . , ym)
dy3
dt= f3(t, y1, y2, y3, . . . , ym)
...dym
dt= fm(t, y1, y2, y3, . . . , ym)
y1(t0) = y10
y2(t0) = y20
y3(t0) = y30
...
ym(t0) = ym0
(2.4)
Neste caso, as quantidades k1, k2, etc. passam a ser vetores com m componentes:
26
k1 =
f1(xn, y
1n, y
2n, . . . , y
mn )
f2(xn, y1n, y
2n, . . . , y
mn )
...
fm(xn, y1n, y
2n, . . . , y
mn )
k2 = ∆t
f1(xn + ∆t/2, y1
n + ∆t · k1[1]/2, y2n + ∆t · k1[2]/2, . . . ymn + ∆t · k1[m]/2)
f2(xn + ∆t/2, y1n + ∆t · k1[1]/2, y2
n + ∆t · k1[2]/2, . . . ymn + ∆t · k1[m]/2)...
fm(xn + ∆t/2, y1n + ∆t · k1[1]/2, y2
n + ∆t · k1[2]/2, . . . ymn + ∆t · k1[m]/2)
k3 = ∆t
f1(xn + ∆t/2, y1
n + ∆t · k2[1]/2, y2n + ∆t · k2[2]/2, . . . ymn + ∆t · k2[m]/2)
f2(xn + ∆t/2, y1n + ∆t · k2[1]/2, y2
n + ∆t · k2[2]/2, . . . ymn + ∆t · k2[m]/2)...
fm(xn + ∆t/2, y1n + ∆t · k2[1]/2, y2
n + ∆t · k2[2]/2, . . . ymn + ∆t · k2[m]/2)
k4 = ∆t
f1(xn + ∆t, y1
n + ∆t · k3[1], y2n + ∆t · k3[2], . . . ymn + ∆t · k3[m])
f2(xn + ∆t, y1n + ∆t · k3[1], y2
n + ∆t · k3[2], . . . ymn + ∆t · k3[m])...
fm(xn + ∆t, y1n + ∆t · k3[1], y2
n + ∆t · k3[2], . . . ymn + ∆t · k3[m])
Com base nestes vetores, a solução para o sistema pode ser avaliada:
yin+1 = yin +∆t
2(k1[i] + 2k2[i] + 2k3[i] + k4[i])
Exemplo 01:) Utilize o método RK4 para resolver os dois primeros passos para o seguinte
sistema, considerando ∆t = 0.2:
dy
dt= −7y sin(x+ 2t) y(0) = 1
dx
dt= 4x− y2 x(0) = 0
De�nindo:
f(t, y, x) = −7y sin(x+ 2t) g(t, y, x) = 4x− y2 t0 = 0 y0 = 1 x0 = 0
Para o primeiro passo de tempo, temos que:
t1 = t0 + ∆t = 0.2
27
Avaliando primeiramente k1 para as duas equações:
ky1 = f(t0, y0, x0) = (−7 · (1) sin(0 + 2 · 0)) = 0
kx1 = g(t0, y0, x0) = (4 · (0)− 12) = −1
Avaliando agora k2
ky2 = f(t0+∆t/2, y0+∆tky1/2, x0+∆tkx1/2) = (−7·(1) sin(−1·0.2/2+2·0.2/2)) = −0.6988339
kx2 = g(t0 + ∆t/2, y0 + ∆tky1/2, x0 + ∆tkx1/2) = (4 · (−1 · 0.2/2)− (1)2) = −1.4
Para k3:
ky3 = ∆tf(t0 + ∆t/2, y0 + ∆tky2/2, x0 + ∆tkx2/2)
= (−7 · (1 + 0.2(−0.6988339)/2) sin((0.2 ∗ (−1.4)/2) + 2 ∗ 0.1)) = −0.3904146283
kx3 = ∆tg(t0 + ∆t/2, y0 + ∆tky2/2, x0 + ∆tkx2/2)
= (4(0.2 · (−1.4)/2)− (1 + 0.2(−0.6988339)/2)2) = −1.42511690
Para k4:
ky4 = f(t1, y0 + ∆tky3 , x0 + ∆tkx3 )
= (−7(1 + 0.2 · (−0.3904146283)) sin((0.2 · (−1.42511690))− 2 · (0.2))) = −0.7403586
kx4 = g(t1, y0 + ∆tky3 , x0 + ∆tkx3 )
= (4(0.2 · (−1.42511690))− (1 + 0.2 · (−0.3904146283))2) = −1.99002461
Agora os pontos y1 e x1 podem ser avaliados:
y1 =∆t
6(ky1 + 2ky2 + 2ky3 + ky4) = 0.9027048
x1 =∆t
6(kx1 + 2kx2 + 2kx3 + kx4 ) = −0.2880086
Repetindo o procedimento para o próximo passo, obtém-se os seguintes coe�cientes:
ky1 = −0.706188 kx1 = −1.9669104 ky2 = −0.6700916 kx2 = −2.6311658
ky3 = −0.285797 kx3 = −2.902888 ky4 = 0.405631 kx4 = −4.1892917
Com isso, obtém-se:
y2 = 0.828960 x2 = −0.8621522
28
2.3.4. Métodos de Runge-Kutta Adaptativos
Como visto anteriormente, o passo de tempo utilizado (∆t) possui uma grande in�uência
no erro associado à resolução de equações diferenciais através de métodos numéricos. Além
disso, o valor de ∆t irá de�nir quantos passos serão necessários para atingir o tempo �nal
desejado. Para reduzir o gasto computacional envolvido, deve-se utilizar o maior passo de
tempo possível que garanta que o erro esteja dentro do limite estabelecido.
Para um grande número de problemas comuns na engenharia, a resolução pode necessitar
de passos muito pequenos em um determinado intervalo de tempo e possibilitar que passos
maiores sejam utilizados para os demais intervalos. Por exemplo, quando a solução representa
um decaimento exponencial, usualmente nos instantes iniciais ocorre uma rápida variação
na solução que tende a estabilizar conforme o tempo avança. Como consequência, passos de
tempo pequenos são necessários nos instantes inicias, porém após um determinado tempo os
passos podem ser aumentados sem que o erro seja superior ao especi�cado.
Para evitar estes inconvenientes, pode-se utilizar um passo de tempo que se adapte a
solução. A implementação deste tipo de algoritmo requer que o erro de truncamento local
seja de alguma forma estimado ao longo da resolução (para cada passo de tempo). Apesar de
isto gerar um gasto computacional extra, normalmente o ganho associado ao uso de passos
adaptativos supera em muito este gasto.
Existem duas formas duas abordagens distintas para implimentar métodos de Runge-
Kutta com passo adaptativo. Uma delas consiste em estimar o erro através da diferença
entre duas predições obtidas com passos diferentes usando um método de mesma ordem.
A outra abordagem consiste em avaliar o erro local de truncamento através da comparação
entre os resultados obtidos com métodos de ordem distintas, utilizado o mesmo passo de
tempo.
Método de Passo Duplo
O método adaptativo mais simples é o método de passo duplo, que consiste em avaliar
cada ponto duas vezes, uma vez através de um único passo e outra através de dois passos
com a metade do tamanho do primeiro. A diferença entre os dois resultados representa
uma estimativa do erro de truncamento local. Por exemplo, quando o ponto yi+1 vai ser
calculado com base no ponto yi, primeiramente calcula-se utilizando um passo ∆t e na
sequência calcula-se novamente yi+1 através de dois passos com tamanho ∆t/2.
29
Suponha que y∗1 representa o valor calculado em t+∆ com o passo ∆t e y∗2 o valor calculado
para o mesmo tempo t+ ∆t mas com o passo ∆t/2, como ilustrado na �gura a seguir.
Por exemplo, considere o método de Runge-Kutta de primeira ordem (Euler explícito).
Como visto anteriormente, neste caso o erro de truncamento local é da ordem de ∆t2.
De�nindo a solução exata em t + ∆t como x, pode-se de�nir a solução exata em termos
de y∗1:
x = y∗1 + ∆t2φ+O(∆t3)
onde φ é um escalar. De forma semelhante, para y∗2 são necessários dois passos, de modo
que:
x = y∗2 + 2
(∆t
2
)2
φ+O(∆t3) = y∗2 +∆t2
2φ+O(∆t3)
Desprezando os termos da ordem de ∆t3, igualando as duas expressões temos que:
y∗1 + ∆t2φ = y∗2 + 2
(∆t
2
)2
φ
De�nindo ∆ = y∗2 − y∗1:
∆ = y∗2 − y∗1 =∆t2
2φ
Assim, o termo ∆ representa uma estimativa do erro associado com y∗2. Para métodos de
maior ordem, o parâmetro ∆ também irá representar uma estimativa do erro associado com
y∗2, porém neste caso o erro não será da ordem de ∆t2 mas irá depender da ordem do método
utilizado. Por exemplo, para o método RK4, o parâmetro ∆ será da orde de ∆t5.
Se o valor ∆ for inferior ao erro local permitido, pode-se aumentar o passo ∆t, enquanto
que se ∆ for maior que o erro permitido, deve-se diminuir ∆t.
30
Considere, por exemplo, que o erro máximo permitido seja emax = 10−6. Se ∆ = |y∗2−y∗1| <
10−6, valida-se o passo e utiliza-se o valor y∗2. Para corrigir o passo de tempo ∆t, diferentes
relações podem ser utilizadas, como por exemplo:
∆t =(emax
∆
)α∆t∗
onde ∆t∗ representa o passo de tempo antigo e ∆t o valor corrigido. Como neste caso
emax > ∆, o valor novo obtido será maior que o anterior. O valor de α recomendado para
este caso (emax > ∆) é de α = 0.2
Considerando agora que ∆ > emax. Neste caso, o passo é invalidado e deve ser recalculado
com um passo menor. A relação utilizada para recalcular o passo é a mesma apresentada
anteriormente, porém neste caso o valor corrigido será menor que o anterior e o valor reco-
mendado para α é de α = 0.25.
Exemplo 02) Utilize o método de Runge-Kutta de quarta ordem juntamente com o
método do passo duplo para estimar o valor de y(1.25). Considere que o erro de truncamento
local não pode ser superior a 10−3 e como estimativa inicial para o passo de tempo considere
∆t = 1:dy
dt= 4e0.8t − 0.5y y(0) = 2
Avaliando o primero passo com ∆t = 1, obtemos y(1) = y∗1 = 6.201. Considerando agora
∆t = 0.5, a resolução de dois passos implica que y(1) = y∗2 = 6.195. Assim:
∆ = y∗2 − y∗1 = 0.006
Como o valor é superior ao erro máximo, deve-se refazer este passo. Recalculando o passo
de tempo:
∆t =
(10−3
0.006
)0.25
1 = 0.638943
Recalculando o primeiro passo de tempo com o novo valor de ∆t, obtém-se que y∗1 = 4.3481
e y∗2 = 4.3475, de modo que ∆ = 6 × 10−4. Como o valor é inferior ao erro máximo, este
ponto pode ser validado:
t1 = t0 + ∆t = 0.638943
y(t1) = y1 = 4.3475
Para avaliar o próximo passo, primeiramente pode-se avaliar um novo ∆t:
∆t =
(10−3
6× 10−4
)0.2
· 0.638943 = 0.707672
31
Utilizando este valor, obtém-se qu para t2 = t1 + ∆t = 0.638943 + 0.707672 = 1.346615,
y∗1 = 8.4886328 e y∗2 = 8.4869326, de modo que ∆ = 0.0017. Como o valor está acima do
erro máximo, deve-se reduzir o passo de tempo:
∆t =
(10−3
0.0017
)0.25
0.707672 = 0.619755
Recalculado os valores, agora para t2 = t1 +∆t = 0.638943+0.619755 = 1.258698, obtém-
se que y∗1 = 7.849350209 e y∗2 = 7.84849122, de modo que ∆ = 8.59 × 10−4. Como o erro
está abaixo do especi�cado, pode-se validar o passo. Assim:
t2 = 1.258698
y(t2) = y2 = 7.84849122
Como deseja-se saber y(1.25), o valor de t2 já está acima do especifado. Para obter o valor
no ponto t = 1.25, pode-se realizar uma interpolação linear entre os valores para t1 e t2:
y(1.25) = y1 +(t2 − 1.25)
(t2 − t1)(y2 − y1) = 7.79936
Método de Runge-Kutta-Fehlberg
Outra alternativa para estimar o erro local de truncamento é avaliar um dado ponto através
de métodos com ordem distinta. Por exemplo, pode-se avaliar a solução em um dado ponto
y(ti+1) através de um método de segunda ordem (y∗i+1) e de um método de terceira ordem
(yi+1). A diferença entre os dois valores (yi+1 − y∗i+1) representa uma estimativa do erro de
truncamento neste ponto. Se os dois valores forem muito próximos, isto indica que a redução
no passo de tempo não irá diminuir signi�cativamente o erro local de truncamento.
A abordagem mais utilizada consiste em avaliar a solução em um ponto através de métodos
de Runge-Kutta de quarta e quinta ordem, sendo que as estimativas de quarta ordem (y∗i+1)
e de quinta ordem (yi+1) podem ser representadas como:
y∗i+1 = yi + ∆t
(37
378k1 +
250
621k3 +
125
594k4 +
512
1771k6
)yi+1 = yi + ∆t
(2825
27648k1 +
18575
48384k3 +
13525
55296k4 +
277
14336k5
1
4+ k6
)onde:
k1 = f(ti, yi)
32
k2 = f
(ti +
∆t
5, yi +
∆t
5k1
)k3 = f
(ti +
3∆t
10, yi +
3∆t
40k1 +
9∆t
40k2
)k4 = f
(ti +
3∆t
5, yi +
3∆t
10k1 −
9∆t
10k2 +
6∆t
5k3
)k5 = f
(ti + ∆t, yi −
11∆t
54k1 +
5∆t
2k2 −
70∆t
27k3 +
35∆t
27k4
)k6 = f
(ti +
7∆t
8, yi +
1631∆t
55296k1 +
175∆t
512k2 +
575∆t
13824k3 +
44275∆t
110592k4 +
253∆t
4096k5
)Da mesma forma como para o método anterior, o erro associado com a estimativa do ponto
yi+1 é avaliado como:
∆ = yi+1 − y∗i+1
Com base neste valor, utiliza-se o mesmo critério apresentado anteriormente para corrigir
o passo de tempo ∆t, ou seja, se o erro for menor que o especi�cado aumenta-se o passo de
tempo e se o erro for menor descarta-se o valor obtido e reduz-se o passo de tempo até que
o erro esteja dentro do valor tolerável. As relações utilizadas para reduzir ou aumentar o
passo são as mesmas apresentadas para o método do passo duplo.
33
2.4. Métodos de Passos Múltiplos
Os métodos de passos múltiplos são aqueles onde a variável yn+1 é determinada com base
nos valores desta variável em mais de um ponto, como por exemplo no método baseado na
regra do ponto médio:
yn+1 = yn−1 + 2∆ty′(tn) → y′(tn) =yn+1 − yn−1
2∆t
Os demais métodos, que determinam yn+1 com base somente nos valores de yn, são cha-
mados de métodos de passo único.
Como visto anteriormente, os métodos de passos múltiplos não iniciam automaticamente,
sendo necessário avaliar alguns pontos iniciais através de outro método. Dentre os métodos
de passos múltiplos, os mais utilizados são os da família dos métodos de Adams, como
apresentado a seguir.
2.4.1. Método de Adams-Bashforth
Considere novamente o seguinte PVI:
dy
dt= f(t, y) y(t0) = y0
Uma maneira de obter o valor em um ponto yn+1 com base em um valor conhecido yn é
integrar a equação desde um ponto tn até um ponto tn+1:∫ tn+1
tn
y′(t)dt =
∫ tn+1
tn
f(t, y)dt
Como o lado esquerdo avalia a integral de uma derivada, podemos escrever a expressão
anterior como:
y(tn+1)− y(tn) = yn+1 − yn =
∫ tn+1
tn
f(t, y)dt
Como a função y(t) não é conhecida, não é possível resolver diretamente a integral do lado
direito da equação. O método de Adams-Bashforth consiste em substituir a função f(t, y)
por um polinômio p(t), permitindo assim a resolução da integral:
yn+1 = yn +
∫ tn+1
tn
p(t)dt
Dependendo da ordem escolhida para o polinômio p(t), diferentes formulações são obtidas.
Os coe�cientes associados ao polinômio são determinados com base nos valores já conhecidos
34
para yn, yn−1, yn−2, etc. Caso um polinômio de ordem k for utilizado, é necessário conhecer
a solução em k + 1 pontos.
Por exemplo, caso um polinômio de primeira ordem for empregado:
p(t) = At+B
é necessário conhecer o valor da função f(x, y) em dois pontos, ou seja, é necessário deter-
minar dois valores (yn e yn−1) para encontar as constantes A e B. Com isso:
Atn +B = f(tn, yn) = fn
Atn−1 +B = A(tn −∆t) +B = f(tn−1, yn−1) = fn−1
Resolvendo para A e B:
A =fn − fn−1
∆tB =
fn−1tn − fntn−1
∆t
Substituindo p(t) = At+B na expressão anterior e resolvendo a integral:
yn+1 = yn +
∫ tn+1
tn
(At+B)dt = yn +A
2(t2n+1 − t2n) +B(tn+1 − tn)
Substituindo as expressões obtidas para A e B e simpli�cando o resultado obtido:
yn+1 = yn +3
2∆t · fn −
1
2∆t · fn−1
Esta relação é conhecida como método de Adams-Bashforth de 2a ordem, e representa um
método explícito onde é necessário conhecer a solução nos dois pontos anteriores.
Usando polinômios de maior ordem, pode-se obter métodos de maior ordem. Um dos mais
utilizados é o de quarta ordem:
yn+1 = yn +∆t
24(55fn − 59fn−1 + 37fn−2 − 9fn−3)
Neste caso, é necessário conhecer 4 pontos anteriores. Em comparação com outros métodos
de quarta ordem, como RK4, o método de Adams-Bashforth costuma apresentar melhores
resultados.
2.4.2. Método de Adams-Moulton
O método de Adams-Moulton é uma modi�cação do método de Adams-Bashftorh que con-
siste em avaliar o polinômio interpolador p(t) em tn+1, tn, tn−1, . . . ao invés de tn, tn−1, tn−2, . . ..
35
Por exemplo, considere novamente um polinômio de primeira ordem p(t) = αt+ β. Neste
caso, as constantes são avaliadas através da relações:
αtn+1 + β = fn+1
αtn + β = fn
De onde se obtém que:
α =fn+1 − fn
∆tβ =
fntn+1 − fn+1tn∆t
De modo que yn+1 pode ser obtido como:
yn+1 = yn +∆t
2fn +
∆t
2fn+1
Considerando que fn+1 = f(tn+1, yn+1), esta relação representa uma fórmula implícita de
segunda ordem.
Utilizando polinômios de ordem maior, obtém-se formulações de ordem superior. O mé-
todo de Adams-Moulton de quarta ordem é dado como:
yn+1 = yn +∆t
24(9fn+1 + 19fn − 5fn−1 + fn−2)
Para evitar as di�culdades associadas ao uso de métodos implícitos, é comum adotar
uma abordagem preditiva-corretiva, utilizando Adams-Bashforth como uma etapa preditiva
e Adams-Moulton como corretora.
Por exemplo, utilizando uma abordagem de segunda ordem, a etapa preditiva é dada
como:
y∗n+1 = yn +3
2∆t · fn −
1
2∆t · fn−1
Sendo este valor posteriormente corrigido com a relação de Adams-Moulton:
yn+1 = yn +∆t
2fn +
∆t
2f ∗n+1
onde f ∗n+1 = f(tn+1, y∗n+1)
Exemplo 03:) Utilizando uma abordagem preditiva-corretiva, utilize os métodos de
Adams de segunda ordem para obter a solução aproximada do seguinte PVI até t = 0.4,
utilizando ∆t = 0.1:dy
dt= −y + 2t y(0) = 2
36
Esta equação foi avaliada anteriormente com o método de Runge-Kutta. Para utilizar
Adams-Bashforth, neste caso é necessário conhecer a solução em dois pontos anteriores.
Pode-se utilizar a condição inicial como um destes pontos:
t0 = 0 y0 = 2
e o valor obtido com o método de Runge-Kutta para y1:
t1 = 0.1 y1 = 1.82
Com isso, pode-se determinar o valor preditivo em t2 = 0.2:
y∗2 = y1 +3
2∆t · f(t1, y1)− 1
2∆t · f(t0, y0) = 1.577
Utilizando Adams-Moulton para corrigir o valor:
y2 = y1 +∆t
2f(t1, y1) +
∆t
2f(t2, y
∗2) = 1.68015
Repetindo o procedimento para o próximo passo de tempo (t3 = 0.3):
y∗3 = y2 +3
2∆t · f(t2, y2)− 1
2∆t · f(t1, y1) = 1.48313
y3 = y2 +∆t
2f(t2, y2) +
∆t
2f(t3, y
∗3) = 1.57199
Finalmente, para o último passo de tempo (t4 = 0.4):
y∗4 = y3 +3
2∆t · f(t3, y3)− 1
2∆t · f(t2, y2) = 1.4902
y4 = y3 +∆t
2f(t3, y3) +
∆t
2f(t4, y
∗4) = 1.4889
37
2.5. Estabilidade dos Métodos de Resolução de EDO's
Na passagem do método de Euler explícito para métodos mais complexos, o principal
objetivo foi obter métodos com melhor precisão (reduzir os erros de truncamento). No
entanto, em muitos casos os resultados obtidos não só possuem uma baixa precisão como
também são catastro�camente distintos da solução exata.
Por exemplo, considere o seguinte PVI:
dy
dt= −12y y(0) = 1
A utilização do método de Euler explícito com ∆t = 0.2 gera o resultado ilustrado a seguir.
Neste caso, o erro aumenta rapidamente conforme se avança na solução. Este problema
está associado com a falta de estabilidade do método empregado.
Antes de avaliar a estabilidade dos métodos, é necessário de�nir dois conceitos relacionados
a análise de estabilidade:
Consistência: Um método de passo único (como os da família de Runge-Kutta) é dito
consistência se o erro de truncamento local ei tende a zero conforme ∆t → 0 para todos os
passos de tempo, ou seja:
lim∆t→0
max1≤i≤N
|ei| = 0
onde N representa o número total de iterações.
Convergência: Um método de passo único é dito convergente com respeito ao pro-blema
de valor inicialdy
dt= f(t, y) y(t0) = y0
se a diferença entre a solução obtida através do método numérico yi e a solução exata no
mesmo ponto φi tender a zero conforme ∆t→ 0:
lim∆t→0
max1≤i≤N
|yi − φi| = 0
38
Dessa forma, um método consistente possui a propriedade de que as equações obtidas para
cada passo se aproximada da própria equação diferencial conforme o passo de tempo tende
a zero, de modo que o erro de truncamento local também tende a zero.
Porém, além do erro de truncamento, deve-se considerar a in�uência do erro de arredonda-
mento, associado ao fato de que os valores numéricos não são representados de forma exata.
Na prática, os parâmetros do sistema, as condições inicias e todo o processo aritmético sub-
sequente possuem erros associados com a precisão numérica �nita dos valores empregados.
Para garantir que um determinado método seja convergente, deve-se então garantir que
além de ser consistente, o erro de arredondamento deve �car limitado a um valor aceitável.
O controle do erro de arredondamento está associado com a estabilidade do método. Este
conceito é muito similar ao conceito de condicionamento de um sistema linear, no sentido
de que um método é dito estável quando pequenas variações nos parâmetros ou condições
iniciais levam a igualmente pequenas variações na solução obtida. A seguir será apresentada
uma análise da estabilidade para os métodos de passo único. Para os métodos de passos-
múltiplos, como existem diversas etapas de aproximação envolvidas em cada passo, a análise
é mais complexa.
2.6. Análise de Estabilidade de Métodos de Passo Único
Para avaliar a estabilidade de um dado método, considere o seguinte PVI:
dy
dt= f(t, y) y(t0) = y0
Considere que a solução exata deste PVI em um ponto tn é dada por φ(tn) e que a solução
obtida com o uso de método numérico neste ponto é yn. Como comentado anteriormente, na
resolução computacional do problema, existem associados erros de arredondamento devido à
precisão limitada da representação numérica. Considere que y∗n é o valor arredondado obtido
em uma máquina real e yn é o valor com precisão in�nita obtido em uma máquina ideal.
O erro total associado será a diferença entre a solução exata e o valor fornecido pelo
computador:
Erro total = φ(tn)− y∗n = (φ(tn)− yn) + (yn − y∗n)
O primeiro termo representa o erro de truncamento e o segundo o erro de arredondamento
associado ao ponto tn.
Para que um dado método numérico seja adequado, são necessárias duas condições:
39
- O erro de truncamento acumulado (global) deve tender a zero conforme ∆t→ 0 (consis-
tência);
- O erro de arredondamento acumulado (que não pode ser eliminado) deve ser pequeno
em comparação com a solução exata (estabilidade).
2.6.1. Equação Modelo: Decaimento de Primeira Ordem
Para ilustrar as características de estabilidade de um método, considere a equação:
dy
dt= −λy y(0) = 1
onde λ > 0. Esta equação representa uma taxa de decaimento de primeira ordem, que
surge, por exemplo, na análise da variação na fração de um dado reagente em uma reação
de primeira ordem.
Considerando novamente que φ(t) é a solução exata da equação, o valor calculado nume-
ricamente pode ser expresso como:
y(t) = φ(t) + ε(t)
Substituindo na equação diferencial:
d(φ+ ε)
dt= −λ(φ+ ε)
Como as solução exata deve satisfazer a equação diferencial, temos que:
dε
dt= −λε
Esta equação serve como uma relação para determinar como o erro varia ao longo do tempo,
ou seja, ao longo dos passos avaliados. Se o erro diminuir com o tempo, o método é estável,
caso contrário será instável.
Por exemplo, considere que o método de Euler explícito seja utilizado para avaliar ε(t):
εn+1 = εn + ∆t(−λεn) = εn(1−∆tλ) → εn+1
εn= 1−∆tλ
Para garantir a estabilidade, é preciso que o erro em tn+1 seja menor que o erro em tn,
assim: ∣∣∣εn+1
εn
∣∣∣ ≤ 1
40
Considerando a relação anterior:
|1−∆tλ| ≤ 1
Desse modo, o passo de tempo para garantir a estabilidade deve ser maior que zero e
menor que 2/λ.
A região contendo os valores de λ∆t que levam a uma solução estável é chamado de domínio
de estabilidade linear. O parâmetro λ pode ser um número complexo, de modo que o domínio
de estabilidade costuma ser representado no plano complexo. De�nindo z = −λ∆t, temos
que para o método de Euler explícito ser estável a seguinte condição deve ser satisfeita:
|1 + z| ≤ 1
Fazendo z = a+ bi:
|1 + (a+ bi)| ≤ 1 → |(a+ 1) + bi| ≤ 1
Para um número complexo qualquer α + βi, o módulo é de�nido como:
|α + βi| =√α2 + β2
Assim, para o caso anterior, temos que:
|(a+ 1) + bi| =√
(a+ 1)2 + b2 ≤ 1 → (a+ 1)2 + b2 ≤ 1
Esta relação representa um círculo de raio menor ou igual a 1, deslocado em uma unidade
para a esquerda no eixo equivalente a a. Como a é a parte real de z e b a parte imaginária
de z, o domínio de estabilidade do método de Euler explícito representa a região indicada
na �gura a seguir.
Considere agora que seja empregado o método de Euler implícito:
εn+1 = εn + ∆t(−λεn+1) → εn+1(1 + ∆tλ) = εn
Com isso:εn+1
εn=
1
1 + ∆tλ
Como λ > 0 e ∆t > 0, esta relação mostra que o método de Euler implícito é estável para
qualquer valor de ∆t, o que é uma característica comum dos métodos implícitos.
41
Considere agora a modi�cação de segunda ordem baseada na regra do ponto médio:
εn+1 = εn +∆t
2(−λεn − λεn+1)
De modo que:
εn+1(1 + ∆tλ/2) = εn(1−∆tλ/2) → εn+1
εn=
1−∆tλ/2
1 + ∆tλ/2
Considerando a condição para estabilidade:∣∣∣1−∆tλ/2
1 + ∆tλ/2
∣∣∣ ≤ 1
Neste caso, a solução também será estável para qualquer ∆t > 0.
Para os demais métodos explícitos, a estabilidade também está condicionada a um domínio
especí�co. Em particular para o caso dos métodos de RK de ordem maior que 1, pode-se
mostrar que o domínio de estabilidade linear engloba o domínio associado ao método de
Euler explícito (RK de primeira ordem). Assim, se a condição ∆t < 2/λ for respeitada, os
métodos de RK de ordem superior serão estáveis. O domínio de estabilidade de métodos de
Runge-Kutta de ordem 1 até 5 é representado na �gura a seguir:
2.7. Problemas Rígidos (Sti�)
Algumas equações ou sistemas de equações diferenciais são classi�cados como rígidos
(sti�). Não existe uma de�nição precisa para classi�car uma EDO como rígida, mas es-
tes tipos de problema compartilham características em comum:
42
- Normalmente existem termos que levam à uma rápida variação na solução;
- Problemas rígidos possuem uma variedade de escalas de tempo associadas, ou seja, em
determinados pontos a solução varia muito mais rapidamente que em outras;
- Métodos explícitos só são estáveis para a resolução de EDO's rígidas se o passo de tempo
utilizado for muito pequeno;
- Usualmente, o passo de tempo utilizado para garantir a estabilidade é menor que o
necessário para garantir a convergência desejada.
A equação teste utilizada anteriormente:
dy
dt= −λy y(t0) = y0
é um exemplo de equação rígida, especialmente para altos valores de λ. Esta equação possui
solução da forma:
y = y0e−λt
43
Como visto anteriormente, se um método explícito for empregado, deve-se usar um valor de
∆t su�cientemente pequeno para garantir a estabilidade da solução. Por exemplo, para o
método de Euler explícito:
∆t <2
λ
Caso for necessário avaliar a solução até um tempo �nal t1, o número de passos necessários
(n) será:
n =t1∆t
→ n >λt12
Assim, o número de passos mínimo necessários é diretamente proporcional ao valor de λ.
2.7.1. Sistemas de Equações Diferenciais Rígidas Lineares
Considere o seguinte PVI:
dy1
dt= −500.5y1 + 499.5y2 y1(0) = 2
dy2
dt= 499.5y1 − 500.5y2 y2(0) = 1
A solução deste problema é a seguinte:
y1(t) = 1.5e−t + 0.5e−1000t
y2(t) = 1.5e−t − 0.5e−1000t
Assim, a solução possui dois termos: um termo e−t que varia lentamente com o tempo e
outro e−1000t que varia rapidamente. Para garantir a estabilidade, é preciso que a solução
se mantenha estável durante toda a resolução, por isso, a estabilidade deve ser assegurada
para o termo e−1000t.
De modo geral, os valores de ∆t empregados devem ser determinados sem o conhecimento
da solução. O PVI anterior pode ser reescrito como:
dY
dt=
−500.5 499.5
499.5 −500.5
Y = AY
O valor de ∆t necessário é de�nido com base nos autovalores da matriz A. Lembrando
que os autovalores λ são de�nidos como os valores onde:
|A− λI| = 0
44
Assim, para este caso:∣∣∣∣∣−500.5− λ 499.5
499.5 −500.5− λ
∣∣∣∣∣ = 0 → (−500.5− λ)2 − 499.52 = 0
Resolvendo a equação de segundo grau, obtém-se as seguintes raízes:
λ1 = −1 λ2 = −1000
Para garantir a estabilidade em um método explícito, pode-se considerar a função teste
usada anteriormente dy/dt = −λy, com λ sendo o maior dos autovalores (em módulo)
associados com o problema. Por exemplo, caso o método de Euler explícito seja empregado,
deve-se garantir que:
∆t ≤ 2
|λmax|Considerando que no exemplo |λmax| = 1000, então temos que o passo de tempo máximo
é de 0.002. Caso o menor autovalor fosse empregado, seria obtido um valor de ∆tmax = 2.
Após os instantes iniciais, o problema passa a ser controlado pelo menor autovalor.
Para medir a importância da rigidez na resolução do problema, pode-se determinar o grau
de rigidez (sti�ness ratio) do problema, de�nido com:
SR =|λmax||λmin|
Normalmente, para SR > 20 o problema já é classi�cado como rígido, sendo necessário
avaliar com cuidado os passos de tempo empregados.
2.7.2. Sistemas Não-Lineares
Para problemas não-lineares, a análise é mais complexa. Para problemas autônomos, da
forma:dY
dt= f(Y )
pode-se linearizar a equação através de uma expansão em série de Taylor em torno do ponto
tn. Desconsiderando os termos de alta ordem:
dY
dt= f(Yn) + J(tn)(Y − Yn)
onde J(tn) é a matriz Jacobiana avaliada em t = tn, de�nida como:
J(tn) = aij =
[∂fi(y)
∂yj
]tn
45
Por exemplo, em um sistema 2× 2 da forma:
dy
dt= f(x, y)
dx
dt= g(x, y)
O Jacobiano em um ponto tn é determinado da forma:
J(tn) =
∂f∂x
∣∣∣tn
∂f∂y
∣∣∣tn
∂g∂x
∣∣∣tn
∂g∂y
∣∣∣tn
Neste caso, a razão de rigidez é de�nida com base nos autovalores da matriz Jacobiana.
Como J(tn) pode depender do tempo, a razão de rigidez pode variar ao longo da solução.
Os valores de ∆t necessários para garantir a estabilidade também podem ser de�nidos com
base nos autovalores da matriz Jacobiana.
Obs.: Caso tn for um ponto �xo, os autovalores do Jacobiano também servem para de�nir
se o ponto é estável ou não, sendo instável caso algum autovalor tenha parte real positiva.
Exemplo: Considere o seguinte PVI:
dy
dt= y2x y(0) = 2
dx
dt= −400yx2 x(0) = 1
Determine a razão de rigidez em t = 0. Caso o método de Euler explícito for empregado,
qual o valor máximo de ∆t?
46
3. Métodos Numéricos para Problemas
de Valor de Contorno
Equações diferenciais de ordem maior que um podem gerar problemas de valor inicial
(PVI) ou problemas de valor de contorno (PVC), dependendo da forma como as condições
conhecidas são especi�cadas. Por exemplo, considere a EDO:
y′′ + p(x)y′ + q(x)y = g(x)
Até o momento, foram analisados principalmente casos onde condições iniciais conhecidas
são especi�cadas da forma:
y(x0) = y0 y′(x0) = y′0
originando desta forma um PVI. Para a resolução numérica de PVI's, pode-se partir da
condição inicial e ir avançando até um tempo �nal arbitrário.
Em muitos casos, os problemas envolvem condições conhecidas em pontos diferentes, sendo
estas chamadas de condições de contorno, podendo ser expressas, por exemplo, como:
y(α) = y0 y(β) = y1
De forma geral, os PVC's envolvem uma coordenada espacial como variável independente.
Assim, a resolução de um PVC's consiste em buscar uma solução que satisfaz a equação
diferencial no intervalo α < x < β juntamente com as condições de contorno especi�cadas.
Isto implica que existem duas condições, em pontos diferentes do domínio de solução, que
deve ser simultaneamente satisfeitas. Por isso, os métodos de marcha (como os de Runge-
Kutta) não podem ser empregados neste caso.
47
3.1. Estratégias de Solução de PVC's
Os métodos para resolver problemas de valor de contorno se dividem em duas categorias:
os baseados em transformar o PVC em um PVI e os baseados em discretizar a equação
utilizando métodos de diferenças �nitas.
Os métodos que transformação PVC's e PVI's consistem basicamente em utilizar uma das
condições de contorno como condição inicial e assumir (chutar) diferentes valores para uma
segunda condição inicial de modo que o resultado obtido satisfaça a segunda condição de
contorno. Por exemplo, considere a seguinte equação utilizada para descrever a variação na
temperatura em uma barra que perde calor para o ambiente por convecção, como apresentado
na �gura a seguir.
Considerando que a barra seja muito �na, com raio muito menor que o comprimento,
pode-se assumir que a equação que descreve a variação na temperatura ao longo de x pode
ser expressa como:d2T
dx2+ h′(Ta − T ) = 0
onde h′ é um coe�ciente de troca térmica e Ta é a temperatura ambiente.
As condições de contorno corresponde a temperatura �xas nas extremidades, em x = 0 e
x = L, de modo que:
T (0) = T1 T (L) = T2
Para transforma a equação em um PVI equivalente, primeiramente é preciso aplicar o
método de redução de ordem. Escrevendo o PVI como:
dT
dx= z T (0) = T1
dz
dx+ h′(Ta − T ) = 0 z(0) = z0
48
O valor de z0 não é conhecido, pois o valor da derivada de T em x = 0 não foi especi�cado.
O método utilizado neste caso consiste em assumir diferentes valores para z0 até que a con-
dição T (L) = T2 seja satisfeita. Para equação lineares, este método funciona razoavelmente
bem, pois pode-se interpolar os valores para T (L) obtidos para diferentes valores de z0 para
encontrar o valor de z0 que corresponde a solução correta do problema. No entanto, de
modo geral este método é pouco utilizado, em particular porque os métodos baseados em
aproximações por diferenças �nitas costumam ser mais e�cientes, como ser apresentado a
seguir.
3.2. Aproximações por Diferenças Finitas
O método de diferenças �nitas é um método de discretização de equações diferenci-
ais. Isto signi�ca que ele transforma uma função contínua em uma representação dis-
creta (pontos). Por exemplo, considere uma função f(x) = 2x de�nida em um intervalo
entre 0 e 1. Esta função pode ser representada de forma discreta como, por exemplo,
f [x] = [0, 0.4, 0.8, 1.2, 1.6, 2], assumindo um espaçamento entre os pontos de 0.2, ou seja, esta
representação está relacionada como um domínio discreto da forma x = [0, 0.2, 0.4, 0.6, 0.8, 1].
As soluções obtidas com a aplicação do método de diferenças �nitas sempre serão discretas.
Caso for necessário obter os valores para algum valor de x que não corresponda exatamente
aos pontos do domínio, pode-se interpolar os valores.
A primeira etapa da aplicação do método de diferenças �nitas consiste exatamente em
de�nir o domínio discreto onde a solução será buscada. Por exemplo, considere o caso
apresentado anteriormente para a distribuição de calor em uma barra estacionária. A região
onde se deseja obter a distribuição de temperatura é no intervalo de x = 0 até x = L. Este
intervalo corresponde ao domínio de solução da equação diferencial. No entanto, ele está
em uma forma contínua e não discreta. Para discretizar o domínio, deve-se dividi-lo em um
determinado número de pontos. Por exemplo, considere que L = 1 e que se deseja dividir o
domínio em pontos com espaçamento ∆x = 0.1, ou seja, deseja-se dividir o domínio em 10
elementos de igual tamanho. Quanto mais elementos forem utilizados, maior será a precisão
do método, porém o gasto computacional também irá aumentar.
O domínio discreto será então:
x = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9.1]
49
Como pode ser visto, este conjunto contém 11 pontos. De maneira geral, o número de pontos
sempre será igual ao número de elementos em que o domínio é dividido mais 1. Este vetor
pode ser representado de uma forma mais simples como:
x[i] = i∆x i = 0, 1, 2, . . . 10
Na resolução de PVI's, não é necessário inicialmente de�nir o domínio de solução, pois
pode-se continuar avançando por quantos passos forem necessários, partindo de um valor
inicial. No entanto, para o caso de PVC's, o domínio de solução é fechado e deve ser
considerado como um todo. Dependendo da formulação utilizada, pode até mesmo ser
necessário obter os valores para todos os pontos ao mesmo tempo.
A estratégia do método de diferenças �ntas consiste em buscar equações algébricas que
aproximem a solução em cada ponto i. Para isso, as derivadas são aproximadas como relações
algébricas envolvendo a solução em diferentes valores de i. Esta aproximação pode ser
realizada de diferentes formas, dependendo da precisão desejada e da natureza do problema e
das condições de contorno. Porém, a origem destas aproximações sempre é uma aproximação
em série de Taylor em torno de cada ponto i, como será discutido a seguir.
3.2.1. Aproximação da Derivada Primeira
Para apresentar o processo de discretização da derivada de uma função contínua y(x) em
um intervalo 0 ≤ x ≤ L, será considerado um domínio discreto como o apresentado na �gura
a seguir, onde o domínio físico contínuo (região entre 0 e L) é dividido em N + 1 pontos.
Lembrando novamente, o objetivo do método de diferenças �nitas é obter aproximações para
o valor da função y(x) em cada um destes N + 1 pontos.
Esta representação discreta do domínio de solução é normalmente chamada de grid nu-
mérico ou malha numérica.
Como visto em aulas anteriores, a expansão em série de Taylor pode ser utilizada para
avaliar o valor de uma função em um dado ponto com base no valor conhecido em outro ponto.
Dependendo da forma como a aproximação é realizada, obtém-se diferentes formulações para
o método de diferenças �nitas, como será apresentado a seguir.
50
Aproximação para Frente (Foward)
Considere que se deseja aproximar o valor em xi+1 com base no valor em xi, ou seja,
deseja-se aproximar y(xi+1) = yi+1 com base em y(xi) = yi.
A expansão em série de Taylor neste caso pode ser expressa como:
yi+1 = yi + (xi+1 − xi)dy
dx
∣∣∣x=xi
+(xi+1 − xi)2
2!
d2y
dx2
∣∣∣x=xi
+(xi+1 − xi)3
3!
d3y
dx3
∣∣∣x=xi
+ . . . (3.1)
Considerando novamente que (xi+1 − xi) seja relativamente pequeno, os termos de alta
proporcionais a (xi+1 − xi)2, (xi+1 − xi)
3, . . . podem ser desprezados. Assim, a derivada
primeira da função y(x) no ponto xi pode ser aproximada como:
dy
dx
∣∣∣x=xi
=yi+1 − yixi+1 − xi
De�nindo ∆x = xi+1 − xi, a expressão pode ser dada por:
dy
dx
∣∣∣x=xi
=yi+1 − yi
∆x
Esta expressão é conhecida como aproximação por diferenças �nitas para frente (ou método
de diferenças �nitas para frente), pois utiliza o valor da função em um ponto a frente xi+1
para estimar o valor da derivada em um ponto anterior xi. Fazendo o limite de ∆x tendendo
a zero:dy
dx
∣∣∣x=xi
= lim∆x→0
yi+1 − yi∆x
obtém-se a própria de�nição da derivada em um ponto. Portanto, esta formulação é con-
sistente.
Como pode ser visto, este método é equivalente ao método de Euler explícito, sendo que
a diferença é que no método de Euler deseja-se aproximar o valor da função em diferentes
pontos com base no conhecimento da derivada dy/dt = f(t, y) e no método de diferenças
�nitas deseja-se aproximar o valor da derivada com base no valor em diferentes pontos. Nos
dois casos, porém, ocorre uma linearização da função em torno de um ponto.
Na aproximação da derivada, foram desprezados os termos O(xi+1− xi)2, assim, o erro de
truncamento local do método de diferenças para frente é da ordem de O(∆x2). De forma
semelhante ao apresentado para o método de Euler, pode-se mostrar que quando aplicado
para avaliar N pontos, o erro associado será da ordem de O(xi+1 − xi) = O(∆x), ou seja, a
aproximação por diferenças para frente é um método de primeira ordem.
51
Aproximação para Trás (Backward)
De forma semelhante ao realizado para obter a derivada em xi com base na expansão para
obter yi+1, pode-se realizar uma expansão para obter a função em xi−1:
yi−1 = yi + (xi−1 − xi)dy
dx
∣∣∣x=xi
+(xi−1 − xi)2
2!
d2y
dx2
∣∣∣x=xi
+(xi−1 − xi)3
3!
d3y
dx3
∣∣∣x=xi
+ . . .
Considerando que o espaçamento entre os pontos é constante, temos novamente que ∆x =
xi − xi−1 = −(xi−1 − xi), assim:
yi−1 = yi − (∆x)dy
dx
∣∣∣x=xi
+(∆x)2
2!
d2y
dx2
∣∣∣x=xi− (∆x)3
3!
d3y
dx3
∣∣∣x=xi
+ . . . (3.2)
Novamente, desprezando os termos de ordem maior ou igual a 2, obtém-se:
dy
dx
∣∣∣x=xi
=yi − yi−1
∆x
esta aproximação é conhecida como aproximação por diferenças �nitas para trás (ou mé-
todo de diferenças �nitas para trás), e também representa uma aproximação de primeira
ordem para a derivada em um dado ponto xi.
A utilização dos métodos para trás ou para frente é, a princípio, equivalente. A única
restrição ocorre nos pontos extremos do domínio. Por exemplo, no ponto x0 não pode ser
aplicado o método para trás pois não existe um ponto anterior a este. De forma semelhante,
no ponto xN a formulação para frente não pode ser utilizada, pois de maneira equivalente
não existe nenhum ponto após este para ser utilizado de base.
Aproximação Central
De forma geral, a utilização de métodos de primeira ordem para a discretização de todos
os pontos do domínio não costuma apresentar bons resultados, especialmente nos casos
envolvendo gradientes aproximadamente simétricos em relação a direção x, como por exemplo
problemas envolvendo condução de calor ou difusão de massa.
Uma aproximação de segunda ordem pode ser obtida fazendo a expansão para yi+1 (Eq.
1) menos a expansão para yi−1 (Eq. 2). Com isso, obtém-se:
yi+1 − yi−1 = 2∆xdy
dx
∣∣∣x=xi
+2(∆x)3
3!
d3y
dx3
∣∣∣x=xi
+ . . .
Desprezando agora os termos da ordem de O(∆x)3, pode-se obter a seguinte expressão para
uma aproximação da derivada primeira em xi:
dy
dx
∣∣∣x=xi
=yi+1 − yi−1
2∆x
52
Esta expressão é conhecida como aproximação pode diferenças �nitas central (ou método de
diferenças �nitas central). Neste caso, o erro de truncamento local é da ordem de O(∆x)3,
de modo que o erro global será da ordem de O(∆x)2. Assim, esta aproximação é de segunda
ordem.
Existem aproximações de ordem superior, porém a aproximação central é a mais empre-
gada, sendo adequada para a maioria dos casos. Novamente, dever-se observar que esta
aproximação não pode ser aplicada nos pontos extremos do sistema. A melhor estratégia
para a discretização da equação é utilizar o método central para os pontos internos, o método
para frente no ponto x0 e o método para trás no ponto xN .
Uma comparação geométrica dos três esquemas de discretização é apresentada na �gura
a seguir.
3.2.2. Aproximação da Derivada Segunda
O método de diferenças �nitas pode ser aplicado também para a discretização de derivadas
de maior ordem. As aproximações para a derivada segunda são obtidas considerando a
de�nição da derivada em um ponto. Da mesma forma que a derivada primeira pode ser
aproximada em termos da variação da função em dois pontos, a derivada segunda pode
ser aproximada em termos da variação na derivada primeira em dois pontos. Por exemplo,
utilizando um esquema para frente:
d2y
dx2
∣∣∣x=xi
=
dy
dx
∣∣∣x=xi+1
− dy
dx
∣∣∣x=xi
xi+1 − xi53
onde a derivada avaliada no ponto xi+1, utilizando novamente um esquema para frente, é
obtida de forma equivalente a derivada no ponto xi como:
dy
dx
∣∣∣x=xi+1
=yi+2 − yi+1
xi+2 − xi+1
Considerando que os pontos sejam igualmente espaçados (grid igualmente espaçado), pode-se
juntar as expressões para de derivada em xi+1 e a expressão para a derivada em xi, obtém-se
a seguinte expressão para a derivada segunda em xi:
d2y
dx2
∣∣∣x=xi
=yi − 2yi+1 + yi+2
(∆x)2
Esta expressão representa o esquema para frente aplicado à derivada segunda. Como ele se
baseia em esquemas de primeira ordem, também é uma aproximação de primeira ordem.
Fazendo um procedimento semelhante utilizando o esquema para trás, obtém-se:
d2y
dx2
∣∣∣x=xi
=yi − 2yi−1 + yi−2
(∆x)2
sendo este o esquema para trás aplicado para a derivada segunda, sendo também um método
de primeira ordem.
Utilizando o esquema central, obtém-se a seguinte aproximação:
d2y
dx2
∣∣∣x=xi
=yi+1 − 2yi + yi−1
(∆x)2
sendo esta uma aproximação de segunda ordem. Novamente, a estratégia que normalmente
resulta em um melhor resultado é utilizar o esquema central para os pontos internos e os
esquemas para frente e para trás para o primeiro e para o último ponto, respectivamente.
3.2.3. Discretização das Condições de Contorno
Além da equação diferencial envolvendo a derivada da função, a resolução de Problemas
de Valor de Contorno necessita que determinadas condições de contorno também sejam em
determinados pontos do domínio.
Para o caso de condições que consistem em de�nir o valor da função em um dado ponto,
basta atribuir este valor para a função discretizada no ponto. Por exemplo, considere o
exemplo anterior de transferência de calor em uma barra, onde a temperatura nas duas
extremidades era conhecida: T = T1 em x = 0 e T = T2 em x = L. Considerando que
54
x0 corresponde ao ponto inicial x = 0 e xN corresponde ao ponto x = L, as condições de
contorno são implementadas de�nindo:
T0 = T1 TN = T2
Em muitos casos, a condição de contorno envolve a própria derivada de função em algum
ponto. Por exemplo, considere que no exemplo da transferência de calor na barra metálica, a
extremidade x = L fosse mantida isolada termicamente. Neste caso, a condição de contorno
é expressa como:dT
dx
∣∣∣x=L
= 0
Neste caso, é necessário discretizar a condição de contorno para transforma-la numa relação
algébrica. Como este ponto corresponde ao último ponto do domínio, os esquemas para
frente e central não pode ser utilizados, pois iriam depender de um ponto yN+1 que está fora
do domínio de solução. Assim, é necessário utilizar o esquema para trás, de modo que a
condição pode ser discretizada como:
dT
dx
∣∣∣x=L
=TN − TN−1
xN − xN−1
= 0 → TN = TN−1
ou seja, a condição de derivada nula na posição x = L implica que a variável em no ponto
xN é igual a variável no ponto anterior (xN−1).
3.2.4. Discretização de PVC's utilizando Diferenças Finitas
A partir das expressões obtidas para aproximar a derivada em um ponto através de ex-
pressões algébricas, pode-se transformar um problema de valor de contorno em um conjunto
de equações algébricas, sendo que para cada ponto do domínio discreto será atribuída uma
equação. Esta é uma característica importante dos métodos de diferenças �nitas: para que a
resolução seja possível, deve-se atribuir uma equação para cada ponto do domínio discreto.
Considere novamente a equação utilizada anteriormente para modelar a transferência de
calor em um barra metálica:d2T
dx2+ h′(Ta − T ) = 0
Para ilustrar os diferentes tipos de condição de contorno, considere agora que na fronteira
x = 0 a barra seja mantida a uma temperatura T = Text e na extremidade x = L a barra
esteja isolada, de modo que as condições de contorno neste caso serão:
T (0) = TextdT
dx
∣∣∣x=L
= 0
55
onde Text, h′, Ta e L são constantes. Para discretizar a equação, deve-se primeiramente
de�nir o grid numérico, ou seja, deve-se de�nir em quantos pontos o domínio de solução
contínuo será dividido e como estes pontos estão distribuídos. Neste caso, será assumido que
o número de pontos será N +1 = 6 e, por simplicidade, será considerado que os pontos estão
igualmente espaçados, portanto o domínio discreto será um conjunto de 6 pontos igualmente
espaçados, como ilustrado na �gura a seguir.
Como pode ser observado, quando 6 pontos são utilizados, o domínio contínuo é dividido
em 5 subdomínios com tamanho ∆x, ou seja, ∆x = L/5.
O objetivo do método de diferenças �nitas é obter uma aproximação para a temperatura
em cada um dos pontos xi, ou seja, deve encontrar 6 valores T0, T1, T2, T3, T4 e T5. Para isso,
é necessário que cada ponto possua uma equação algébrica associada. O primeiro passo para
a resolução é discretizar a equação. Neste caso, a equação possui somente uma derivada
segunda que deve ser discretizada. Para isso, será utilizar o esquema de diferenças central:
d2T
dx2=Ti+1 − 2Ti + Ti−1
∆x2
Além da derivada segunda, a equação também envolve o termo h′(Ta − T ). Como h′ e
Ta são constantes, basta utilizar os seus valores na equação discretizada. Como a expressão
anterior é utilizara para avaliar a derivada segunda no ponto xi, a temperatura T deve ser
substituída pela sua equivalente discreta Ti. Assim, a equação discretizada será:
Ti+1 − 2Ti + Ti−1
∆x2+ h′(Ta − Ti) = 0
Multiplicando por ∆x2 e agrupando os termos:
Ti+1 + Ti−1 − (2 + h′∆x2)Ti + ∆x2h′Ta = 0
56
ou ainda:
(2 + h′∆x2)Ti − Ti+1 − Ti−1 = ∆x2h′Ta
Esta relação pode ser utilizada para obter equações para a temperatura em todos os pontos
internos, ou seja, para todos os pontos exceto x0 e x5. Assim, para o ponto x1 temos que:
(2 + h′∆x2)T1 − T2 − T0 = ∆x2h′Ta
De forma semelhante, para o ponto x2:
(2 + h′∆x2)T2 − T3 − T1 = ∆x2h′Ta
e para os pontos x3 e x4:
(2 + h′∆x2)T3 − T4 − T2 = ∆x2h′Ta
(2 + h′∆x2)T4 − T5 − T3 = ∆x2h′Ta
Para estes pontos x0 e x5, deve-se utilizar as condições de contorno. Considerando a
primeira condição T (0) = Text, isto resulta em:
T0 = Text
A segunda condição de contorno é dada em termos de derivada nula. Neste caso, é preciso
discretizar a condição. Como comentado anteriormente, neste caso somente uma aproxima-
ção para trás pode ser utilizada, de modo que:
dT
dx
∣∣∣x=L
= 0 → T5 − T4
∆x= 0 → T5 − T4 = 0
Isto forma um conjunto de 6 equações lineares que podem ser utilizadas para a obtenção das
6 variáveis T0, T1, T2, . . . Na forma matricial, estas equações podem ser expressas como:
1 0 0 0 0 0
−1 (2 + h′∆x2) −1 0 0 0
0 −1 (2 + h′∆x2) −1 0 0
0 0 −1 (2 + h′∆x2) −1 0
0 0 0 −1 (2 + h′∆x2) −1
0 0 0 0 −1 1
T0
T1
T2
T3
T4
T5
=
Text
∆x2h′Ta
∆x2h′Ta
∆x2h′Ta
∆x2h′Ta
0
Assim, obtém-se um sistema linear tridiagonal que pode ser resolvido para obter a tempe-
ratura em cada um dos pontos do domínio discreto. Para ilustrar, considere um caso onde
57
L = 1 cm (de modo que ∆x = L/5 = 0.2 cm), Text = 100◦C, Ta = 25◦C e h′ = 0.1 cm−2.
Com isso, o sistema linear pode ser avaliado como:
1 0 0 0 0 0
−1 2.004 −1 0 0 0
0 −1 2.004 −1 0 0
0 0 −1 2.004 −1 0
0 0 0 −1 2.004 −1
0 0 0 0 −1 1
T0
T1
T2
T3
T4
T5
=
100
0.1
0.1
0.1
0.1
0
Para resolver este sistema linear, pode-se utilizar o algoritmo de Thomas (TDMA), como
visto em aulas anteriores. Relembrando, este método irá transforma a matriz dos coe�cientes
em uma matriz triangular superior. Os elementos abaixo da diagonal principal serão zerados
e os acima da diagonal principal não serão afetados. Os elementos da diagonal principal (a
partir da linha 2) são reavaliados como:
a′i,i = ai,i − (ai,i−1/ai−1,i−1)ai−1,i
De forma semelhante, os termos do lado direito são reavaliados como:
b′i = bi − (ai,i−1/ai−1,i−1)bi−1
Assim, o sistema pode ser reescrito como:
1 0 0 0 0 0
0 2.004 −1 0 0 0
0 0 1.505 −1 0 0
0 0 0 1.33955 −1 0
0 0 0 0 1.2575 −1
0 0 0 0 0 0.204759
T0
T1
T2
T3
T4
T5
=
100
100.1
50.0501
33.35
25.0008
19.8814
58
Resolvendo o sistema, obtém-se:
T [0] = T0 = 100◦C
T [0.2] = T1 = 98.83◦C
T [0.4] = T2 = 97.96◦C
T [0.6] = T3 = 97.38◦C
T [0.8] = T4 = 97.09◦C
T [1] = T5 = 97.09◦C
Para comparação, a solução exata em cada um destes pontos é:
T [0] = T0 = 100◦C
T [0.2] = T1 = 98.697◦C
T [0.4] = T2 = 97.689◦C
T [0.6] = T3 = 96.97◦C
T [0.8] = T4 = 96.54◦C
T [1] = T5 = 96.4◦C
Na �gura a seguir é apresenta uma comparação entre a solução exata (linha) e a solução
aproximada obtida com o método de diferenças �nitas (pontos). Pode-se observar que o
desvio é relativamente alto, sendo que um resultado melhor pode ser obtido aumentado-se o
número de pontos.
Neste exemplo, a resolução do problema envolveu um sistema linear tridiagonal. Quando
o método de diferenças �nitas é aplicado a um PVC linear, este sempre será o caso. Quando
59
aplicado em equações não-lineares, o sistema de equações algébricas obtido também será não-
linear e deverá ser resolvido com métodos adequados (método de Newton, por exemplo).
60
4. Método das Linhas
O método das linhas é um método semi-discreto para a resolução de EDP's que consiste em
discretizar as variáveis espaciais e manter uma das varáveis contínua (usualmente o tempo),
de modo a transformar a EDP em um sistema de EDO's que pode então ser resolvido através
dos métodos vistos anteriormente para a resolução de PVI's (como os métodos de Runge-
Kutta). A abordagem utilizada para a discretização das variáveis espaciais usualmente é o
método de diferenças �nitas, por isso o método das linhas é muitas vezes chamado de método
de diferenças �nitas semi-discreto.
Este método é aplicado principalmente para equações parabólicas, pois sua aplicação em
equações elípticas origina um conjunto de PVC's, o que por sua vez também precisam ser
resolvido por métodos de discretização. Quando aplicado em equações parabólicas, a va-
riável que possui um caminho característico é mantida contínua enquanto as demais são
discretizadas.
Por exemplo, considere a equação do calor:
∂T
∂t= α
∂2T
∂x2
Para obter uma solução particular para esta equação, é preciso especi�car duas condições de
contorno e uma condição inicial. Considere, por exemplo, as seguintes condições:
T (0, t) = Ta T (L, t) = Tb T (x, 0) = sin(x)
Neste caso, pode-se discretizar a derivada em relação à direção x. Usando um esquema
central:Ti+1 − 2Ti + Ti−1
∆x2
Substituindo esta forma discreta na EDP, obtém-se um sistema de EDO's para avaliar a
variação temporal das variáveis Ti:
dTidt
=α
∆x2(Ti+1 − 2Ti + Ti−1)
61
A aplicação das condições de contorno vão resultar em valores especí�cos para a variável
Ti nas extremidades x = 0 e x = L que serão válidos para qualquer tempo, visto que estas
condições são �xas. A condição T (0, t) = Ta vai resultar em T0 = Ta, enquanto que a
condição T (L, t) = Tb vai resultar em TN = Tb, onde N + 1 é o número total de pontos
utilizados para discretizar o domínio de solução na direção x.
Para resolver este sistema de EDO's para os i pontos, é preciso especi�car condições iniciais
para cada valor Ti. Como a variável Ti representa a temperatura na posição xi, o valor inicial
pode ser obtido diretamente da condição inicial especi�cada anteriormente aplicada no ponto
xi. A condição inicial era da forma:
T (x, 0) = sin(x)
Assim, para cada variável Ti teremos uma condição inicial associada da forma:
Ti(0) = sin(xi)
Observe que enquanto a temperatura é uma função da posição e do tempo (T (x, t)), as
variáveis Ti são funções apenas do tempo Ti(t), pois representam a temperatura em um
ponto �xo xi. A aplicação do método das linhas vai originar uma série de curvas (daí o nome
método das linhas), contínuas em relação ao tempo, que representam como a temperatura
varia em cada ponto xi. A �gura a seguir ilustra a forma da solução de uma EDP com duas
variáveis independentes obtida com o método das linhas.
62
Para ilustrar a aplicação do método das linhas, considere que se deseje obter a variação de
temperatura ao longo de uma barra metálica com uma extremidade isolada e outra mantida
a uma temperatura Text e que perde calor para o meio externo por convecção, da mesma
forma que analisado na aula anterior. Neste caso, porém, considere que se deseje obter como
a temperatura varia ao longo do tempo a partir de um estado inicial T (x, 0) = Tini.
A equação que descreve a variação na temperatura ao longo da posição x e do tempo t
neste caso será:∂T
∂t= α
∂2T
∂x2− h′α(T − Ta)
As condições de contorno associadas a este problema são:
T (0, t) = Text∂T
∂x
∣∣∣x=L
= 0
Além disso, a condição inicial utilizada é da forma:
T (x, 0) = Tini
De forma geral, as condições de contorno podem ser função do tempo, da mesma forma que
a condição inicial pode ser uma função de x.
Para resolver esta equação com o método das linhas, é preciso discretizar a equação na
direção espacial e manter a função contínua no tempo. Assim, deve-se de�nir um domínio
discreto na direção x. Neste caso, será considerado o mesmo domínio discreto utilizada
anteriormente, com N + 1 = 6 pontos:
Neste caso, a discretização da EDP na direção x irá resultar um conjunto de valores
(T0, T1, T2, T3, T4, T5) que representam a temperatura nos pontos respectivos (x0, x1, x2, x3, x4, x5).
Neste caso, porém, estes valores Ti não são necessariamente constantes, mas são uma função
do tempo.
A EDP avaliada neste exemplo envolve a derivada segunda em relação a direção x, então
deve-se discretizar esta derivada. Considerando um esquema central, a derivada segunda nos
63
pontos xi são dadas por:∂2T
∂x2=Ti+1 − 2Ti + Ti−1
(∆x)2
Substituindo na EDP, obtém-se:
dTidt
= α
(Ti+1 − 2Ti + Ti−1
(∆x)2
)− h′α(Ti − Ta)
Esta equação pode ser aplicada nos pontos i = 1, 2, 3, 4 para obter uma EDO para a tem-
peratura em cada um destes pontos. Para resolver este sistema de EDO's, é preciso de�nir
uma condição inicial para cada ponto. Como a temperatura inicial é considerada constante
e igual a Tini, basta associar este valor com a temperatura em cada ponto:
Ti(0) = Tini
Para fechar o sistema de equações, é preciso ainda de�nir equações para T0 e T5, que são
obtidas através da aplicação das condições de contorno. A condição T (0, t) = Text resulta
diretamente em:
T0 = Text
Assim, obtém-se uma equação para a temperatura na extremidade x = 0. Na extremidade
x = L, a condição é de derivada nula. Como visto anteriormente, neste caso pode-se aplicar
um esquema de discretização para trás, de onde se obtém que T5 = T4. Com estas duas
condições extras, o sistema de EDO's pode ser resolvido.
64
De forma resumida, as equações para os 6 pontos são:
T0 = Text
dT1
dt= α
(T2 − 2T1 + T0
(∆x)2
)− h′α(T1 − Ta) T1(0) = Tini
dT2
dt= α
(T3 − 2T2 + T1
(∆x)2
)− h′α(T2 − Ta) T1(0) = Tini
dT3
dt= α
(T4 − 2T3 + T2
(∆x)2
)− h′α(T3 − Ta) T1(0) = Tini
dT4
dt= α
(T5 − 2T5 + T3
(∆x)2
)− h′α(T4 − Ta) T1(0) = Tini
T5 = T4
As equações para i = 1, 2, 3, 4 formam um sistema de PVI's que pode ser resolvido por
algum dos métodos vistos anteriormente para a resolução de PVI's, como os métodos de
Runge-Kutta.
Considere, por exemplo, novamente que um caso onde L = 1 cm (∆x = 0.2 cm), Text =
100◦C, Ta = 25◦C e h′ = 0.1 cm−2. Além disso, considere que α = 0.01 cm2/s e Tini = 25◦C.
Neste caso, a barra metálica está inicialmente na mesma temperatura que o ambiente externo.
Em um dado instante, a temperatura na extremidade x = 0 é aumentada para 100◦C.
Resolvendo o sistema de equações utilizando Runge-Kutta de quarta ordem, obtém-se as
curvas apresentadas na �gura a seguir.
65
Pode-se observar que para altos valores de tempo as temperaturas tendem a um valor es-
pecí�co, ou seja, tendem para um valor estacionário. Estes valores correspondem exatamente
aos obtidos na aula anterior para a resolução do caso onde o comportamento transiente foi
desprezado.
66
