Métodos Numéricos na resolução de Equações Diferenciais Ordinárias (EDO)

KELTON JORGE GARCIA SANTOS

Métodos Numéricos na Resolução de PVI’s

Aplicações com MATLAB

Licenciatura em ensino de Matemática

UNICV, 2009

ii

KELTON JORGE GARCIA SANTOS

Métodos Numéricos na Resolução de PVI’s

Aplicações com MATLAB

Trabalho científico apresentado na UNICV

para a obtenção do grau de Licenciado em

Ensino de Matemática sob a orientação do

Engenheiro, Aurélio Vicente.

Licenciatura em ensino de Matemática

UNICV, 2009

iii

O Júri

O Orientador

O Arguente

O Presidente

UNICV, Praia aos ---------/----------/2009

iv

Dedico este trabalho à minha família que sempre

me apoiou e me incentivou durante todo esse

tempo.

v

Agradecimentos

Minha imensa gratidão a Deus por ter conseguido chegar até essa etapa.

Agradeço ao meu orientador Aurélio Vicente pela confiança, paciência e

disponibilidade que demonstrou para comigo durante a realização desse trabalho.

Finalmente, agradeço a todos os meus colegas e professores que de forma directa ou

indirectamente contribuíram para que eu chegasse até essa etapa.

vi

ÍNDICE

INTRODUÇÃO/CONTEXTUALIZAÇÃO ........................................................................................ 8

1. Nota Introdutória .......................................................................................................................................... 8

2. Objectivos ......................................................................................................................................................... 9

3. Importância do Estudo ............................................................................................................................... 9

4. Metodologias ................................................................................................................................................ 10

EQUAÇÕES DIFERENCIAIS ORDINÁRIAS ................................................................................ 11

1.1 Introdução .................................................................................................................................................. 11

1.2 Solução de uma EDO .............................................................................................................................. 13

1.2.1 Solução Geral e Solução particular duma EDO ............................................................................... 14

1.3 Métodos Analíticos para a resolução de EDO de 1ª ordem................................................. 15

1.3.1 Equações Diferenciais de Variáveis Separáveis ............................................................................. 15

1.3.2 Diferencial Exacta ...................................................................................................................................... 16

1.3.3 Equações diferenciais Lineares de 1ª ordem.................................................................................. 17

1.4 Problema de Valores Iniciais e Problemas de Valores na Fronteira .............................. 17

1.4.1 Interpretação Geométrica ...................................................................................................................... 18

1.4.2 Existência e Unicidade de Solução de PVI’s ..................................................................................... 20

1.5 Soluções Numéricas de PVI’s ............................................................................................................. 21

MÉTODOS DE PASSO SIMPLES ................................................................................................... 23

2.1 Métodos da série de Taylor ................................................................................................................ 24

2.1.1 Metodos de Euler ................................................................................................................................. 25

2.1.1.1 Método de Euler Explícito (ou Progressivo) .......................................................................... 25

2.1.1.2 Metodo de Euler Implícito (ou Regresivo) .............................................................................. 27

2.1.1.3 Erros dos métodos de Euler .......................................................................................................... 29

2.1.2 Método de Heun ................................................................................................................................... 30

2.1.2.1 Erros do método de Heun ................................................................................................................... 31

2.1.3 Métodos de Taylor de ordem mais elevada ............................................................................ 31

2.1.3.1 Erros ........................................................................................................................................................ 32

2.1.4 Análise da Consistência dos métodos da série de Taylor ................................................ 32

2.2 Métodos de Runge Kutta (RK) .......................................................................................................... 34

2.2.1 Métodos de Runge Kutta de primeira ordem ............................................................................ 37

2.2.2 Métodos de Runge Kutta de terceira ordem .............................................................................. 37

Métodos Numéricos na Resolução de PVI’s. 7

2.2.3 Método de Runge Kutta de 4ª ordem ............................................................................................ 38

2.2.4 Erros dos métodos de Runge Kutta .................................................................................................... 39

2.3 Análise da Convergência dos métodos de passo simples .................................................... 41

MÉTODOS DE PASSO MÚLTIPLO .............................................................................................. 43

3.1 Métodos de Adams Bashforth (AB) ................................................................................................ 45

3.2 Métodos de Adams Moulton (AM) .................................................................................................. 48

3.3 Processo preditor corrector............................................................................................................... 51

3.3.1 Considerações .............................................................................................................................................. 51

CONCLUSÃO ..................................................................................................................................... 55

Referências Bibliográficas ......................................................................................................... 57

Bibliografia ...................................................................................................................................... 58

ANEXO ................................................................................................................................................ 59

vii


INTRODUÇÃO/CONTEXTUALIZAÇÃO

1. Nota Introdutória

Vários fenómenos nos campos da Física, Engenharia, Economia, Biologia etc., são

modelados e resolvidos matematicamente por meio das equações diferenciais.

Existem vários métodos analíticos para solucionar Equações Diferenciais Ordinárias

(EDO’s), mas esses métodos não são aplicáveis para todos os problemas, daí a

necessidade de recorrer a técnicas numéricas para obter soluções aproximadas para tais

problemas.

É nessa óptica que vamos trabalhar o tema acima referido, onde vamos apresentar

os principais métodos numéricos utilizados para solucionar (aproximadamente) EDO’s

de primeira ordem com valores iniciais (PVI’s).

Alguns dos métodos apresentados serão acompanhados por rotinas desenvolvidas

em MATLAB, que para além de ser um software que utiliza uma linguagem de

programação muito acessível, possibilita o estudo e análise dos resultados obtidos pelos

vários métodos, bem como ajuda a minimizar os erros cometidos nas aproximações


(pois no MATLAB temos a oportunidade de utilizar números com várias casas decimais,

o que ajuda e diminuir os erros de arredondamento).

2. Objectivos

Pretendemos com a realização deste trabalho:

Apresentar técnicas numéricas para a resolução de PVI’s;

Analisar os erros de aproximação;

Fazer Confrontação entre os vários métodos, identificando as vantagens e

desvantagens de cada um;

Desenvolver algoritmos em MATLAB para alguns métodos apresentados.

3. Importância do Estudo

Ao estudar as EDO’s, observamos que uma grande parte delas não podem ser

resolvidas pelos métodos analíticos conhecidos, e algumas que conseguimos resolver

apresentam soluções extremamente complexas.

Deste modo, na ausência de soluções analíticas, o problema é solucionado,

recorrendo a soluções aproximadas, obtidas numericamente.

Essas técnicas apresentam grande potencial de utilização para a disciplina de

Análise Numérica, podendo ser utilizadas por estudantes de várias áreas tais como:

Engenharia, Física, Biologia, Química etc.


4. Metodologias

Para um trabalho deste género como é normal, a metodologia utilizada foi

basicamente de pesquisas bibliográficas, mas também a internet foi um auxílio

importante para a realização do referido trabalho.

O software MATLAB serviu de suporte não só para a implementação e execução dos

algoritmos, mas também para a comparação dos resultados obtidos pelos diferentes

métodos (os que foram implementados).


CAPÍTULO I

EQUAÇÕES DIFERENCIAIS ORDINÁRIAS

1.1 Introdução

A origem do estudo das Equações Diferenciais coincide com a do Cálculo Diferencial

e Integral introduzidas no séc. XVII por Isaac Newton e Gottfried Leibnitz. As Equações

Diferenciais tais como as que conhecemos hoje foram introduzidas posteriormente por

Gottfried Leibnitz.

Do séc. XVIII até o séc. XX as Equações Diferenciais eram utilizadas Somente no

estudo de fenómenos Físicos, mas precisamente na Mecânica Clássica.

Actualmente existem várias áreas tais como a Física, Biologia, Ecologia, Química,

Economia etc., que fazem uso das equações Diferenciais como meio para solucionar

vários problemas inerentes às referidas áreas.

Problemas relacionados com a Dinâmica populacional (desde seres humanos até

bactérias), fluxo sanguíneo no sistema vascular, volumes de líquidos celulares,


desintegração de substâncias radioactivas, elementos da electricidade (Intensidade da

corrente eléctrica, capacitância, indutância etc.), estudos epidemiológicos etc., são alguns

dos inúmeros problemas que podem ser modelados e resolvidos matematicamente por

meio das Equações Diferenciais.

Em termos conceptuais, uma Equação Diferencial é uma relação entre uma função e

suas derivadas. Entretanto, existem dois tipos de Equações Diferenciais:

Equação Diferencial Ordinária (EDO). Equação que envolve uma só função

(incógnita) -normalmente designada por y - de uma só variável (normalmente

designada por x ou por t ), e um número finito de derivadas dessa função.

(AGUDO, 1992)

Equação Diferencial Parcial (EDP). Equação que contém funções com mais do

que uma variável e derivadas parciais dessas funções em relação a essas

variáveis. (AGUDO, 1992)

Para esse Trabalho vamos estudar somente as Equações Diferenciais Ordinárias. Para

mais informações sobre as EDP, consultar qualquer manual que trata da teoria das

EDP’s, como por exemplo [AGU92].

Definição 1.1 Seja 2:

nF IR IR

+Ω ⊆ → , uma função real definida num aberto Ω . À

equação

( )( ), , , ´ ,..., 0n

F x y y y y = (1.1)

com y função real de variável x (definida em certo intervalo I IR⊂ ) dá-se o nome de

Equação Diferencial Ordinária (EDO) de ordem n . (AGUDO, 1992)

Pode-se dizer que se trata de uma equação que envolve uma só função (que denotamos

por y ) de uma só variável (que denotamos por x ) e um número finito de derivadas

dessa função.

Exemplo 1.1 ´ 4 0y y+ =


Exemplo 1.2 ´ 3 ´ xy yy e+ =

Exemplo 1.3 8 3 0dy

xdx

+ − =

Exemplo 1.4 2

1 1

4

dy

y dx x x=

−

Definição 1.2 A ordem de uma EDO é a ordem da derivada de maior ordem da variável

dependente que aparece na EDO.

Nesse Trabalho a nossa atenção vai estar centrada unicamente para as equações

diferenciais ordinárias de 1ª ordem.

Definição 1.3 O grau de uma EDO é o expoente da derivada de maior ordem da

variável dependente que aparece na EDO.

Exemplo 1.5 ( )9 2´ý xy x− = , Tem Ordem 9 e Grau 1;

Exemplo 1.6 ( ) ( )3 10

´ ´ 6y y y x+ + = , Tem Ordem 2 e Grau 3;

Exemplo 1.7 2

2

23 2

d y dyy x

dx dx− + = , Tem Ordem 2 e Grau 1;

Exemplo 1.8 4 52

2

20

d y dyx

dx dx

+ − =

, Tem Ordem 2 e Grau 4;

1.2 Solução de uma EDO

Definição 1.4 Chama-se solução da equação diferencial (1.1) no intervalo I IR⊂ a toda

a função ( ) :y y x I IR IR= ⊂ → com derivadas até a ordem n tal que

( ) ( ) ( ) ( ) ( ), , ´ , ´ ,..., 0, n

F x y x y x y x y x x I = ∀ ∈ (1.2)

Por outras palavras, a solução da EDO (1.1) é uma função ( )y y x= , com derivadas até a

ordem n que satisfaça (1.1). (AGUDO, 1992)


Exemplo 1.9 ( ) 2y x x c= + é solução da EDO ´ 2y x=

Exemplo 1.10 ( ) ( )y x sen x= é solução da EDO ´ 0y y+ =

Exemplo 1.11 ( ) 2 1y x x= − não é solução da EDO ( )4 2

´ 1y y+ = −

Definição 1.5 Uma equação diferencial diz-se na forma normal quando for

“explicitada” em relação à derivada de maior ordem, i.e., quando for escrita na forma

( ) ( )( )1, , , ´ ,...,

n ny f x y y y y

−= (1.3)

Por exemplo uma EDO de 1ª ordem na forma normal é dada por ( )´ ,y f x y=

1.2.1 Solução Geral e Solução particular duma EDO

Definição 1.6 O conjunto de todas as soluções de uma equação diferencial (num dado

intervalo diz-se a solução geral (ou integral geral) da equação (nesse intervalo).

Resolver ou integrar uma equação diferencial é determinar o seu integral geral.

Exemplo 1.123

1 26

xy c x c= − + + ( 1c e 2c constantes reais arbitrários) é Solução Geral da

EDO ý x= − .

Exemplo 1.13 xy ce

−= é Solução Geral da EDO ´ 0y y+ = .

Definição 1.7 Chama-se Solução Particular ou Integral Particular duma equação

diferencial a toda a solução obtida atribuindo valores às constantes arbitrárias da

Solução Geral.

Exemplo 1.14 2 xy e−= é uma Solução Particular da EDO ´ 0y y+ = .

Exemplo 1.15 3

2 36

xy x= − + + é uma Solução Particular da EDO ý x= − .


Deve-se dizer que nem sempre é possível determinar a solução geral duma equação

diferencial. Em muitos problemas basta-nos saber que existem soluções e que estas se

comportam de determinada forma.

1.3 Métodos Analíticos para a resolução de EDO de 1ª ordem

Apresentamos em seguida alguns tipos de EDO de 1ª ordem da forma ( ),′ =y f x y (ou

redutíveis a esta forma) e os respectivos métodos analíticos utilizados na sua resolução.

Em particular, uma equação

( ) ( ), , 0N x y y M x y′ + = (1.4)

será redutível à forma normal sempre que ( ), 0N x y =/ . Os pontos que anulam

( ),N x y merecem uma atenção especial (AGUDO, 1992). O leitor pode ler [AGU92] para

melhor esclarecimento.

1.3.1 Equações Diferenciais de Variáveis Separáveis

Definição 1.8 Uma EDO de 1ª ordem ( ),dy

f x ydx

= diz-se de variáveis separáveis se

poder tomar a forma:

( ) ( )h y y g x′ = (1.5)

que também se escreve

( ) ( )h y dy g x dx= ,

pois dy y dx′=


Método de resolução.

Para encontrar a solução geral desta equação diferencial, basta primitivar as funções de

ambos os membros, ou seja,

( ) ( )h y dy g x dx=∫ ∫ .

1.3.2 Diferencial Exacta

Definição 1.9 Uma equação diferencial ( ) ( ), , 0M x y dx N x y dy+ = , com M e N contínuas

nalgum aberto Ω de 2IR diz-se diferencial exacta, se satisfazer a condição: (RILEY,

HOBSON, & BENCE, 2007)

( ) ( ), ,M x y N x y

y x

∂ ∂=

∂ ∂ (1.6)

Método de Resolução

Para encontrar a solução da equação diferencial exacta ( ) ( ), , 0M x y dx N x y dy+ = ,

temos de encontrar uma função ( ),F x y de classe 1C em Ω tal que:

( )

( )( )

( )∂ ∂

= =∂ ∂

e F x y F x y

M x y N x yx y

, ,, , (1.7)

A solução da equação diferencial exacta será dada por:

( ),F x y C=

Na maioria das vezes a equação diferencial ( ) ( ), , 0M x y dx N x y dy+ = não é exacta, mas

normalmente é possível transformá-lo numa equação diferencial exacta por meio da sua

multiplicação por um factor integrante ( ),x yµ , tal que a equação

( ) ( ) ( ) ( ), , , , 0x y M x y dx x y N x y dyµ µ+ = seja uma equação diferencial exacta, isto é,

( ) ( ) ( ) ( ), , , ,x y M x y x y N x y

y x

µ µ∂ ∂=

∂ ∂ .


1.3.3 Equações diferenciais Lineares de 1ª ordem

Definição 1.10 Uma equação diferencial linear de 1ª ordem é uma equação da forma

( ) ( )ý p x y q x+ = , com ( ) ( ) e p x q x contínuas nalgum intervalo I IR⊂ .

• Se ( ) 0q x ≠ denomina-se equação completa.

• Se ( ) 0q x = , a equação é de variáveis separáveis. (RILEY, HOBSON, & BENCE,

2007)

Método de resolução

Para resolver equações desse tipo, temos de encontra o factor integrante ( )( )p x dx

x eµ ∫= e

multiplicá-lo pela equação ( ) ( )ý p x y q x+ = . Pelo que a solução geral dessa equação

será:

( ) ( )

( ) ( )

x q x dx Cy

x x= +∫ µ

µ µ.

É Importante salientar que a maioria das Equações diferenciais de 1ª ordem não se

enquadram em nenhuma das classificações acima referidas, nem podem ser

transformadas em nenhuma delas.

1.4 Problema de Valores Iniciais e Problemas de Valores na Fronteira

Como foi dito anteriormente nem sempre é possível encontrar a solução geral duma

equação diferencial. Na maior parte das aplicações não interessa conhecer a solução

geral (mesmo quando for possível determina-lo) mas sim soluções particulares que

satisfaçam determinadas condições fixadas previamente. Assim temos:

Definição 1.11 Chama-se Problema de Valores Iniciais (PVI) ou Problema de Cauchy

a todo o problema que consiste em determinar a solução de uma EDO de 1ª ordem

satisfazendo condições iniciais dadas num ponto (do Intervalo I IR⊂ considerado). O

problema é dado da seguinte forma:


( )( )

( )0 0

, , y f x y x x I

y x y

′ = ∈

= (1.8)

Exemplo 1.16 A função 2 ( constante)kxy e k= é solução do problema de valores

iniciais( )

´

0 2

y ky

y

=

=

Observação 1.1 Como acontece no exemplo acima, muitas vezes não se explicita, na

prática, o intervalo em que estamos a considerar a equação. Neste exemplo poderia ser

qualquer intervalo que contivesse a origem, nomeadamente ] [,I = −∞ +∞ .

Exemplo 1.17 A função 2sin( ) cos( )y x x= + é solução do problema de valores iniciais

( )

( )

´´ 0

0 1 .

´ 0 2

y y

y

y

+ =

=

=

De uma forma geral tem-se o seguinte

Se f é contínua com respeito a x , em 2IRΩ ⊂ , então a solução de (1.8) satisfaz

( ) ( )( )0

0,

x

xy x y f s y s ds= + ∫ (1.9)

Da mesma maneira, se y for definido por (1.9), então y é contínua em I e ( )0 0y x y=

(QUORTERONI, SACCO, & SALERI, 2007).

1.4.1 Interpretação Geométrica (AGUDO, 1992)

Uma equação da forma ( ),y f x y′ = , com f definida nalgum aberto 2IRΩ ⊂ , associa

a cada ponto ( ),x y de Ω uma recta de coeficiente angular ( ),y f x y′ = e, portanto, uma

direcção. Obtêm-se, assim, o que se chama um campo de direcções associada à equação

diferencial dada.


Ao gráfico de uma solução, ( )y y x= , da equação diferencial dá-se o nome de curva

integral ou curva de soluções e, é obvio que uma tal curva é em cada um dos pontos

( )( ),x y x tangente à recta de coeficiente angular ( ),y f x y′ = .

Esta propriedade permitirá obter, embora de forma aproximada, a curva integral de uma

equação diferencial ( ),y f x y′ = que passa por um ponto dado (ou seja a solução gráfica

de um PVI) substituindo, em pequenos intervalos, a curva solução por seguimentos

tangentes à curva. Para facilitar a construção podemos começar por desenhar algumas

linhas com equações da forma ( ), .tef x y c= , que são camadas por isoclínicas (ou seja

“linhas com a mesma inclinação”) da equação diferencial dada, pois, em todos os pontos

da linha ( ), .tef x y c= as rectas associadas a ( ),y f x y′ = têm um mesmo coeficiente

angular. A figura 1.1 apresenta as soluções particulares do exemplo 1.9 para

0, 1 2 e c c c= = =

Figura 1.1. Curvas soluções do exemplo 1.9 para c=0, c=1 e c=2

Definição 1.12 Chama-se Problema de Valores na Fronteira (PVF) ou Problemas de

Contorno a todo o problema que consiste em determinar a solução (ou soluções) de uma

equação diferencial satisfazendo condições iniciais em dois ou mais pontos do Intervalo

I IR⊂ considerado. O problema é dado pela seguinte forma:


( )( )( )

( )0 0

1 1

, ,

y f x y x x I

y x y

y x y

′ = ∈

=

=

(1.10)

Exemplo 1.18 Tendo em conta o exemplo (1.12), a função 3 1

6 6

xy = − + é a solução do

problema de contorno ( ) ( )

( ) ( )

´ 0

10 ´ 1

3

1 ´ 0 0

y x

y y

y y

+ =

+ = − + =

(onde 1 0c = e 2

1

6c = ).

A partir daqui vamos tratar unicamente de PVI. Se o leitor estiver interessado em

aprofundar o seu conhecimento sobre PVF deve procurar qualquer livro que trata das

equações diferenciais.

Antes de tentar resolver um PVI, precisamos de nos assegurar de que esta tem

solução e que esta solução é única. Para isso o teorema seguinte relembra-nos as

condições suficientes para que tal se verifique.

1.4.2 Existência e Unicidade de Solução de PVI’s

Teorema 1.1 (de existência e unicidade)

Seja 2:f IR IRΩ ⊆ → uma função satisfazendo as seguintes condições:

f é contínua e ( )0 0,x y um ponto de Ω ( aberto).

f é Lipschitziana relativamente ao seu segundo argumento, i.e., existe uma

constante L com 0 L≤ ≤ ∞ (a constante de Lipchitz), tal que

( ) ( )2 1 2 1 1, 2, , , , f x y f x y L y y x y y IR− ≤ − ∀ ∈ Ω ∀ ∈

Então existe um intervalo [ ]0 0, , 0I x h x h h= − + > , onde o problema de Cauchy (1.8), tem

uma e uma só solução. (AGUDO, 1992)


Demonstração. Ver (AGUDO, 1992).

Para todos os métodos numéricos a serem apresentados nos capítulos posteriores,

consideraremos o PVI (1.8) com a hipótese que f satisfaz as condições do teorema 1.1.

1.5 Soluções Numéricas de PVI’s

Para algumas equações diferenciais existem métodos analíticos para a obtenção da

solução exacta. No entanto quando a função for muito “complexa”, ou for conhecida

apenas através de valores tabelados esses métodos se mostram ineficazes ou então

apresentam soluções muito complexas, dificultando a análise dos mesmos. Daí a

necessidade de recorrer aos métodos numéricos que oferecem a possibilidade de

encontrar não as soluções exactas, mas sim soluções aproximadas (em que as

aproximações variam de método para método) desses problemas.

Com o advento dos computadores, o desenvolvimento de técnicas computacionais,

particularmente a criação de softwares tais como o MAPLE, MATLAB e o MATHEMATICA

possibilitou trabalhar com os métodos numéricos com uma grande precisão e

cometendo uma percentagem de erro muito reduzido.

Nesse trabalho vamos utilizar o software MATLAB para a implementação dos vários

métodos a serem apresentados, pois trata-se de um software de fácil utilização e que é

muito utilizado no campo da engenharia. O software utiliza uma linguagem de

programação muito simples possibilitando ao programador a implementação rápida dos

algoritmos, além de trazer um grande número de funções pré implementados. (Ver

anexo para mais esclarecimentos sobre o software MATLAB)

Utilizando métodos numéricos na resolução de EDO temos a possibilidade de obter

soluções aproximadas num conjunto discreto de pontos nx do intervalo [ ]0 ,I x X= .

Desse modo, em vez da solução exacta, ficaremos a conhecer os valores aproximados

desta função num conjunto discretos de pontos nesse intervalo (diz-se que o problema

foi discretizado).


O processo de discretização pode ser apresentado da seguinte maneira:

(QUORTERONI, SACCO, & SALERI, 2007)

Consideremos [ ]0 ;I x X IR= ⊂ o intervalo de integração. Vamos considerar uma

sucessão de nós de discretização de I, ( )0 0,1, 2,......n

x x nh n N= + = , com 0h > , e onde N

é o máximo inteiro tal que Nx X= . Essa sucessão n

x define subintervalos [ ]1,n n n

I x x += .

Em geral fazemos 0X xh

N

−= ( N não é mais do que o número de valores da variável

independente onde queremos calcular os valores aproximados da solução). No processo

acima descrito, os nós de discretização formam uma malha uniforme do intervalo I que

normalmente são designados por nós da malha. À distância h dá-se o nome de

tamanho de passo da malha e N é o número de passos

Em termos de notações vamos representar a solução exacta no nó ( )( ),n n

x y x , por ( )ny x

e a solução numérica aproximada no mesmo nó será representada por n

y , isto é,

( )ny x (solução exacta),

( )n ny y x≅ (solução aproximada).

Também, considerarmos que ( ),n n n

f f x y≡ e obviamente ( )0oy y x= .

Nosso objectivo é então determinar aproximações ny da solução verdadeira ( )ny x nos

pontos da malha. Portanto a solução numérica será uma tabela de valores dos pares

( ),n nx y , onde como já foi dito anteriormente ( )n ny y x≅ .


CAPÍTULO II

MÉTODOS DE PASSO SIMPLES

Definição 2.1 Um método numérico para a aproximação do PVI (1.8) é chamado

Método de Passo Simples (ou passo único) se para 0n∀ ≥ , a solução aproximada 1n

y +

depender somente de ny . Caso contrário é designado método de passo múltiplo (ou

multipasso). (QUORTERONI, SACCO, & SALERI, 2007)

Por outras palavras, o método de passo simples é um método numérico para resolver

PVI’s em que apenas o resultado do passo anterior é utilizado para determinar o

resultado do passo seguinte.

Exemplo 2.1 Método de Euler Progressivo ( ou Explícito)

1n n ny y hf+ = + (2.1)

Exemplo 2.2 Método de Euler Regressivo (ou implícito)

1 1n n ny y hf+ += + (2.2)


Exemplo 2.3 Método de Crank-Nickolson

[ ]1 12

n n n n

hy y f f+ += + + (2.3)

Exemplo 2.4 Método de Heun

( )1 1,2

n n n n n n

hy y f f x y hf+ += + + + (2.4)

Definição 2.2 Um método numérico para a aproximação do PVI (1.8) é chamado método

explícito se 1ny + poder ser calculado directamente em termos de valor (ou valores) de

comk

y k n≤ . Um método diz-se implícito se 1ny + depender implicitamente de si

próprio através de f . (QUORTERONI, SACCO, & SALERI, 2007)

Os métodos dos exemplos (2.1) e (2.4) são métodos explícitos, enquanto os métodos dos

exemplos (2.2) e (2.3) são métodos Implícitos.

2.1 Métodos da série de Taylor

Os métodos da série de Taylor para a aproximação do PVI (1.8) são os métodos

obtidos usando o desenvolvimento em série de Taylor centrada em nx (i.e., admitindo

que ( ) [ ]( )1

0 ,p

y x C x X+∈ ) da solução do PVI (1.8),

( ) ( ) ( ) ( ) ( ) ( )

( )( ) ( ) [ ]

2

1

11

1 1 0

2! !

, , ,1 !

pp

n n n n n

pp

n n n n n n

h hy x y x hy x y x y x

p

hy x x x x x X

p

+

++

+ +

′ ′′= + + + + +

+ < < ∈+

K

ξ ξ

(2.5)

em que 1n nx x h+ = + .

Dentre os métodos da série de Taylor para a aproximação do PVI (1.8) temos:


2.1.1 Metodos de Euler

Os primeiros métodos a serem apresentados são os métodos de Euler. Por serem os

mais simples, são estudados com maior profundidade, o que facilita a abordagem aos

restantes métodos. Dentre os métodos de Euler temos:

2.1.1.1 Método de Euler Explícito (ou Progressivo)

Consideremos o PVI (1.8)

( )( )

( )0 0

, ,

y f x y x x I

y x y

′ = ∈

=

Pretendemos obter a aproximação de ( )y x no intevalo [ ]0 ,I x X= . Discretizando (1.8),

vamos determinar as soluções aproximadas nos pontos , 1, 2,......n

x n N= .

Recorendo ao teorema de Taylor obtemos:

( ) ( ) ( ) ( )2

1 12!

n n n n n n n

hy x y x hy x y x x+ +

′ ′′= + + < <ξ ξ (2.6)

Atendendo a que 1n nh x x+= − e ( )( ) ( ),

n n nf x y x y x′= . Fazendo 0n = e desprezando os

termos de ordem 2≥ em (2.6)* obtemos a solução aproximada

( )1 0 0 0,y y hf x y= + (2.7)

Seja 1y a solução aproximada de ( )1y x , substituido ( )1y x por 1

y em (2.6) para 1n = e

procedendo de maneira análoga obtemos para ( )2y x a aproximação

( )2 1 1 1,y y hf x y= + (2.8)

Desta forma conseguimos gerar uma aproximação numérica para o PVI (1.8) nos pontos

nx fazendo

* Os termos de ordem 2≥ em (2.6) normalmente são representados por ( )2

O h , que corresponde aos

termos truncados em (2.6)).


( )1 , 0,1, 2,..., 1n n n n

y y hf x y n N+ = + = − (2.9)

O método de aproximação numérica (2.9) para obter a solução aproximada ny da

solução exacta ( )ny x é conhecido como método de Euler Explícito (ou Progressivo).

(RAINER, 1998)

Exemplo 2.5 Utilizando o método de Euler progressivo encontrar a solução do PVI

[ ]2

; 0,3(0) 1

y x yI

y

′ = −=

=. Compare soluções para

1 1 11, , e

2 4 8h = .

A solução exacta do PVI acima é dado por ( ) 2 2 2xy x e x x−= − + − + . Para 1

4h = , utilizando

o método de Euler progressivo (2.9) teremos

( )

( )1

2

1.0 0.25 0 1 0.75

0.75 0.25 0.0625 0.75 0.5781, etc...

y

y

= + − =

= + − =

Essa iteração continua até atingir o último passo (que corresponde a 3 0

120.25

N−

= = )

( )2

12(3) 3.7808 0.25 2.75 3.7808 4.7262y y≈ = + − =

Na tabela 2.1 estão representadas as soluções aproximadas pelo método de Euler

Progressivo (que corresponde ao PROGRAMA 1 do anexo) e a solução exacta do PVI

dado.

Avaliando a tabela 2.1 vê-se que tomando o tamanho do passo ( h ) cada vez menor, a

solução aproximada se aproxima cada vez mais da solução exacta, i.e., quanto menor o

valor do tamanho do passo melhor é a solução aproximada obtida. Por conseguinte, a

diminuição do tamanho do passo tem como desvantagem o aumento de iterações, i.e.,

diminuindo o tamanho do passo, aumenta-se o número de nós no intervalo dado, onde

temos de calcular a solução aproximada.


nx

ny ( ) Exacta

ny x

= 1h .= 0 5h .= 0 25h .= 0 125h

.

.

.

.

.

.

.

.

.

.

0

0 125

0 25

0 375

0 50

0 75

1 00

1 50

2 00

2 50

3 00

.

.

.

.

1 0000

0 0000

1 0000

4 0000

.

.

.

.

.

.

.

1 0000

0 5000

0 3750

0 6875

1 4688

2 7344

4 4922

.

.

.

.

.

.

.

.

.

1 0000

0 7500

0 5781

0 4961

0 5127

0 8665

1 6749

2 9578

4 7262

.

.

.

.

.

.

.

.

.

.

.

1 0000

0 8750

0 7676

0 6794

0 6121

0 5448

0 5743

0 9488

1 7717

3 0644

4 8395

.

.

.

.

.

.

.

.

.

.

.

1 0000

0 8831

0 7836

0 7033

0 6434

0 5901

0 6321

1 0268

1 8446

3 1679

4 9502

Tabela 2.1 Resultados do exemplo 2.5

2.1.1.2 Metodo de Euler Implícito (ou Regresivo)

Procedendo ao desenvolvimento em série de Taylor da solução exacta de (1.8) em torno

de nx , em que ao contrário do método progressivo em que usamos uma diferença

Progressiva, agora vamos utilizar uma diferença regressiva, i.e.,

( ) ( ) ( ) ( ) ( )2

3

2!n n n n

hy x h y x hy x y x O h′ ′′− = − + − (2.10)

Desprezando os termos de ordem 2≥ em (2.10), podemos obter para cada n

x , a

equação de diferenças:

1 1, 0,1,2,..., 1n n n

y y hf n N+ += + = − (2.11)

que corresponde ao método de Euler implícito (ou Regressivo).

Notemos que o termo que queremos determinar ( 1ny + ) está presente em ambos os

membros de (2.11) (pois ( )1 1 1,

n n nf f x y+ + += ), pelo que ele é determinado de forma

implícita , daí a designação de método implícito.


Como a função f é geralmente não linear, para determinar 1n

y + devemos resolver a

equação algébrica

( ) ( ) ( )1 1 1 1 1, com ,n n n n n n

y y y y hf x y+ + + + += ϕ ϕ = + (2.12)

A forma da equação acima sugere a utilização do método de ponto fixo, i.e., iterar de

acordo com

( ) ( )( )1

1 1 , com 0,1, 2,...k k

n ny y k

+

+ += ϕ = (2.13)

a partir de uma estimativa inicial ( )0

1ny + (que pode ser obtido utilizando o método de Euler

explícito).(PINA, 1995)

Erros de Aproximação

Raramente é possível, com realismo, seguir a propagação do erro num algorítimo de

aproximação de um PVI, mas é possível obter algumas estimativas desse erro, se bem

que grosseiras.

O Erro de truncatura local (ou erro de discrtização local) ( )1nh+Τ é o erro que se

comete ao passar do nón

x para o nó 1n

x + , quando substituimos um processo infinito por

um processo finito (por exemplo no desenvolvimento em série de taylor , quando

retivermos os termos até à ordem 2 e negligenciarmos os restantes termos) (KINCARD,

2002). Define-se o erro de discretização local por unidade de passo ( 1( )

nt h+ ) como o

quociente entre o erro de truncatura local e o tamanho do passo h ( i.e., 11

( )n

n

ht

h

++

Τ= ).

A diferença entre a solução exacta e a solução aproximada ( ( ) , 0,1,...,n ny x y n N− = ) é

denominado por Erro de truncatura global (ou erro discretização global) ( )( )ne h .Os

erros de arredondamento ocorrem porque os computadores trabalham com um

número fixo de dígitos e nas multiplicações e divisões perdem-se os dígitos que

ultrapassam esse número.


Devido a dificuldade da análise dos erros de arredondamento, nesse trabalho

daremos atenção somente para os erros de truncatura (com especial atenção para o erro

de truncaturra local). Indicamos [PIN95], [QUO07], ou [KRE98] se o leitor estiver

interessado numa análise mais profunda do assunto.

2.1.1.3 Erros dos métodos de Euler

Em ambos os métodos de Euler, no desenvolvimento em série de Taylor da solução

exacta ( )ny x nós negligenciamos os termos de ordem 2≥ , então, o erro de truncatura

local será a soma de todos os termos (ou resíduos) que não incluímos na solução

aproximada n

y (esses resíduos são representados por ( )1pO h + †.

O termo

( ) ( ) ( )2

2

12

n n

hh y O h+

′′Τ = =ξ (2.14)

Corresponde a termos truncados em (2.9) no desenvolvimento em série de Taylor e

portanto é o erro de truncatura local para o método de Euler progressivo. Para o método

de Euler regressivo o erro de truncatura local será

( ) ( ) ( )2

2

1 1, 2

n n n n n

hh y O h x xξ ξ+ +

′′Τ = − = < <

Vimos anteriormente que o erro de truncatura local para os métodos de Euler é

dado por ( )2O h . Se esse foi o erro cometido em cada passo, então no final do intervalo

0x X, , depois de N passos, o erro de truncatura global ( ( ) , 0,1,...,n ny x y n N− = ) para

o método de Euler progressivo será aproximadamente: (MATHEWS & FINK, 2004)

( ) ( ) ( )( ) ( ) ( )

2 20 1

1 2 2 2 2

N

k

k

X x yh h hNy Ny y h h O h

=

′′−′′ ′′ ′′≈ = = =∑

ξξ ξ ξ

† Relembremos que ( )p

T O h= se existirem constantes positivos p e h tais que para

0h > suficientemente pequeno , 0pT Ch C≤ > .


Analogamente encontramos que para o método de Euler regressivo ( ) ( )1

ne h O h=

Em relação ao erro de truncatura global ( ( ) , 0,1,...,n ny x y n N− = ), raramente é

necessário saber o seu valor num ponto arbitrário do intervalo. Normalmente interessa-

nos conhecer o erro de truncatura global somente no final do intervalo de integração

(LINZ & WANG, 2003).

Dum modo geral, se o erro de truncatura local para um determinado método for

( )1pO h + , o seu erro de truncatura global ( ) ( )n n ne h y x y= − será dado por ( )pO h , pois,

enquanto que ( )+Τ1n

h é o erro cometido quando passamos do ponto ( )n nx y, para o

ponto ( )+ +1 1n nx y, , o erro de truncatura global ( )n

e h é a acumulação dos n erros (de

truncatura local) cometidos nas sucessivas passagens desde o ponto ( )0 0x y, até o ponto

( )n nx y, (tendo em conta que, ( )=

0 0y y x ). (KINCARD, 2002)

Os métodos de Euler têm como vantagem a fácil aplicação, mas como são métodos de

primeira ordem, temos de recorrer a valores de tamanho de passo muito pequenos para

podermos obter uma aproximação razoável (como é visível no exemplo 2.5). Por isso,

para construir a solução aproximada num dado intervalo [ ]0 ,I x X= , são necessários

muitos passos, o que torna o processo lento, por isso, são métodos pouco utilizados na

prática.

2.1.2 Método de Heun

O método de Heun‡ pode ser obtido pela combinação dos métodos de Euler explícito

e implícito

Somando (2.9) e (2.11) vem:

( )1 1 1

1 12 2

n n n n n n

n n n n

y y y y hf hf

y y h f f

+ + +

+ +

+ = + + + ⇔

⇔ = + +

‡ Também conhecido como método de Euler Modificado ou ainda método trapezoidal, foi desenvolvido palo matemático alemão Karl Heun (1859-1929).


Dividindo por 2 obtemos

( )1 12

n n n n

hy y f f+ += + + (2.15)

Tal como acontece com método de Euler implícito, também no método (2.15) o termo

1ny + está presente em ambos os membros da igualdade, daí tratar-se também de um

método implícito. Substituindo em (2.15) ( )1 1 1,n n n

f f x y+ + +≡ por ( )1,n n n

f x y hf+ + , pela

utilização do método de Euler explícito (i.e., 1n n n

y y hf+ = + ) obtemos

( )1 1( ,2

n n n n n n

hy y f f x y hf+ += + + + (2.16)

que é denominado método de Heun.

2.1.2.1 Erros do método de Heun O método de Heun também pode ser deduzido a partir de (1.9), em que o integral

presente no lado direito da igualdade pode ser aproximado pela regra do trapézio entre

os pontos 1 e n nx x + originando um erro de truncatura local (MATHEWS & FINK, 2004)

( ) ( ) ( )3

3

1 1, 12

n n n n n

hh y O h x xξ ξ+ +

′′′Τ = − = < < (2.17)

Como se vê acima, o erro de truncatura local para o método de Heun é de maior ordem

que os métodos de Euler, daí esse método produzir melhores aproximações que os

métodos de Euler.

2.1.3 Métodos de Taylor de ordem mais elevada

Um meio de obter melhores aproximações é utilizando mais termos no

desenvolvimento em série de Taylor da solução exacta ( )y x de (1.8), originando assim

um erro de truncatura de ordem mais elevada.


Por exemplo, para conseguir um método de Taylor de 2ª ordem consideremos o

seguinte desenvolvimento

( ) ( ) ( ) ( ) ( )2

3

2n n n n

hy x h y x hy x y x O h′ ′′+ = + + + (2.18)

Contudo não temos directamente ( )y x′′ . A partir de ( ),y f x y′ = temos de derivar f em

ordem a x .

Define-se o método de Taylor de segunda ordem pela fórmula

( ) ( )2

1 , , , 0,1,.., 12

n n n x n n y n n n

hy y hf f x y f x y f n N+

′ ′ = + + + = − (2.19)

onde

( ),n n n

f f x y= e ( ) ( ) ( ) ( ), , ,x n n y n n n

df x y f x y f f x y y x

dx′ ′ ′′ + = = .

De modo Semelhante definem-se os métodos de Taylor de ordem superior a segunda. Os

métodos de Taylor de ordem mais elevada (superior a primeira), são muito pouco

utilizados, pois à medida que aumentamos a ordem, as derivadas tornam-se mais

complexas.

2.1.3.1 Erros

No método de Taylor, quando negligenciamos os termos de ordem maior que p , então o

erro de truncatura local será a soma de todos os termos que nós não incluímos na

solução aproximada.

Pela série de Taylor esses termos são

( )( )

( ) ( ) ( )1

1 1

11 !

pn p

n

hh y O h

pξ

++ +

+Τ = =+

(2.20)

2.1.4 Análise da Consistência dos métodos da série de Taylor

Definição 2.3 (Consistência) Um método diz-se consistente se (PINA, 1995)


( )10

lim 0nh

t h+→

= (2.21)

onde ( ) ( )0

1 1

1

N

n n

x n x

t h Sup hh

+ +

≤ ≤

= Τ

e diz-se que a sua ordem de consistência é 0p > se

1 , 0p

nt ch c+ ≤ ≤ ≤ ∞ (2.22)

Para h suficientemente pequeno e c não dependente de h , i.e., ( ) ( )1

p

nt h O h+ = .

Em geral, um método diz-se consistente se o erro de discretização local por unidade de

passo for um infinitésimo com respeito a h .

Ambos os métodos de Euler são métodos consistentes com ordem de consistência igual a

1, enquanto o método de Heun tem ordem de consistência igual a 2. Os métodos de

Taylor de ordem mais elevada têm ordem de consistência igual a 2,3,4,... se são

respectivamente de 3ª ordem, 4ª ordem, 5ª ordem, …. (PINA, 1995).

Antes de introduzirmos os métodos de Runge Kutta vamos descrever brevemente a

generalização dos métodos de passo simples. No geral esses métodos podem ser escritos

na forma (QUORTERONI, SACCO, & SALERI, 2007)

( )1 , , ;n n n n ny y h x y f h+ = + φ (2.23)

Onde φ é normalmente designado por função incremental. Podemos escrever

( ) ( ) ( ) ( )1

1, , ; 0 1p

n n n n ny x y x h x y f h O h n N

+

+ = + + ≤ ≤ −φ (2.24)

Onde ( ) ( )1

1

p

nO h h+

+= Τ é o resíduo no nó 1nx + quando no esquema numérico é utilizada

a solução exacta.

Por exemplo nos métodos acima apresentados temos que

• Para o método de Euler progressivo ( )( , , ; ) , ;n n n n n

x y f h f x yφ =


• No método de Heun ( ) ( )( )1

, , ,( , , ; )

2

n n n n n n

n n n

f x y f x y hf x yx y f hφ

++ +=

• Para o método de Taylor de 2ª ordem

( ) ( ) ( )( , , ; ) , , ,2

n n n n n x n n y n n n

hx y f h f x y f x y f x y fφ ′ ′ = + + .

É fácil ver que para os métodos acima apresentados, o erro de truncatura local é dado

por

( ) ( ) ( )1 1 ( , , ; )n n n n n nh y x y x h x y f h+ +Τ = − − φ (2.25)

2.2 Métodos de Runge Kutta (RK)

A quando da apresentação dos métodos de Taylor de ordem mais elevada ( )2p ≥ , vimos

que o método de Taylor de 2ª ordem ( )2p = tomava a seguinte forma

( ) ( )2

1 , ,2

n n n x n n y n n n

hy y hf f x y f x y f+

′ ′ = + + + (2.26)

Esse método tem a grande inconviniência de necessitar do cálculo das derivadas. Para

evitar as derivações de f recorre-se normalmente aos chamados métodos de Runge

Kutta §.

Os métodos de Runge Kutta de 2ª ordem têm a seguinte forma

( )

( )

( )

1 1 1 2 2

1

2 2 21 1

,

,

n n

n n n

n n

y y h F F

F f f x y

F f x h y h F

ω ω

α β

+ = + +

= =

= + +

(2.27)

§ Esses métodos foram desenvolvidos por volta de 1900 pelos matemáticos Germânicos Carl. Runge (1856-1927) e Martin Wilhelm. Kutta (1867-1944)


As constantes 1 2 2 21, , ,ω ω α β devem ser determinadas de modo em que a diferença entre

os 2os membros de (2.26) e da primeira equação de (2.27) seja O ( )3h quando 0h → .

Para isso, vamos desenvolver 2

F em série de Taylor para função de duas variáveis

( ) ( )2 2 21 1, ,n x n n y n nF f hf x y hF f x yα β′ ′= + + (2.28)

Substituindo (2.28) na 1ª equação de (2.27) obtemos

( ) ( ) ( )2

1 1 2 2 2 2 21, ,

n n n x n n y n n ny y h f h f x y f x y fω ω ω α ω β+

′ ′ = + + + + (2.29)

Comparando (2.29) com (2.26) concluimos que

1 2 2 2 2 21

11

2ω ω ω α ω β+ = = =e (2.30)

As equações (2.30) não determinam os parâmetros 1 2 2 21, , ,ω ω α β de forma única, pelo

que existem várias possibilidades de métodos de Runge Kutta de 2ª ordem. (PINA, 1995)

Vejamos alguns exemplos

Tomando 1 2 2 21

1 e 1

2ω ω α β= = = = obtemos o método de Heun.

( )

( )

( )

1 1 2

1

2 1

2

,

,

+ = + +

= =

= + +

n n

n n n

n n

hy y F F

F f f x y

F f x h y hF

(2.31)

Outra possibilidade é tomar 1 2 2 21

10, 1 e

2ω ω α β= = = = , obtendo deste modo o

método conhecido como o método de Euler modificado (ou método de Euler Cauchy).

( )

1 2

12

1

,2 2

,

+ = +

= + +

=

n n

n n

n n

y y hF

hFhF f x y

F f x y

(2.32)


No geral a família dos métodos de Runge Kutta é dada pelas seguintes fórmulas (PINA,

1995)

( )

1

1

1

,

ω

α

β

+=

=

= +

= +

= +

=

∑

∑

s

n n i i

i

i n i

s

i n im m

m

i i i

y y h F

X x h

Y y h F

F f X Y

(2.33)

Nas fórmulas acima s indica o número de estágios do método (ou seja número de

funções a calcular), os valores , e i i im

ω α β são entendidos como parâmetros próprios de

cada método. Normalmente estes parâmetros são organizados no quadro da figura 2.1,

que é conhecido como quadro de Butcher.

1α

11β

12β …

1sβ

2α

21β

22

β … 2s

β

…

…

… …

sα

1sβ

2sβ …

ssβ

1

ω 2ω …

sω

Figura 2.1 Quadro de Butcher

Um método de Runge Kutta diz-se explícito se a matriz dos “ β ” da figura 2.1 for

estritamente triangular inferior (figura 2.2), i.e. ( )0 para j , , 1,...,ij

i i j sβ = ≥ = e implícito

se, esta matriz não for triangular inferior, i.e., existe pelo menos um j i> tal que 0ijβ =/ .

Nesse trabalho somente estudaremos os métodos de Runge Kutta explícitos, indicando

[BUT08] caso ao leitor estiver interessado no estudo dos métodos implícitos.

0 0 0 … 0

2α

21β

0 … 0

…

…

… …

sα 1s

β 2sβ … 0

1

ω 2ω …

sω

Figura 2.2


Como aconteceu com os métodos da série de Taylor, nós dizemos que (2.33) é um

método com ordem de consistência ( )1p p ≥ se ( ) ( )1

1

p

n h O h+

+Τ = , com 0h → .

2.2.1 Métodos de Runge Kutta de primeira ordem ( )1p =

Fazendo 1s = nas fórmulas (2.33) obteremos o método de Runge Kutta de primeira

ordem que é dado por

( )1 1 1 1, ,ω+ = + =n n n ny y h F F f x y (2.34)

Nesse caso temos apenas um parâmetro 1ω para determinar. Procedendo de forma

semelhante ao método de Runge Kutta de 2ª ordem visto anteriormente obteremos

11ω = . O método resultante será

( )

1 1

1 ,

+ = +

=

n n

n n

y y hF

F f x y (2.35)

que corresponde ao método de Euler progressivo.

Se aumentarmos o número de estágios ( )s obteremos novos métodos.

2.2.2 Métodos de Runge Kutta de terceira ordem 3p =

Se procedermos de forma análoga ao método de 2ª ordem ( )2p = podemos obter os

métodos de Runge Kutta de 3ª ordem ( )3p = .

A fórmula geral de um de um método de Runge Kutta de 3ª ordem é dada por

( )

( )

( )

( )

1 1 1 2 2 3 3

1

2 2 21 1

3 3 31 1 32 2

,

,

,

ω ω ω

α β

α β β

+ = + + +

=

= + +

= + + +

n n

n n

n n

n n

y y h F F F

F f x y

F f x h y h F

F f x h y h F h F

(2.36)

Como exemplo de métodos de Runge Kutta de 3ª ordem temos


Método de Kutta: (PINA, 1995)

( )

( )

( )

1 1 2 3

1

12

3 2 1

46

,

,2 2

, 2

n n

n n

n n

n n

hy y F F F

F f x y

hFhF f x y

F f x h y hF hF

+ = + + +

=

= + +

= + + −

(2.37)

Método de Heun de terceira ordem (PINA, 1995)

( )

( )

1 1 3

1

12

23

34

,

,3 3

22,

3 3

n n

n n

n n

n n

hy y F F

F f x y

hFhF f x y

hFhF f x y

+ = + +

=

= + +

= + +

(2.38)

2.2.3 Método de Runge Kutta de 4ª ordem (RK4)

A formula geral dos métodos de Runga Kutta de 4ª ordem é dada por

( )

( )

( )

( )

( )

1 1 1 2 2 3 3 4 4

1

2 2 21 1

3 3 31 1 32 2

4 4 41 1 42 2 43 3

,

,

,

,

n n

n n

n n

n n

n n

y y h F F F F

F f x y

F f x h y h F

F f x h y h F h F

F f x h y h F h F h F

+ = + + + +

=

= + +

= + + +

= + + + +

ω ω ω ω

α β

α β β

α β β β

(2.39)

Como exemplo vamos apresentar a fórmula de Runge Kutta mais utilizada (PINA, 1995)

( )

( )

( )

1 1 2 3 4

1

12

23

4 3

2 26

,

,2 2

,2 2

,

n n

n n

n n

n n

n n

hy y F F F F

F f x y

hFhF f x y

hFhF f x y

F f x h y hF

+ = + + + +

=

= + +

= + +

= + +

(2.40)


Por ser de entre todos os métodos de Runge Kutta a mais utilizada, ela é designada

simplesmente por método de Runge Kutta. Essa popularidade é devido ao bom

equlíbrio entre o esforço computacional (em cada passo efectuamos quatro cálculos de

f ) e a elevada precisão iminente a este método. (PINA, 1995)

A partir da quinta ordem, são necessárias mais cálculos de f (ou seja mais estágios

s ) do que a ordem p , o que se traduz numa perda de eficácia. Por exemplo, para o

método de Runge Kutta de 5ª ordem precisamos efectuar seis cálculos de f , i.e., a

ordem 5p = só se obtém com seis estágios ( )6s = (PINA, 1995). A tabela seguinte

mostra a relação entre a ordem e o número mínimo de estágios para os métodos de

Runge Kutta (QUORTERONI, SACCO, & SALERI, 2007)

Ordem

p

1 2 3 4 5 6 7 8

mins 1 2 3 4 6 7 9 11

Tabela 2.2 Relação entre a ordem e o número de estágios

A partir da 5ª ordem os métodos de Runge Kutta perdem eficácia, e por isso tornam-se

menos atractivos e económicos, daí a razão da grande popularidade dos métodos com

ordem 4, pois requerem apenas 4 estágios, ou seja, quatro cálculos de f por passo.

Daquilo que acabamos de ver, podemos concluir que quanto maior for a ordem dos

métodos de Runge Kutta maior será o número de vezes em que é calculado ( ),f x y (ou

seja maior o número de estágios) num passo. Para remediar esse inconveniente

apresentaremos no próximo capítulo uma outra classe de métodos.

2.2.4 Erros dos métodos de Runge Kutta

Como os métodos de Runge Kutta correspondem aos métodos da série de Taylor, um

método de Runge Kutta de ordem p tem erro de truncatura local


( ) ( )1

1

p

nh O h

++Τ = (2.41)

Em termos de consistência, todos os métodos de Runge Kutta apresentados

anteriormente são consistentes. (PINA, 1995)

Para melhores explicações sobre o erro dos métodos de Runge Kutta indicamos [QUO07]

ou [BUT08], que apresentam esse assunto de forma mais detalhada.

Exemplo 2.6 Utilizando o método de Runge Kutta de 4ª ordem resolva o PVI

( ) ( ) [ ]2, , com 0 1, em 0,3f x y x y y I= − = = .Compare os resultados para

1 1 11, , ,

2 4 8h = .

A solução exacta do PVI é ( ) 2 2 2xy x e x x−= − + − + . A tabela (2.3) mostra a solução exacta

do PVI e as soluções aproximadas para os diferentes valores do tamanho de passo

(utilizando o PROGRAMA 4 do anexo). Por exemplo para 1

8h = temos

( )( )

( ) ( ) ( ) ( )( )

2

1

2

2

2

3

2

4

1

0 1 1

0.0625 (1 0.125*0.5*( 0.5)) 0.9648

0.0625 (1 0.125*0.5*( 0.9648)) 0.9358

0.125 1 0.125* 0.9358 0.8674

0.1251 1 2 0.9648 2 0.9358 0.8674 0.8831

6

= − = −

= − + − = −

= − + − = −

= − + − = −

= + − + − + − + − =

F

F

F

F

y

nx n

y ( )ny x Exacta

h = 1 .h = 0 5 .h = 0 25 .= 0 125h

.

.

.

.

.

.

.

.

.

.

0

0 125

0 25

0 375

0 50

0 75

1 00

1 50

2 00

2 50

3 00

.

.

.

.

1 0

0 6458

1 8880

4 9788

.

.

.

.

.

.

.

1 0

0 6439

0 6329

1 0279

1 8659

3 1693

4 9517

.

.

.

.

.

.

.

.

.

1 0

0 7837

0 6435

0 5902

0 6322

0 0269

1 8647

3 1680

4 9503

.

.

.

.

.

.

.

.

.

.

.

1 0

0 8831

0 7837

0 7033

0 6435

0 5901

0 6321

1 0269

1 8647

3 1679

4 9502

.

.

.

.

.

.

.

.

.

.

.

1 0

0 8831

0 7836

0 7033

0 6434

0 5901

0 6321

1 0268

1 8446

3 1679

4 9502



Comparando os resultados do PVI dada no exemplo acima, é possível ver que o

método de Runge Kutta de 4ª ordem oferece uma aproximação da solução exacta muito

melhor que o método de Euler Explícito (ver exemplo 3.5). Essa precisão, como já foi

dito anteriormente é devido ao facto do método de Runge Kutta de 4ª ordem utilizar

mais termos no desenvolvimento em série de Taylor da solução exacta do que o método

de Euler explícito.

No exemplo anterior precisamos de efectuar vários cálculos de f (quatro ao todo)

para obtermos a primeira ( 1y ) de 24 aproximações. Essa tarefa seria ainda bem mais

cansativa se a nossa função fosse mais complexa. O elevado número de cálculos da

função f por passo, é uma das grandes desvantagens dos métodos de Runge Kutta, mas

com o advento da informática, a penosa tarefa dos cálculos foi substituída pelos

computadores.

2.3 Análise da Convergência dos métodos de passo simples

Terminado a apresentação dos métodos de passo simples, vamos agora analisar a

convergência desses métodos, i.e., investigar em que condições o erro de truncaturra

global tende para zero com o parâmetro h .

Para o estudo da convergência dos métodos de passo simples, vamos iniciar com a

seguinte definição.

Definição 2.4 Diz-se que um método de passo simples satisfaz a condição de Lipchitz, se

a função incremental φ verificar a relação

( ) ( ) [ ]0, , , ,

hx v x L v x I x X− ≤ − ∀ ∈ =φ φ ω ω (2.42)

para h suficientemente pequeno. (PINA, 1995)

O seguinte teorema estabelece as condições de convergência dos métodos de passo

simples


Teorema 2.1 Se um método de passo simples satisfazer a condição de Lipchitz, a

consistência é uma condição necessária e suficiente para convergência da solução

aproximada n

y para a solução exacta ( )ny x , i.e., (PINA, 1995)

( )( ) ( )10 0

lim 0 lim 0+→ →

= ⇔ − =n n nh h

t h y y x (2.43)

Demonstração. Ver (PINA, 1995)

A partir do teorema acima, fica claro que para demonstrar a convergência dos métodos

de Taylor e dos métodos de Runge Kutta estudados anteriormente, é preciso provar

apenas que as respectivas funções incrementais satisfazem a condição de Lipchitz, pois,

como foi dito anteriormente esses métodos são todos consistentes.

Todos os métodos de passo simples apresentados são métodos consistentes e

convergentes. Indicamos [PIN95], [QUO07] ou [BUT08] para uma análise mais

aprofundada do estudo da consistência e convergência dos métodos de passo simples.


CAPÍTULO III

MÉTODOS DE PASSO MÚLTIPLO

Como vimos anteriormente, nos métodos de passo simples o valor 1ny + é

calculado a partir do valor de ny utilizando apenas informações no intervalo [ ]1,n nx x + .

Quando efectuamos o passo 1ny + , já dispomos dos valores aproximados 0 1

, ,...,n

y y y . A

ideia dos métodos de passo múltiplo (ou métodos multipasso) é recorrer às

informações obtidas em vários passos anteriores para obter uma melhor aproximação

de 1ny + (PINA, 1995). Se recorrermos à informação obtida nos dois passos anteriores

temos os métodos de 2 passos, se a informação for obtida nos 3 passos anteriores temos

os métodos de 3 passos, …, se a informação for obtida nos q passos anteriores temos os

métodos de q passos ( 1q ≥ ). Geralmente temos

Definição 3.1 Um método de q passos ( 1q ≥ ) é um método que para 11,

nn q y +∀ ≥ −

depende dos valores 1,...,n n qy y + − , mas não depende dos valores ky com 1k n q< + − .

(PINA, 1995)


Exemplo 3.1 Um método explícito de 2 passos muito conhecido pode ser obtido

utilizando a seguinte aproximação para a primeira derivada de (1.8) (QUORTERONI,

SACCO, & SALERI, 2007)

( ) 1 1

2

n ny y

y xh

+ −−′ ≈

originando

1 1 ,2 n 1,n n ny y hf+ −= + ≥

Onde ( )0 0 1, y y x y= devem ser determinados por um qualquer método de passo simples

(normalmente utiliza-se o método de Runge Kutta de quarta ordem) e ( ),n n nf f x y≡ .

Esse método é conhecido como método do ponto médio.

Exemplo 3.2 Um método implícito de dois passos é o método de Simpson

(QUORTERONI, SACCO, & SALERI, 2007), obtido de (1.9), com 0 1 1

e n n

x x x x− += = e

utilizando o método de Cavalieri- Simpson para aproximar o integral de f , i.e.,

( )( ) ( )1

11 1

, 4 ,3

n

n

x

n n nx

hf s y s ds f f f

+

−− +≈ + +∫

originando o método

[ ]1 1 1 14 , 1,

3n n n n n

hy y f f f n+ − − += + + + ≥

onde tal como no exemplo anterior ( )0 0y y x= e 1y têm de ser determinados por um

qualquer dos métodos de passo simples.

Observação 3.1 Dos exemplos que acabamos de ver fica claro que um método de

q passos− precisa de q valores iniciais 0 1 1, ,..., qy y y − para poder ser inicializado (esses

valores são determinados pelos métodos de passo simples, sendo o método de Runge Kutta

de quarta ordem a mais utilizada).


Existem, vários tipos de métodos multipassos, mas nós vamos abordar somente a

família de métodos multipassos conhecidas como métodos de Adams, indicando [QUO07]

ou [PIN95] para o estudo de outros métodos.

3.1 Métodos de Adams Bashforth (AB)

Relembremos que a solução exacta do PVI (1.8) é dada por

( ) ( ) ( )( )1

1 ,+

+ = + ∫n

n

x

n nx

y x y x f s y s ds

Nos métodos de Adams Bashforth** a aproximação do integral da direita da equação

acima é feita recorrendo a um polinómio interpolador da integranda ( )( ),f x y x nos nós

n mx − , 1

,...n m

x − + , 1,

n nx x− i.e., um polinómio mp de grau m≤ . Seja

( ) ( )0

m

m n j n j

j

p x L x f− −=

=∑ (3.1)

a representação de Lagrange desse polinómio interpolador, então (RAINER, 1998)

( )( ) ( )1 1

,n n

n n

x x

mx x

f x y x dx p x dx+ +

≈∫ ∫ (3.2)

em que

( ) ( )1 1

0 0

n n

n n

m mx x

m n j n j j n jx x

j j

p x dx L x f dx h f+ +

− − −= =

= =∑ ∑∫ ∫ γ (3.3)

Os coeficientes jγ são definidos por (PINA, 1995)

( )11 +

−γ = ∫n

n

x

j n jx

L x dxh

(3.4)

A expressão geral dos métodos de Adams Bashforth é

** Método desenvolvido por John Cauch Adams (1819-1892) para resolver uma equação diferencial criada por Francis Bashforth (1819-1912)


1

0

+ −=

= + γ

∑

m

n n j n j

j

y y h f (3.5)

que são métodos de 1m + passos.

Da fómula acima concluimos que estes métodos são explícitos.

Exemplo 3.3 Deduzir o método de Adams Bashforth para 1m = e o seu erro de truncatura

local.

Para 1m = temos que ( ) 11 1

1 1

n nn n

n n n x

x x x xp x f f

x x x x−

−− −

− −= +

− − é o polinómio que interpola f

nos pontos ( ) ( )1 1, , e n n n nx f x f− − .

Efectuando os cálculos obtemos

10

1

1

1

1 3

2

1 1

2

n

n

n

n

x hn

xn n

x hn

xn n

x xdx

h x x

x xdx

h x x

+−

−

+

−

−γ = =

−

−γ = = −

−

∫

∫

pelo que temos

1 1

3 1, 1

2 2n n n ny y h f f n+ −

= + − ≥

Para obtermos o erro de truncatura local, notemos que o erro de interpolação

polinomial é dado por (PINA, 1995)

( )( ) ( )( )

( ) ( )( ) ( )11

1, , ,

1 !m

m m n nf x y x p x f y W x x xm

+−− = ξ ξ < ξ <

+ ,

onde ( )( ) ( )1 .....m n n n mW x x x x x x− −= − − − .

De acordo com que foi dito acima, para 1m = temos

( )( ) ( ) ( )( )( )( )1 1 1

1, , ,

2!n n n nf x y x p x f y x x x x x x− −′′− = ξ ξ − − < ξ < ,

Pelo que

( )( ) ( ) ( )( ) ( )( )1 1

1, , ,

2!

n n

n n

x h x h

n nx x

f x y x p x f y x x x x dx+ +

−′′− = ξ ξ − −∫ ∫


Depois do cálculo do integral acima obtemos

( ) ( )( )31

5,

12n h h f y+ ′′Τ = ξ ξ .

Com ordem de consistência igual a 2.

Apresentamos a seguir alguns métodos de Adams Bashforth para diferentes valores de

m e os respectivos erros de truncatura local (PINA, 1995).

Para 0m =

( ) ( )

1

2

1

0

1

2

n n n

n n

y y hf n

h h f

+

+

= + ≥

′Τ = ξ (3.6)

Para 1 m =

( )

( ) ( )

1 1

3

1

3 12

5

12

n n n n

n n

hy y f f n

h h f

+ −

+

= + − ≥

′′Τ = ξ

(3.7)

Para 2m =

( )

( ) ( )

1 1 2

4

1

23 16 5 212

3

8

n n n n n

n n

hy y f f f n

h h f

+ − −

+

= + − + ≥

′′′Τ = ξ

(3.8)

Para 3m =

( )

( ) ( ) ( )

1 1 2 3

5

1

55 59 37 9 324

251

720

n n n n n n

iv

n n

hy y f f f f n

h h f

+ − − −

+

= + − + − ≥

Τ = ξ

(3.9)

Para calcular o valor 1ny + , os métodos de AB precisam dos valores 1

, ,...,n m n m n

y y y− − + .

Isto quer dizer que o primeiro valor a ser calculado é my . Os valores 0 1 1

, ,...,m

y y y − devem

ser calculados por um método auxiliar que deverá ser de passo simples (normalmente é

utilizado o método de Runge Kutta de 4ª ordem com igual ordem de consistência).


Uma das desvantagens dos métodos de passo múltiplo é a de que não são “auto-

iniciáveis”, i.e., precisam de alguns valores pré computados para poderem iniciar, mas ao

contrário dos métodos de Runge Kutta esses métodos têm a vantagem de exigirem

somente um cálculo de f por passo, visto que os restantes valores são aproveitados dos

passos anteriores (PINA, 1995).

3.2 Métodos de Adams Moulton (AM)

Para os métodos de AB, o polinómio ( )mp x interpolador nos nós ,...,

n m nx x− é usada

para aproximar a integrada do integral em (1.9) entre os nós n m

x − , 1,...

n mx − + ,

1,

n nx x− . Para

uma melhor aproximação, os métodos de Adams Moulton†† surgiram com a ideia de

estender o polinómio interpolador até o nó 1n

x + , i.e., utilizar os nós

1 1 1, ,..., , ,

n m n m n n nx x x x x− − + − + para construir um polinómio interpolador ( )1mp x+ de grau

1m≤ + .

Utilizando o processo adoptado nos métodos de AB chegamos a expressão geral dos

métodos de Adams Moulton. (PINA, 1995)

1

1 1

0

+

+ − +=

= + γ

∑m

n n j n j

j

y y h f (3.10)

que são métodos de 2m + passos.

Os métodos de AM são implícitos, pois, no segundo membro de (3.10) figura o valor

( )1 1 1,n n nf f y+ + += .

Exemplo 3.4 Deduzir o método de Adams Moulton para 1m = e o respectivo erro de

truncatura local.

†† Tal como os métodos de Adams Bashforth, esses métodos também foram desenvolvidos por John Couch Adams. O nome de Forest Ray Moulton (1872-1952) ficou associado a esses métodos, pois ele observou que eles podiam ser usados em conjunto com os métodos de Adams Bashforth, originando os pares preditores correctores.


Para 1m = temos que ( )( )( )

( )( )( )( )

( )( )1 1 1

2 1

1 1 1 1 1

n n n n

n n

n n n n n n n n

x x x x x x x xp x f f

x x x x x x x x+ − +

−− − + − +

− − − −= + +

− − − −

( )( )( )( )

11

1 1

n n

n

n n n n

x x x xf

x x x x−

++ −

− −+

− − é o polinómio que interpola f nos pontos

( ) ( ) ( )1 1 1 1, , , , e n n n n n nx f x f x f− − + + .

Efectuando os cálculos obtemos

( )( )( )( )

( )( )( )( )

( )( )( )( )

10

1 1 1

11

1 1 1

1 12

1 1

1 5

12

1 8

12

1 1

12

n

n

n

n

n

n

x hn n

xn n n n

x hn n

xn n n n

x hn n

xn n n n

x x x xdx

h x x x x

x x x xdx

h x x x x

x x x xdx

h x x x x

+−

+ + −

++

− − +

+− +

− +

− −γ = =

− −

− −γ = =

− −

− −γ = = −

− −

∫

∫

∫

pelo que temos

( )1 1 15 812

n n n n n

hy y f f f+ + −= + + −

Para obtermos o erro de truncatura local, notemos que o erro de interpolação

polinomial é dado por (PINA, 1995)

( )( ) ( )( )

( ) ( )( ) ( )21 1 1 1

1, , ,

2 !m

m m n nf x y x p x f y W x x xm

++ + − +− = ξ ξ < ξ <

+ ,

então, para 1m = temos

( )( ) ( ) ( )( )( )( )( )2 1 1

1, , ,

3!n n nf x y x p x f y x x x x x x− +′′′− = ξ ξ − − − ,

Pelo que

( )( ) ( ) ( )( ) ( )( )( )1 1 1

1, , ,

6

n n

n n

x h x h

n n nx x

f x y x p x f y x x x x x x dx+ +

− +′′′− = ξ ξ − − −∫ ∫

Depois do cálculo do integral acima obtemos

( ) ( )( )41

1,

24n h h f y+ ′′′Τ = − ξ ξ ,

i.e. tem ordem de consistência igual a 3.


Apresentamos a seguir alguns métodos de Adams Moulton para diferentes valores de m

e seus respectivos erros de truncatura local. (PINA, 1995)

Para 0m =

( )

( ) ( )

1 1

3

02

1

12

+ += + + ≥

′′Τ = − ξ

n n n n

n n

hy y f f n

h h f

(3.11)

Para 1m =

( )

( ) ( )

1 1 1

4

5 8 112

1

24

+ + −= + + − ≥

′′′Τ = − ξ

n n n n n

n n

hy y f f f n

h h f

(3.12)

Para 2m =

( )

( ) ( ) ( )

1 1 1 2

5

9 19 5 224

19

720

+ + − −= + + − + ≥

Τ = − ξ

n n n n n n

iv

n n

hy y f f f f n

h h f

(3.13)

Como os métodos de AM são implícitos, temos de calcular ( )1 1 1,n n nf f y+ + +≡

recorrendo a uma equação algébrica não linear

( )1 1n ny y+ += ϕ (3.14)

Com ( ) ( )1

1 1 0 1 1

1

,m

n n j n j n n

j

y y h f h f x y+

+ − + + +=

ϕ = + γ + γ

∑ (PINA, 1995)

Normalmente utiliza-se o método iterativo linear do ponto fixo em que a aproximação

inicial de 1ny + (designemos por ( )

1

0ny + ) é dada por um dos métodos de AB.

A combinação dos métodos acima é designada por par preditor corrector, pois, o

método de AB funciona como fórmula preditora de 1ny + , enquanto as iterações de ponto

fixo do método implícito de AM actuam como correctoras de 1ny + . (PINA, 1995)


3.3 Processo preditor corrector

Como vimos anteriormente, nos métodos multipassos implícitos é necessário

utilizar um método auxiliar para calcular a estimativa inicial ( )0

n my + . Esse método utilizado

(geralmente explícito) para estimar ( )0

n my + é conhecido como preditor, sendo o método

linear propriamente dito designado por corrector.

O método explícito produz um valor aproximado de n my + (designemos por ( )0

n my + ), logo a

seguir o valor ( )0

n my + é utilizado no método implícito, gerando um valor melhorado (ou

corrigido) de n my + .

Como exemplo de método preditor corrector temos

Método de Adams-Bashforth-Moulton (ABM) (MILLER, 1991)

0

1 2 3

é dado;

, , são determinados por um método de passo simples

y

y y y

( ) ( )

( )( )( )

1

1

*

1 2 3

*

1 1 1 2

para 3,4,..., 1, temos

(Preditor) 55 59 37 924

(Corrector) 9 , 19 524

n

n

n n n n n

n n n n n n

n N

hy y f f f f

hy y f x y f f f

+

+

− − −

+ + − −

= −

= + − + −

= + + − +

(3.15)

O termo do erro de truncatura local do preditor e do corrector de (3.15) são de ordem

( )5O h .

3.3.1 Considerações

O corrector de (3.15) usa a aproximação ( )( )*

1 1 1,

n n nf f x y+ + +≈ no cálculo de 1n

y + . Como

1ny + é também uma aproximação de ( )1ny x + , ele pode ser usado novamente no corrector

de (3.15) para obter uma nova aproximação para 1n

f + , que por sua vez gera uma nova

aproximação para 1ny + . Na prática, invocando o corretor uma ou duas vezes é em geral


suficiente para produzir uma boa precisão, desde que se utilize um preditor da mesma

ordem. (PINA, 1995)

Sendo esse trabalho um trabalho limitado (em relação ao número de folhas), nós

não vamos fazer a análise da consistência e convergência dos métodos de passo

múltiplo, pois são análises muito complexas e longas, indicando para isso [QUO07] ou

[PIN95], caso o leitor estiver interessado.

Exemplo 3.5 Resolva pelos vários métodos apresentados o PVI dos exemplos (2.5) e (2.6)

com espaçamento 0.125h = .

A solução exacta do PVI ( ) ( ) [ ]2, , com 0 1, em 0,3f x y x y y I= − = = é

( ) 22 2

xy x e x x

−= − + − + .

A tabela 3.1 apresenta as soluções do PVI do exemplo acima utilizando os métodos de

Euler (explícito e Implícito), Heun, RK4, AB, e o par preditor corrector ABM.

A resolução desse exemplo foi feita recorrendo aos Programas 1 a 6 do anexo 2).

nx .0 125h = ( ) Exacta

ny x

Euler EXP Euler IMP Heun RK4 AB PC-ABM

.

.

.

.

.

.

.

.

.

.

.

0 000

0 125

0 25

0 375

0 50

0 75

1 00

1 50

2 00

2 50

3 00

.

.

.

.

.

.

.

.

.

.

.

1 0000

0 8750

0 7676

0 6794

0 6121

0 5448

0 5743

0 9488

1 7717

3 0644

4 8395

.

.

.

.

.

.

.

.

.

.

.

1 0000

0 8926

0 8025

0 7313

0 6804

0 6439

0 7012

1 1216

1 9788

3 2963

5 0887

.

.

.

.

.

.

.

.

.

.

.

1 0000

0 8838

0 7850

0 7052

0 6459

0 5935

0 6363

1 0324

1 8711

3 1750

4 9577

.

.

.

.

.

.

.

.

.

.

.

1 000000

0 883128

0 783699

0 703336

0 643470

0 590135

0 632123

1 026873

1 864668

3 167919

4 950217

.

.

.

.

.

.

.

.

.

.

.

1 000000

0 883128

0 783699

0 703336

0 643461

0 590115

0 632097

1 026844

1 864642

3 167897

4 950199

.

.

.

.

.

.

.

.

.

.

.

1 000000

0 883128

0 783699

0 703336

0 643471

0 590136

0 632124

1 026873

1 864667

3 167917

4 950214

.

.

.

.

.

.

.

.

.

.

.

1 000000

0 883128

0 783699

0 703335

0 643469

0 590133

0 632120

1 026869

1 864664

3 167915

4 950212



Da tabela acima, vemos que os métodos que produzem resultados mais próximos da

solução exacta são os métodos de Runge Kutta de 4ª ordem, os métodos de Adams

Bashforth e o preditor corrector ABM. Destes três métodos, os métodos de RK4 e ABM

produzem resultados melhores. A partir da tabela acima, também podemos ver que

utilizando o preditor AB sem a companhia do corrector AM obtém-se aproximações

menos interessantes do que se tivesse sido utilizado o corrector.

Ainda do exemplo acima podemos ver que o método RK4 tem precisão comparável

com o preditor corrector ABM, mas esse método (RK4) tem a desvantagem de requerer

um maior número de cálculos de ( , )f x y por passo, do que os métodos preditores

correctores e por isso são aconselháveis para funções mais simples de computar.

Para funções muito complexas, é aconselhável utilizar o método ABM. Por exemplo,

na tabela 3.2 são apresentados as soluções do PVI ' 1y

y xx

= − + com (1) 0y = com

[ ]1,4x ∈ cuja solução é ( ) ( )( )ln 1y x x x x= − + , utilizando os métodos RK4 e ABM com

0.2h = .

Observando a Tabela 3.2, vê-se que utilizando o preditor corrector de Adams

Bashforth Moulton obtém-se resultados sensivelmente mais próximas da solução exacta,

do que utilizando o método de Runge Kutta de 4ª ordem. Esse melhor resultado deve-se

ao facto de que no método RK4, efectuamos 4 cálculos de f por passo, o que faz

aumentar o erro de arredondamento para funções mais complexas (que é o caso do PVI

dado), enquanto para o preditor corrector ABM efectuamos apenas um cálculo da função

f por passo (visto que os restantes cálculos são aproveitados dos passos anteriores).


Tabela 3.2 Soluções pelos métodos ABM e RK4

do PVI ' 1= − +y

y xx

com (1) 0=y .

nx .h = 0 5 ( )ny x Exacta

ABM RK4

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

1 0

1 2

1 4

1 6

1 8

2 0

2 2

2 4

2 6

2 8

3 0

3 2

3 4

3 6

3 8

4 0

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

0 0000000

0 021212

0 088934

0 207987

0 381977

0 613700

0 905389

1 258871

1 675667

2 157062

2 704160

3 317915

3 999161

4 748636

5 566994

6 454821

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

0 0000000

0 021212

0 088934

0 207987

0 381975

0 613694

0 905381

1 258860

1 675654

2 157047

2 704143

3 317896

3 999140

4 748613

5 566970

6 454795

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

0 0000000

0 021214

0 088938

0 207994

0 381984

0 613705

0 905393

1 258875

1 675670

2 157065

2 704163

3 317917

3 999163

4 748638

5 566995

6 454822


CONCLUSÃO

Depois da elaboração deste trabalho provou-se que os métodos numéricos apresentam-

se como alternativa plausível para a resolução de PVI’s, na medida que oferecem a

oportunidade de obter soluções aproximadas muito próximos da solução exacta

(solução essa que nem sempre é possível obter-se pelos métodos analíticos).

Neste trabalho abordamos dois conjuntos de métodos numéricos para a resolução de

PVI’s: os métodos de passo simples (que utilizam a informação do passo anterior para

computar o passo seguinte) e os métodos de passo múltiplo (que utilizam as

informações dos vários passos computados anteriormente para computar o passo

seguinte). Dentre os métodos de passo simples, observa-se que os métodos de Euler são

os métodos que produzem piores resultados, pelo facto de utilizarem somente o termo

de grau 1 no desenvolvimento em série de Taylor da solução exacta. Por contrapartida

os métodos de Runge Kutta (em particular o método de Runge Kutta de 4ª ordem)

produzem melhores resultados, pois utilizam mais termos no desenvolvimento em série

de Taylor da solução exacta, embora possuam a desvantagem de necessitarem de vários

cálculos de f por passo.

Em relação aos métodos de passo múltiplo, observa-se que a utilização de um

método de Adams Bashforth (como por exemplo o método de Adams Bashforth de 4ª

ordem, que foi o exemplo apresentado) “sozinho”(i.e., não acompanhado dum método de


Adams Moulton para servir de corrector) produzem resultados não muito plausíveis

(por exemplo comparando com RK4), enquanto que quando utilizados em conjunto com

os métodos de Adams Moulton (da mesma ordem) formando assim um par preditor

corrector, produzem resultados muito próximos da solução exacta.

Vale também frisar o contacto com o software MATLAB como ferramenta auxiliar,

que foi deveras gratificante dada a sua potencialidade no tratamento do conteúdo do

tema em questão.

Com a utilização do software MATLAB, foi possível evitar cálculos de difíceis

execuções (caso usássemos apenas lápis e papel) inerentes aos vários métodos

apresentados (principalmente os métodos de Runge Kutta de 4ª ordem), obtendo os

resultados de forma rápida e eficiente.

Finalmente, os objectivos preconizados foram atingidos, aprofundei conhecimentos

não só do tema em questão, mas de diferentes áreas da matemática, que é um dos

propósitos de uma monografia. Infelizmente o tema é muito vasto, daí a impossibilidade

de abordar e confrontar outros métodos, também interessantes, que podem ser

encontradas nas várias bibliografias utilizadas para a confecção desta monografia, tais

como [BUT08], [QUO07] entre outras.


Referências Bibliográficas

AGUDO, F. D. (1992). Análise Real (Vol. III). Escolar Editora.

ASCHER, U. M., & PETZOLD, L. R. (1998). Computer Methods for Ordinary Differential

Equations and Differential-Algebraic Equations. Philadelphia: Siam.

BUTCHER, J. C. (2008). Numerical Methods for Ordinary Differential Equations

(Second Edition ed.). John Wiley & Sons Ltd.

KINCARD, D. (2002). Numerical Analysis: Mathematics of Scientific Computing. Word

Cheney Brooks/Cole.

LINZ, P., & WANG, R. L. (2003). Numerical Methods. JONES AND BARTLETT

MATHEMATICS.

MATHEWS, J. H., & FINK, K. D. (2004). Numerical Methods Usin MATLAB (Fourth

Edition ed.). International Edition.

MILLER, R. K. (1991). Introduction to Differential Equations (second edition ed.). New Jersey: Prentice Hall.

MORAIS, V., & VIEIRA, C. (2006). MATLAB 7 & 6 Curso Completo. FCA.

PINA, H. (1995). Métodos Numéricos. Mc Graw-Hill.

QUORTERONI, A., SACCO, R., & SALERI, F. (2007). Texts in Applied Mathematics.

Numerical Mathematics (2ª ed.). Springer.

RAINER, K. (1998). Numerical Analysis. Springer.

RILEY, K. F., HOBSON, M. P., & BENCE, S. J. (2007). Mathematical Methods for Physics

and Engineering (3ª ed.). Cambridge.


Bibliografia

[AGU92] AGUDO, F. D. (1992). Análise Real (Vol. III). Escolar Editora.

[ASC98] ASCHER, U. M., & PETZOLD, L. R. (1998). Computer Methods for Ordinary Differential

Equations and Differential-Algebraic Equations. Philadelphia: Siam.

[BUT08] BUTCHER, J. C. (2008). Numerical Methods for Ordinary Differential Equations (Second Edition ed.). John Wiley & Sons Ltd.

[CAR09] CARVALHO, N. T. (2009). Aproximações Numéricas e Aplicações com MATLAB

de:Equações não Lineares & Equações Diferenciais. Praia: Monografia (Licenciatura em Matemática)- Universidade de cabo verde.

[KIN02] KINCARD, D. (2002). Numerical Analysis: Mathematics of Scientific Computing. Word Cheney Brooks/Cole.

[LIN03] LINZ, P., & WANG, R. L. (2003). Numerical Methods. JONES AND BARTLETT MATHEMATICS.

[MAT04] MATHEWS, J. H., & FINK, K. D. (2004). Numerical Methods Usin MATLAB (Fourth Edition ed.). International Edition.

[MIL91] MILLER, R. K. (1991). Introduction to Differential Equations (second edition ed.). New Jersey: Prentice Hall.

[MOR06] MORAIS, V., & VIEIRA, C. (2006). MATLAB 7 & 6 Curso Completo. FCA.

[PIN95] PINA, H. (1995). Métodos Numéricos. Mc Graw-Hill.

[QUO07] QUORTERONI, A., SACCO, R., & SALERI, F. (2007). Texts in Applied Mathematics.

Numerical Mathematics (2ª ed.). Springer.

[RAI98] RAINER, K. (1998). Numerical Analysis. Springer.

[Spr09] Springer Online Reference Works. Cauchy problem, numerical methods for ordinary

differential equations. Disponível em < http://eom.springer.de/C/c020960.htm >. Acesso em 15 de Setembro de 2009.

[RIL07] RILEY, K. F., HOBSON, M. P., & BENCE, S. J. (2007). Mathematical Methods for Physics

and Engineering (3ª ed.). Cambridge.


ANEXO I— INTRODUÇÃO AO MATLAB

O MATLAB (MATrix LABoratory) é um poderoso software matemático de computação numérica, de análise e visualização de dados que trabalha a base de matrizes. Sua vantagem está na manipulação e cálculos matriciais, como por exemplo, resolução de sistemas lineares, cálculos de autovalores e autovectores, resolução de equações diferenciais analiticamente e numericamente, entre outros. Também este software dispõe da possibilidade de visualização e manipulação de gráficos de 2D e 3D bem como a possibilidade de implementação de nossos próprios programas, para além dos que já

vêm implementados.

Figura 1. Janela inicial do MATLAB

Os comandos do MATLAB são muito próximo da forma como escrevemos expressões algébricas, tornando mais simples o seu uso. Apresentamos em seguida alguns

comandos básicos do MATLAB.

Operadores aritméticos

+ Adição - Subtracção * Multiplicação / Divisão ^ Potenciação pi, e Constantes (π, Euler)


Funções abs(#) Valor absoluto cos(#), sin(#),tan(#) Funções trigonométricas( seno, cosseno, tangente) exp(#) Exponencial (e(#)) log(#) Logaritmo na base e (ln(#)) log10(#) Logaritmo na base 10 sqrt(#) Raiz quadrada Operadores Lógicos and “e” (&) not negação lógica(~) or “ou” (|) xor “ou exclusivo” Operadores relacionais == Igual ~= Diferente < Menor <= Menor ou igual > Maior >= Maior ou igual Valores Booleanos 1 Verdadeiro 0 Falso

Definindo Funções

No software MATLAB o usuário pode definir sua própria função, construindo um M-file (um ficheiro de formato .m) no M-file Debugger (Fig 2). Uma vez definida a função, o usuário pode “chamar” a sua função da mesma maneira que se chamam as funções

(abs(#), sin(#), cos(#), etc...).


Figura 2. Janela do Editor Debugger

Exemplo

Vamos implementar a função π= 2Area r r( ) no M-file que vamos nomear por area.m. No

editor/Debugger vamos escrever o seguinte

Figura 3. Escrevendo a função area

Uma vez guardada essa função como um M-file chamado area.m, ela pode ser

chamado no comando Window do MATLAB, como mostra a figura seguinte.


Figura 4. Chamando a função area para r=3.

O M-file do exemplo acima foi guardado no ambiente de trabalho, por isso, ele deve ser

chamado do ambiente de trabalho (ver parte entre linhas negras na figura acima).

Controle de fluxo num programa

O controlo de fluxo num programa escrito no MATLAB é muito semelhante a aqueles

usados noutras linguagens de programação.

for Repete instruções um número específico de vezes; while Repete instruções um número indefinido de vezes; if Execução condicional de instruções; switch Escolhe uma opção, entre várias, através da avaliação de uma expressão; end Termina estruturas for, while, switch, e if;

break Termina execução do ciclo for ou while;

Como exemplo vamos implementar um programa que calcula o factorial de um número dado.

function F=factorial(n) % F função que calcula o factorial de um número % n-número>0 %inicialização F=1; %Ciclo for na forma inicio:fim for i=1:n F=F*i; end


O símbolo “%”, é usado quando queremos fazer comentários durante a implementação dum programa. Inserindo o símbolo “%”numa linha de comandos, tudo o que

escrevermos a seguir não será incluído no algoritmo.

Gráficos

Com o software MATLAB pode-se construir gráficos de 2 e 3 dimensões de curvas e superfícies. O comando plot é utilizado para construir gráficos de funções 2D. O seguinte exemplo

mostra como construir os gráficos das funções ( )cosy x= e 2z x= no intervalo [ ]0,π .

Exemplo >>x=0:0.1:pi; >> y=cos(x); >> z=x.^2; >> plot(x,y,x,z,'+')

Na primeira linha do exemplo especificou-se o domínio com espaçamento 0.1. Nas duas linhas seguintes estão definidas as funções. Notemos que as três primeiras linhas terminam com um “ponto e vírgula”. O “ponto e vírgula” é utilizado para evitar que as matrizes x, y e z sejam apresentados na janela de comandos logo após serem escritos. A quarta linha contém o comando plot que executa o gráfico. Os dois primeiros termos no comando plot (x e y) produzem o gráfico da função y=cos(x), a terceira e quarta linha(x e z) produzem o gráfico de z=x2 ( na construção do gráfico de uma função, quando quisermos elevar a variável a um expoente temos de colocar um ponto “.” Antes de introduzir o símbolo “^”). O último termo “+” serve para apresentar a função z como sucessivos “+”, como mostra a figura seguinte.

Figura 5


Os gráficos 3D são obtidos pela especificação de um rectângulo do domínio da função através do comando mesgrid e utilizando os comandos mesh ou surf para obter o

gráfico.

Exemplo

>> x=-pi:0.1:pi; >> y=x; >> [x,y]=meshgrid(x,y); >> z=sin(cos(x+y)); >> mesh(z) Cujo gráfico será

Figura 6

MATLAB & PVI’s

Como foi dito anteriormente o software MATLAB dispõe de vários programas pré implementados. Em relação aos PVI´s, o software possui alguns programas tais como: ode45 (Runge Kutta de 4ª e 5ª ordem), ode23 (Runge Kutta de 2ª e 3ª ordem) e o ode113

(Adams)

As funções ode23, ode45 e ode113 apresentam a seguinte sintaxe

[T,Y] = solver(odefun,tspan,y0), onde solver pode ser ode23, ode45 ou ode113

Na função acima o argumento odefun é a nossa função ( ),y f t y′ = ,

tspan [ ]0 ,t T= corresponde ao intervalo de integração e y0 representa a condição inicial

do PVI.


Considerações

Apresentamos essa pequena introdução do MATLAB, apenas para o leitor ter uma pequena ideia deste poderoso software, mas para uma melhor compreensão indicamos [MOR06], ou de uma forma mais prática, que procure a ajuda que vem incorporada no programa. O leitor também pode aceder ao site da MATH WORKS (que é o site oficial do software MATLAB) http://www.mathworks.com para mais esclarecimentos.


ANEXO II---Programas Implementados e utilizados para a apresentação dos exemplos.

PROGRAMA 1 %Método de Euler explícito (ou progressivo) para a aproximação %numérica da solução do PVI y´=f(x,y), com x pertencente %ao intervalo [x0, X], com a condição inicial y(x0)=y0. % y(n+1)=y(n)+h*f(x(n),y(n)), n=0,1,2,...,N-1 function E=eulerexplicito(f,x0,X,y0,N) % Entradas % f- função f(x,y) (tipo string); % x0- limite lateral esquerdo do intervalo; %X- limete lateral direito do intervalo % y0- condição inicial; % N- número de passos; % Saída %E=[x y] onde x corresponde ao vector das abcissas e % y corresponde ao vector das ordenadas. %definir h h=((X-x0)/N); % Definir vectores coluna de y e x x=zeros(1,N+1); y=zeros(1,N+1); y(1)=y0; %definir os nós x=x0:h:X; %cálculo de y for n=1:N y(n+1)=y(n)+h*feval(f,x(n),y(n)); end E=[x' y'];


PROGRAMA 2 %Método de Euler implícito (ou regressivo) para a aproximação %numérica da solução do PVI y´=f(x,y), com x pertencente %ao intervalo [x0,X], com a condição inicial y(x0)=y0. % y(n+1)=y(n)+h*[f(x(n+1),y(n)+h*f(x(n),y(n))], n=0,1,2,...,N-1 %0nde o método de euler explicito é utilizado para a obtenção da primeira %aproximação de y(n+1) (pois estamos na presença dum método implícito). function I=eulerimplicito(f,x0,X,y0,N) % Entradas % f- função f(x,y) (tipo string); % x0- limite lateral esquerdo do intervalo; %X-limite lateral direito do intervalo % y0- condição inicial; % N- número de passos; % Saída %I=[x y] onde x correspode ao vector das abcissas e % y corresponde ao vector das ordenadas. %definir h h=((X-x0)/N); % Definir vectores coluna de y e x x=zeros(1,N+1); y=zeros(1,N+1); y(1)=y0; %definir os nós x=x0:h:X; %cálculo de y for n=1:N y(n+1)=y(n)+h*feval(f,x(n+1),y(n)+h*feval(f,x(n),y(n))); end I=[x' y'];


PROGRAMA 3 %Método de Heun para a aproximação %numérica da solução do PVI y´=f(x,y), com x pertencente ao %ao intervalo [x0,X], com a condição inicial y(x0)=y0. % y(n+1)=y(n)+(h/2)*[f(x(n),y(n))+f(x(n+1),y(n+1))], n=0,1,2,...,N-1 %0nde o método de euler explicito é utilizado para a obtenção da primeira %aproximação de y(n+1) (pois estamos na presença dum método implícito). function H=heun(f,x0,X,y0,N) % Entradas % f- função f(x,y) (tipo string); % x0- limite lateral esquerdo do intervalo; %X-limite lateral direito o intervalo % y0- condição inicial; % N- número de passos; % Saída %H=[x y] onde x corresponde ao vector das abcissas e % y corresponde ao vector das ordenadas. %definir h h=((X-x0)/N); % Definir vectores coluna de y e x x=zeros(1,N+1); y=zeros(1,N+1); y(1)=y0; %definir os nós x=x0:h:X; %cálculo de y for n=1:N F1=feval(f,x(n),y(n)); F2=feval(f,x(n+1),y(n)+h*F1); y(n+1)=y(n)+(h/2)*(F1+F2); end H=[x' y'];


PROGRAMA 4 %Método de Runge Kutta de 4ª ordem para a aproximação %numérica da solução do PVI y´=f(x,y), com x pertencente %ao intervalo [x0,X], com a condição inicial y(x0)=y0. % y(n+1)=y(n)+(h/6)*(F1+2F2+2F3+F4), n=0,1,2,...,N-1 %onde %F1=f(x(n),y(n)); %F2=f(x(n)+(h/2),y(n)+(hF1/2)); %F3=f(x(n)+(h/2),y(n)+(hF2/2)); %F4=f(x(n)+h, y(n)+h*F3) function R=rk4(f,x0,X,y0,N) % Entradas % f- função f(x,y) (tipo string); % x0- limite lateral esquerdo do intervalo; %X-limite lateral direito o intervalo % y0- condição inicial; % N- número de passos; % Saída %R=[x y] onde x corresponde ao vector das abcissas e % y corresponde ao vector das ordenadas. %definir h h=((X-x0)/N); % Definir vectores coluna de y e x x=zeros(1,N+1); y=zeros(1,N+1); y(1)=y0; %definir os nós x=x0:h:X; %cálculo de y for n=1:N F1=feval(f,x(n),y(n)); F2=feval(f,x(n)+(h/2),y(n)+(h*F1/2)); F3=feval(f,x(n)+(h/2),y(n)+(h*F2/2)); F4=feval(f,x(n)+h, y(n)+h*F3); y(n+1)=y(n)+(h/6)*(F1+2*F2+2*F3+F4); end R=[x' y'];


PROGRAMA 5 %Processo preditor corrector de Adams-Bashforth-Moulton (de quarta ordem) para a aproximação %numérica da solução do PVI y´=f(x,y), com x pertencente %ao intervalo [x0,X], com a condição inicial y(x0)=y0. function A=abm(f,x0,X,y0,N) % Entradas % f- função f(x,y) (tipo string); % x0- limite lateral esquerdo do intervalo; % X-limite lateral direito o intervalo % y0- condição inicial; % N- número de passos; % Saída %A=[x y] onde x corresponde ao vector das abcissas e % y corresponde ao vector das ordenadas. %definir h h=((X-x0)/N); % Definir vectores coluna de y e x x=zeros(1,N+1); y=zeros(1,N+1); y(1)=y0; %definir os nós x=x0:h:X; %cálculo de y % Lembrar que os 4 primeiros coordenadas de x e y devem ser obtidos por RK4 for n=1:3 F1=feval(f,x(n),y(n)); F2=feval(f,x(n)+(h/2),y(n)+(h*F1/2)); F3=feval(f,x(n)+(h/2),y(n)+(h*F2/2)); F4=feval(f,x(n)+h, y(n)+h*F3); y(n+1)=y(n)+(h/6)*(F1+2*F2+2*F3+F4); for n=4:N %preditor G1=feval(f,x(n),y(n)); G2=feval(f,x(n-1),y(n-1)); G3=feval(f,x(n-2),y(n-2)); G4=feval(f,x(n-3),y(n-3)); p=y(n)+(h/24)*(55*G1-59*G2+37*G3-9*G4);


%corrector y(n+1)=y(n)+(h/24)*(9*feval(f,x(n+1),p)+19*G1-5*G2+G3); end end A=[x' y']; PROGRAMA 6 %Método de Adams-Bashforth (de quarta ordem) para a aproximação %numérica da solução do PVI y´=f(x,y), com x pertencente ao % intervalo [x0,X], com a condição inicial y(x0)=y0. function B=ab(f,x0,X,y0,N) % Entradas % f- função f(x,y) (tipo string); % x0- limite lateral esquerdo do intervalo; %X- limite lateral direito o intervalo % y0- condição inicial; % N- número de passos; % Saída %B=[x y] onde x corresponde ao vector das abcissas e % y corresponde ao vector das ordenadas. %definir h h=((X-x0)/N); % Definir vectores coluna de y e x x=zeros(1,N+1); y=zeros(1,N+1); y(1)=y0; %definir os nós x=x0:h:X; %cálculo de y % Lembrar que os 4 primeiros coordenadas de x e y devem ser obtidos por RK4 for n=1:3 F1=feval(f,x(n),y(n)); F2=feval(f,x(n)+(h/2),y(n)+(h*F1/2));


F3=feval(f,x(n)+(h/2),y(n)+(h*F2/2)); F4=feval(f,x(n)+h, y(n)+h*F3); y(n+1)=y(n)+(h/6)*(F1+2*F2+2*F3+F4); for n=4:N %preditor G1=feval(f,x(n),y(n)); G2=feval(f,x(n-1),y(n-1)); G3=feval(f,x(n-2),y(n-2)); G4=feval(f,x(n-3),y(n-3)); y(n+1)=y(n)+(h/24)*(55*G1-59*G2+37*G3-9*G4); end end B=[x' y'];

Documents

Métodos Numéricos na resolução de Equações Diferenciais Ordinárias (EDO)