[Cuminato] Discretização de Equações Diferenciais Parciais

Prévia do material em texto

DISCRETIZAC¸A˜O DE EQUAC¸O˜ES
DIFERENCIAIS PARCIAIS
Te´cnicas de Diferenc¸as Finitas
Jose´ Alberto Cuminato
Instituto de Cieˆncias Matema´ticas e de Computac¸a˜o
Universidade de Sa˜o Paulo – Campus de Sa˜o Carlos
Messias Meneguette Junior
Faculdade de Cieˆncias e Tecnologia
Universidade Estadual Paulista – Presidente Prudente
13 de outubro de 2006
2
Suma´rio
1 Conceitos Fundamentais 7
1.1 Introduc¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Aproximac¸o˜es de Derivadas por Diferenc¸as Finitas . . . . . . . . . . . . . . . 9
1.3 Problema de Valor Inicial em Equac¸o˜es Ordina´rias . . . . . . . . . . . . . . . 15
1.3.1 Convergeˆncia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4 Problema de Valor de Fronteira em Equac¸o˜es Ordina´rias . . . . . . . . . . . . 29
1.4.1 Me´todos de Diferenc¸as Finitas . . . . . . . . . . . . . . . . . . . . . . 30
1.5 Exerc´ıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2 Introduc¸a˜o a`s Equac¸o˜es Diferenciais Parciais 45
2.1 Introduc¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.2 Classificac¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.3 Caracter´ısticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.4 Discretizac¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.5 Domı´nios Gene´ricos e Transformac¸a˜o de Varia´veis . . . . . . . . . . . . . . . . 53
2.5.1 Transformac¸a˜o de Coordenadas Cartesianas em Polares . . . . . . . . 55
2.5.2 Transformac¸a˜o de Coordenadas Cartesianas em Obl´ıquas . . . . . . . 57
2.5.3 Transformac¸a˜o de Coordenadas Cartesianas em Triangulares . . . . . 60
2.6 Exerc´ıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3 Equac¸o˜es Parabo´licas 69
3.1 Introduc¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.2 Me´todos de Diferenc¸as Finitas . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.3 Problemas Na˜o Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
3.4 Outros Me´todos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
3.5 Equac¸o˜es Parabo´licas em Duas Dimenso˜es . . . . . . . . . . . . . . . . . . . . 97
3.6 Exerc´ıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4 Equac¸o˜es El´ıpticas 113
4.1 Introduc¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.2 Me´todos de Diferenc¸as Finitas . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
4.3 Erro de Truncamento Local . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
4.4 Condic¸o˜es de Fronteira em Domı´nios Gerais . . . . . . . . . . . . . . . . . . . 123
4.5 Condic¸a˜o de Fronteria de Neumann . . . . . . . . . . . . . . . . . . . . . . . . 128
3
4 SUMA´RIO
4.6 Diferenc¸as Finitas em Coordenadas Polares . . . . . . . . . . . . . . . . . . . 130
4.7 Me´todos Iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
4.8 Convergeˆncia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
4.9 Me´todo dos Gradientes Conjugados . . . . . . . . . . . . . . . . . . . . . . . . 143
4.10 Exerc´ıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
5 Equac¸o˜es Hiperbo´licas 157
5.1 Introduc¸a˜o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
5.1.1 Equac¸a˜o escalar de advecc¸a˜o . . . . . . . . . . . . . . . . . . . . . . . 157
5.1.2 Equac¸a˜o da onda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
5.2 Diferenc¸as Finitas para a Equac¸a˜o de Advecc¸a˜o . . . . . . . . . . . . . . . . 190
5.3 Diferenc¸as Finitas para a Equac¸a˜o da Onda . . . . . . . . . . . . . . . . . . . 206
5.4 Me´todos Nume´ricos para Leis de Conservac¸a˜o . . . . . . . . . . . . . . . . . . 211
5.5 Exerc´ıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 234
SUMA´RIO 5
PREFA´CIO
O universo de te´cnicas de soluc¸a˜o nume´rica de equac¸o˜es diferenciais parciais e´ bastante vari-
ado e inclue, entre outros, os me´todos de diferenc¸as finitas, volumes finitos, elementos finitos,
elementos de contorno, me´todos espectrais e me´todos de colocac¸a˜o. O objetivo deste texto
e´ apresentar uma introduc¸a˜o a` soluc¸a˜o nume´rica de equac¸o˜es diferenciais parciais, enfati-
zando os principais conceitos com vistas a`s aplicac¸o˜es pra´ticas, utilizando discretizac¸a˜o por
diferenc¸as finitas. O conteu´do sera´ desenvolvido atrave´s de exemplos-modelos de equac¸o˜es
parabo´licas, el´ıpticas e hiperbo´licas. Pretende-se introduzir de maneira bastante natural os
conceitos de estabilidade e convergeˆncia, dependeˆncia da discretizac¸a˜o das condic¸o˜es de fron-
teira, condic¸a˜o CFL e me´todo das caracter´ısticas. Incluimos, sempre que poss´ıvel, uma certa
dose de teoria e rigor no tratamento desses conceitos sem no entanto relegarmos o lado ex-
perimental, especialmente o lado ilustrativo atrave´s de exemplos e gra´ficos. Nosso objetivo e´
o de apresentar uma variedade de te´cnicas de soluc¸a˜o sobre as quais o leitor devera´ adquirir
conhecimentos teo´ricos e tambe´m experimental. Por essa raza˜o e´ imprescind´ıvel que o leitor
resolva os exerc´ıcios pra´ticos, que envolvem programac¸a˜o em MATLAB.
Na literatura e´ bastante popular a ide´ia, especialmente em equac¸o˜es parciais, de que
existem tantos me´todos nume´ricos quantos sa˜o os problemas a serem resolvidos. Em geral
entretanto, alguns conceitos e te´cnicas fundamentais esta˜o presentes em todos os me´todos.
O propo´sito deste texto e´ o de chamar a atenc¸a˜o para esses conceitos esclarecendo-os atrave´s
da ana´lise de exemplos. E´ quase sempre poss´ıvel classificar um problema em quatro grandes
tipos: Equac¸o˜es Parabo´licas, Equac¸o˜es El´ıpticas, Equac¸o˜es Hiperbo´licas e Problemas do
tipo misto. Do ponto de vista dida´tico esta divisa˜o parece ser recomendada para melhor
associac¸a˜o entre as te´cnicas nume´ricas e a natureza do problema a ser resolvido. O conteu´do
aqui desenvolvido segue essa lo´gica.
Este livro foi escrito com o propo´sito de servir como texto ba´sico de disciplinas de
ca´lculo nume´rico e matema´tica aplicada ministradas em cursos de graduac¸a˜o de engenharia,
matema´tica e a´reas afins em cieˆncias exatas. O material mais avanc¸ado pode ainda servir
para cursos introduto´rios de po´s-graduac¸a˜o nas a´reas acima citadas. Para cumprir esse ob-
jetivo foi introduzido uma lista bastante completa de exerc´ıcios que envolve desde a fixac¸a˜o
dos conceitos ba´sicos ate´ a implementac¸a˜o de projetos pra´ticos em computador. A maioria
dos me´todos nume´ricos apresentados foram programados em MATLAB e podem ser utiliza-
dos pelo professor para ilustrar os conceitos ministrados em sala de aula. A programac¸a˜o
em MATLAB inclue uma interface gra´fica padronizada que permite a soluc¸a˜o de problemas
razoavelmente complexos sem a necessidade de um grande esforc¸o de programac¸a˜o por parte
do aluno. Isto viabiliza a utilizac¸a˜o de exerc´ıcios pra´ticos para a ilustrac¸a˜o e reconhecimento
de conceitos como convergeˆncia, estabilidade, precisa˜o, etc.
O conteu´do foi dividido em 5 cap´ıtulos, os dois primeiros apresentando os conceitos
ba´sicos de discretizac¸a˜o de problemas de valor inicial e de contorno para equac¸o˜es diferenciais
ordina´rias e a teoria elementar das equac¸o˜es diferenciais parciais. Nos outros treˆs cap´ıtulos
apresentamos separadamente as te´cnicas de discretizac¸a˜o de cada um dos tipos de equac¸o˜es
parabo´licas, hiperbo´licas e el´ıpticas. Todos os cap´ıtulos inclue uma lista bastante extensa
de exerc´ıcios, muitos dos quais abordam temas mais avanc¸adosdo que aqueles tratados no
texto. Essa atitude tem o propo´sito de estimular o aluno que pretende aprofundar-se nessa
a´rea do conhecimento a pesquisar as refereˆnciais e aguc¸ar sua curiosidade. No entanto,
6 SUMA´RIO
alertamos os professores que pretendem utilizar o texto como guia de um curso, para o fato
de que uma selec¸a˜o dos exerc´ıcios pode ser recomendada para na˜o desestimular os estudantes
com problemas muito dif´ıceis. A ordem de apresentac¸a˜o dos exerc´ıcios na˜o segue um crite´rio
de dificuldade crescente e sim a ordem dos to´picos no texto.
O material deste livro deriva em boa parte de cursos de graduac¸a˜o que os autores teˆm
ministrado em suas respectivas unidades de ensino. Agradecimento especial nos e´ devido aos
autores da apostila [28], aos alunos do ICMC pela leitura dos originais e o apoio financeiro
do Instituto de Cieˆncias Matema´ticas e de Computac¸a˜o - USP Sa˜o Carlos, Faculdade de
Cieˆncias e Tecnologia - UNESP Presidente Prudente, CNPq e Pro´-Reitoria de Graduac¸a˜o
da USP atrave´s dos Programa SIAE 97 e SIAE 98.
Cap´ıtulo 1
Conceitos Fundamentais
1.1 Introduc¸a˜o
Um modelo matema´tico e´ uma representac¸a˜o, idealizada e muitas vezes simplificada,
da natureza. Quando derivado de maneira criteriosa, sua soluc¸a˜o simula propriedades dos
processos naturais envolvidos, tais como velocidade e pressa˜o no escoamento de um fluido,
temperatura, humidade, direc¸a˜o dos ventos e precipitac¸a˜o na previsa˜o do tempo, trajeto´ria
de um sate´lite artificial, entre outras muitas aplicac¸o˜es. Assim, as soluc¸o˜es das equac¸o˜es de
um modelo, devem captar os aspectos relevantes no comportamento do problema modelado,
na˜o sendo poss´ıvel, na maioria dos casos, justificar a utilizac¸a˜o de hipo´teses simplificado-
ras que modificam a natureza do problema (tal como linearidade) que tornariam poss´ıvel
a determinac¸a˜o de uma soluc¸a˜o exata. Em outra palavras, quando estamos resolvendo
um problema real, em geral, na˜o e´ poss´ıvel forc¸a´-lo a satisfazer hipo´teses que permitam a
obtenc¸a˜o de uma soluc¸a˜o exata. Da´ı a necessidade da procura de soluc¸o˜es nume´ricas, ou
aproximadas.
A importaˆncia da modelagem matema´tica cresceu muito nos u´ltimos anos porque a
ana´lise detalhada de um modelo e suas propriedades, viabilizada pelo avanc¸o tecnolo´gico
dos computadores, permite um melhor entendimento do evento modelado e, mais do que
isso, permite a simulac¸a˜o de mudanc¸as nos paraˆmetros do modelo e a ana´lise da respec-
tiva resposta, que na˜o seriam poss´ıveis na situac¸a˜o real. Por exemplo, no projeto de um
carro, alterac¸o˜es na forma resultam em modificac¸o˜es no comportamento aerodinaˆmico, de
cuja ana´lise obte´m-se um indicador dos ganhos e/ou perdas em performance; a localizac¸a˜o
ideal da asa de um avia˜o em relac¸a˜o ao casco pode ser obtida pela resposta a` simulac¸a˜o
matema´tica das equac¸o˜es da aerodinaˆmica, e ate´ a melhor pol´ıtica de vacinac¸a˜o contra
doenc¸as transmiss´ıveis, tipo sarampo, pode ser decidida com base em modelos matema´ticos.
A economia de tempo gerada por esta maneira de projetar um produto, ou tomar deciso˜es,
e´ clara, diminuindo sensivelmente o nu´mero de proto´tipos ou modelos em tamanho reduzido
a serem constru´ıdos e ensaiados. No projeto de equipamentos com tecnologia avanc¸ada, a
simulac¸a˜o atrave´s da modelagem matema´tica e´ uma ferramenta indispensa´vel a ser utilizada
na eliminac¸a˜o de casos triviais ou imposs´ıveis, fornecendo um guia seguro para a selec¸a˜o dos
casos a serem ensaiados em modelos de escala reduzida ou para a construc¸a˜o de proto´tipos.
Como virtualmente todas as a´reas em matema´tica aplicada, o principal objetivo das
7
8 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
equac¸e˜s diferenciais e´ o de modelar matematicamente eventos da natureza nas mais diversas
a´reas de aplicac¸a˜o. Um modelo com va´rias varia´veis independentes envolve necessariamente
derivadas com respeito a cada umas dessas varia´veis ou derivadas parciais, sendo portanto
denominado uma equac¸a˜o diferencial parcial ou um sistema de equac¸o˜es diferenciais parciais.
Resumindo as ide´ias apresentadas acima, o processo de simulac¸a˜o e´ considerado na lit-
eratura como sendo constitu´ıdo de treˆs partes nitidamente distintas. A fase de modelagem,
isto e´, a construc¸a˜o de um conjunto de equac¸o˜es matema´ticas que reputamos representar
os fenoˆmenos e os processos modelados. Uma segunda fase, de soluc¸a˜o desse conjunto de
equac¸o˜es, normalmente utilizando te´cnicas de discretizac¸a˜o nume´rica e um computador e
finalmente, a fase de interpretac¸a˜o dos resultados face a`s caracter´ısticas do problema origi-
nal. Esse e´ um processo complexo que exige do profissional um conjunto bastante amplo de
habilidades; como por exemplo: um bom conhecimento de engenharia, que sera´ necessa´rio
para os ajustes finos do modelo desprezando complicac¸o˜es que na˜o sa˜o fundamentais, um
bom conhecimento de me´todos nume´ricos para selecionar aquele que melhor se adapta ao
problema e finalmente um bom faro de detetive para analisar os resultados e interpreta´-los a`
luz das restric¸o˜es e caracter´ısticas do problema. Este texto trata exclusivamente da segunda
fase desse processo.
E´ claro portanto, que a soluc¸a˜o do modelo matema´tico, ou seja, das equac¸o˜es represen-
tantes desse modelo, e´ fundamental para a compreensa˜o do fenoˆmeno modelado. Este e´ o
papel da discretizac¸a˜o das equac¸o˜es parciais, uma vez que, como ja´ dissemos, uma soluc¸a˜o
anal´ıtica nem sempre esta´ dispon´ıvel: ou tem uma forma na˜o pra´tica ou e´ imposs´ıvel de
ser obtida. Assim os me´todos nume´ricos sa˜o amplamente usados para a aproximac¸a˜o dessas
soluc¸o˜es. A esseˆncia dos me´todos nume´ricos esta´ na representac¸a˜o discreta (finita) do prob-
lema que, em geral, e´ originalmente modelado como um cont´ınuo. Essa discretizac¸a˜o e´ que
viabiliza o uso dos computadores no tratamento nume´rico das equac¸o˜es diferenciais.
Na literatura, especialmente em equac¸o˜es parciais, e´ comum a ide´ia, de que existem tan-
tos me´todos nume´ricos quantos sa˜o os problemas a serem resolvidos. Em geral entretanto,
algumas te´cnicas e caracter´ısticas fundamentais esta˜o presentes em todos os me´todos. O
propo´sito deste texto e´ o de chamar a atenc¸a˜o para esses conceitos definindo-os precisamente
e ilustrando-os atrave´s da ana´lise de exemplos. E´ quase sempre poss´ıvel classificar um prob-
lema em equac¸o˜es diferenciais em quatro grandes classes: Equac¸o˜es Parabo´licas, Equac¸o˜es
El´ıpticas, Equac¸o˜es Hiperbo´licas e Problemas do tipo misto. Do ponto de vista dida´tico
esta divisa˜o parece ser recomendada para melhor associac¸a˜o entre as te´cnicas nume´ricas e a
natureza do problema a ser resolvido. Seguiremos neste texto essa premissa lo´gica.
O objetivo deste texto e´ o de desenvolver uma introduc¸a˜o a` soluc¸a˜o nume´rica de equac¸o˜es
diferenciais parciais, enfatizando os principais conceitos com vistas a`s aplicac¸o˜es pra´ticas,
atrave´s da discretizac¸a˜o por diferenc¸as finitas. O conteu´do sera´ desenvolvido por meio de
exemplos modelos de equac¸o˜es parabo´licas, el´ıpticas e hiperbo´licas. Pretende-se introduzir
de maneira bastante natural os conceitos de estabilidade e convergeˆncia, dependeˆncia da
discretizac¸a˜o das condic¸o˜es de fronteira, condic¸a˜o CFL e me´todo das caracter´ısticas.
Os conceitos fundamentais e as principais dificuldades, bem como a notac¸a˜o e nomen-
clatura ba´sicas sera˜o introduzidos o caso de problemas de uma dimensa˜o, ou seja, por meio
de equac¸o˜es diferenciais ordina´rias. Na sequeˆncia, sera˜o tratados cada um dos principais
tipos de equac¸o˜es parciais.
1.2. APROXIMAC¸O˜ES DE DERIVADAS POR DIFERENC¸AS FINITAS 9
1.2 Aproximac¸o˜es de Derivadas por Diferenc¸as Finitas
“A esseˆncia dos me´todos nume´ricos esta´ na discretizac¸a˜o do cont´ınuo. E´ esta
discretizac¸a˜o que torna finito o problemae portanto viabiliza sua soluc¸a˜o atrave´s
dos computadores [19].”
A ide´ia geral do me´todo de diferenc¸as finitas e´ a discretizac¸a˜o do domı´nio e a substi-
tuic¸a˜o das derivadas presentes na equac¸a˜o diferencial, por aproximac¸o˜es envolvendo somente
valores nume´ricos da func¸a˜o. Na pra´tica substitui-se as derivadas pela raza˜o incremental
que converge para o valor da derivada quando o incremento tende a zero. Dizemos enta˜o
que o problema foi discretizado. Quando o domı´nio tem mais de uma varia´vel, a ide´ia acima
e´ aplicada para cada uma das varia´veis separadamente.
Nesta sec¸a˜o trataremos apenas o problema unidimensional, pois isto facilita o entendi-
mento. No entanto, a generalizac¸a˜o para outras dimenso˜es pode ser obtida sem muitas
dificuldades, como veremos oportunamente. Seja x0 um nu´mero real qualquer, e h um
nu´mero positivo. Definimos malha de passo h associada a x0 como o conjunto de pontos
xi = x0 ± ih, i = 1, 2, . . . , N.
Nos pontos dessa malha sera˜o calculadas aproximac¸o˜es para uma func¸a˜o y(x) e suas
derivadas. A ferramenta matema´tica ba´sica no ca´lculo de aproximac¸o˜es para as derivadas e´
a se´rie de Taylor que relaciona valores da func¸a˜o e suas derivadas, num ponto x, com valores
dessa mesma func¸a˜o numa vizinhac¸a de x, ou seja, com y(x+ h). Se y(x) tem derivadas ate´
a ordem n+ 1 em x, podemos escrever:
y(x+ h) = y(x) + hy′(x) +
h2
2!
y′′(x) + · · ·+ h
n
n!
y(n)(x) +
+
hn+1
(n+ 1)!
y(n+1)(ξ), x < ξ < x+ h. (1.1)
O u´ltimo termo da expressa˜o (1.1) representa o erro da aproximac¸a˜o de y(x + h) pelo
polinoˆmio (na varia´vel h) de grau n:
Pn(h) = y(x) + y
′(x)h+
y′′(x)
2!
h2 + . . .+
y(n)(x)
n!
hn.
Se n = 1 em (1.1) obtemos uma aproximac¸a˜o para a derivada y′(x), conhecida como fo´rmula
progressiva, que e´ dada por:
y′(x) =
y(x+ h)− y(x)
h
− h
2
y′′(ξ).
Utilizaremos a notac¸a˜o:
∆y(x) = y(x+ h)− y(x)
para representar a diferenc¸a progressiva de y(x). O termo −hy′′(ξ)/2 representa o erro dessa
aproximac¸a˜o. Portanto teremos:
10 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
y′(x) =
1
h
∆y(x)− h
2
y′′(ξ).
De modo semelhante, tomando −h em (1.1), ainda com n = 1, obtemos a fo´rmula
regressiva que utiliza a diferenc¸a regressiva (∇y(x)) e seu erro, ou seja:
y′(x) =
y(x)− y(x− h)
h
+
h
2
y′′(ξ)
=
1
h
∇y(x) + h
2
y′′(ξ).
Tomando n = 2 em (1.1), e reescrevendo (1.1) para h e −h, respectivamente, obtemos
y(x+ h) = y(x) + hy′(x) +
h2
2!
y′′(x) +
h3
3!
y′′′(ξ1)
e
y(x− h) = y(x)− hy′(x) + h
2
2!
y′′(x) − h
3
3!
y′′′(ξ2).
Subtraindo a u´ltima expressa˜o da penu´ltima obtemos a fo´rmula centrada ou fo´rmula de
diferenc¸as centrais
y′(x) =
y(x+ h)− y(x− h)
2h
− h
2
3!
y′′′(ξ) (1.2)
onde ξ ∈ (x− h, x+ h) e foi utilizado o teorema do valor intermedia´rio va´lido para func¸o˜es
cont´ınuas:
1
2
(y′′′(ξ1) + y′′′(ξ2)) = y′′′(ξ), para algum ξ ∈ [min{ξ1, ξ2},min{ξ1, ξ2}].
A notac¸a˜o utilizada para denotar a diferenc¸a central e´ δhy(x) e a equac¸a˜o (1.2) nos fornece
tambe´m o erro dessa fo´rmula, ou seja:
y′(x) =
1
2h
δhy(x)− h
2
3!
y′′′(ξ).
Observac¸a˜o 1.2.1 Quando se fizer necessa´rio deixar claro a importaˆncia do passo h para
um operador que esta´ sendo utilizado, escreveremos: δh, δ
2
h,∆h,∇h, etc..
Erro e Ordem de Aproximac¸a˜o de Uma Fo´rmula de Diferenc¸a
Seja F(x, h) uma fo´rmula de diferenc¸a para aproximac¸a˜o da derivada de ordem q de uma
func¸a˜o y(x) com erro E(x, h). Enta˜o:
y(q)(x) = F(x, h) + E(x, h).
1.2. APROXIMAC¸O˜ES DE DERIVADAS POR DIFERENC¸AS FINITAS 11
Dizemos que a fo´rmula F(x, h) e´ de ordem p se E(x, h) = hpR(x), ondeR(x) na˜o depende de
h. Nesse caso usamos a notac¸a˜o E(x, h) = O(hp). Essa notac¸a˜o significa que limh→0 E(x,h)hp
e´ uma constante finita. Por exemplo no caso da fo´rmula centrada temos que:
F(x, h) = y(x+ h)− y(x− h)
2h
e E(x, h) = −h
2
3!
y′′′(ξ).
de forma que essa fo´rmula e´ de segunda ordem.
Seguindo as mesmas ide´ias podemos estabelecer uma expressa˜o para o ca´lculo aproxima-
do da segunda derivada.Tomando n = 3 em (1.1) com h e −h obtemos:
y(x+ h) = y(x) + hy′(x) +
h2
2!
y′′(x) +
h3
3!
y′′′(x) +
h4
4!
y(4)(ξ1)
e
y(x− h) = y(x)− hy′(x) + h
2
2!
y′′(x)− h
3
3!
y′′′(x) +
h4
4!
y(4)(ξ2).
Somando estas duas u´ltimas expresso˜es e isolando y′′(x) obtemos:
y′′(x) =
y(x+ h)− 2y(x) + y(x− h)
h2
− h
2
12
y(4)(ξ),
=
1
h2
δ2h
2
y(x)− h
2
12
y(4)(ξ).
onde ξ ∈ (x− h, x+ h). Neste caso,
F(x, h) = y(x+ h)− 2y(x) + y(x− h)
h2
, E(x, h) = −h
2
12
y(4)(ξ)
e a fo´rmula e´ de ordem 2.
Observac¸a˜o 1.2.2 O operador δ2h
2
deve ser interpretado como a composic¸a˜o do operador
δh
2
com ele mesmo. O leitor deve observar que nessa composic¸a˜o utilizamos o operador δ
com passo h2 . Isto e´ feito para evitar a utilizac¸a˜o de pontos muito distantes daquele em que a
derivada esta´ sendo aproximada. Observe tambe´m que δ2y(x) = ∆∇y(x), e que utilizaremos
as notac¸o˜es simplificadas δ2 = δ2h
2
e δ = δh
2
.
Outras expresso˜es podem ser obtidas considerando-se mais pontos. Por exemplo, a com-
binac¸a˜o linear
y′(x) = c1y(x+ h) + c2y(x) + c3y(x− h)
impo˜e, por expansa˜o em se´rie de Taylor em torno de x e por comparac¸a˜o dos termos de
mesma ordem em h, as seguintes condic¸o˜es sobre os coeficientes c1, c2 e c3:
c1 + c2 + c3 = 0
c1 − c3 = 1h
c1 + c3 = erro,
12 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
onde erro = coeficiente multiplicando a segunda derivada.
A soluc¸a˜o e´ simples quando c2 =
α
h , ou seja,
c1 = (1− α) 1
2h
, c3 = −(1 + α) 1
2h
, erro = −α
h
.
Assim
y′(x) =
1
2h
{(1− α)y(x + h) + 2αy(x) − (1 + α)y(x − h)} − αh
2
y′′(x)
− h
2
3!
y′′′(x) + 0(h3).
Note que apenas α = 0 fornece erro = O(h2), e nesse caso a fo´rmula resultante e´
a fo´rmula centrada. Os casos α = −1 e α = 1 , fornecem, respectivamente, as fo´rmulas
progressiva e regressiva. Essa te´cnica permite a construc¸a˜o das mais diversas fo´rmulas. Por
exemplo:
2hy′(x) = −3y(x) + 4y(x+ h)− y(x+ 2h) +O(h2).
2hy′(x) = 3y(x)− 4y(x− h) + y(x− 2h) +O(h2).
2hy′(x) = −y(x+ 2h) + 8y(x+ h)− 8y(x− h) + y(x− 2h) +O(h4).
O uso de operadores de diferenc¸as finitas e´ muito u´til para se obter aproximac¸o˜es de
ordem mais elevada para as derivadas. Com a intenc¸a˜o de ilustrar a importaˆncia dos op-
eradores na obtenc¸a˜o de fo´rmulas de diferenc¸as finitas, apresentamos a seguir a derivac¸a˜o
de algumas dessas fo´rmulas. Informac¸o˜es mais detalhadas podem ser obtidas, por exemplo,
em [Ames 1992] e [31].
Alguns desses operadores mais utilizados no contexto de diferenc¸as finitas sa˜o listados
abaixo:
Deslocamento Ey(x) = y(x+ h)
Derivada Dy(x) = y′(x)
Diferenc¸a Centrada δy(x) = y(x+ h2 )− y(x− h2 )
Operador Me´dia µy(x) = 12 (y(x +
h
2 ) + y(x− h2 ))
Diferenc¸a Progressiva ∆y(x) = y(x+ h)− y(x)
= (E − 1)y(x) =⇒ E = 1 +∆.
(1.3)
Pela fo´rmula de Taylor tem-se
Ey(x) = (1 + hD +
h2
2
D2 + · · ·)y(x) = ehDy(x) (1.4)
logo, formalmente podemos escrever:
hD = lnE = ln(1 + ∆). (1.5)
1.2. APROXIMAC¸O˜ES DE DERIVADAS POR DIFERENC¸AS FINITAS 13
Utilizando agora a expansa˜o em se´rie de Taylor da func¸a˜o ln(1+x) = x− 12x2+ 13x3−· · ·,
va´lida para |x| < 1, podemos definir o operador ln(1 + ∆) como:
ln(1 + ∆)y(x) = (∆− 1
2
∆2 +
1
3
∆3 − · · ·)y(x).
Assim a equac¸a˜o (1.5) pode ser reescrita na forma:
hy′(x) = (∆− 1
2
∆2 +
1
3
∆3 − · · ·)y(x)
constituindo-se numa fo´rmula que expressa a primeira derivada de uma func¸a˜o em termos
de suas diferenc¸as progressivas. Portanto, aproximac¸o˜es de qualquer ordem podem ser obti-
das considerando-se na se´rie acima o nu´mero adequado de termos, pagando-se o prec¸o do
aumento do nu´mero de pontos envolvidos na fo´rmula.Por exemplo uma aproximac¸a˜o de
ordem 2, que utiliza treˆs pontos, e´ obtida considerando-se apenas os 2 primeiros termos da
soma. Observe o desenvolvimento abaixo:
hy′(x) =
(
∆− 1
2
∆2
)
y(x)
= y(x+ h)− y(x)− 1
2
(y(x+ 2h)− y(x+ h)− y(x+ h) + y(x))
= −1
2
y(x+ 2h) + 2y(x+ h)− 3
2
y(x).
que utiliza treˆs pontos, para aproximar a derivada.
Outras fo´rmulas que podem ser obtidas pela manipulac¸a˜o de operadores sa˜o:
hy′(x+ h) = (∆ +
1
2
∆2 − 1
6
∆3 + · · ·)y(x) (1.6)
h2y′′(x) = (∆2 −∆3 + 11
12
∆4 − 5
6
∆5 + · · ·)y(x) (1.7)
h2y′′(x+ h) = (∆2 − 1
12
∆4 +
1
12
∆5 + · · ·)y(x) (1.8)
hy′(x) = µ(δ − 1
6
δ3 +
1
30
δ5 − · · ·)y(x) (1.9)
e
h2y′′(x) = (δ2 − 1
12
δ4 +
1
90
δ6 − 1
560
δ8 + · · ·)y(x). (1.10)
Observac¸a˜o 1.2.3 E´ importante notar que, na pra´tica, os valores da func¸a˜o efetivamente
utilizados numa fo´rmula, para cada ponto da malha, na˜o sera˜o os valores exatos e sim val-
ores calculados. Dessa forma ale´m do erro de truncamento dessa fo´rmula estara˜o tambe´m
presentes, os erros de arredondamento do computador. Consequentemente, quando aprox-
imamos uma equac¸a˜o diferencial, as derivadas sa˜o substitu´ıdas por aproximac¸o˜es como as
14 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
que podem ser deduzidas a partir de (1.6), (1.7), (1.8), (1.9) e (1.10) mas com os valores
efetivamente calculados, ou seja,
hy′(xi) ≃ yi+1 − yi
hy′(xi) ≃ yi − yi−1
2hy′(xi) ≃ yi+1 − yi−1
e
h2y′′(xi) ≃ yi+1 − 2yi + yi−1,
onde yi e´ uma aproximac¸a˜o de y(xi).
Generalizac¸a˜o Para Duas Varia´veis
A generalizac¸a˜o para mais de uma varia´vel e´ bastante simples. Como ilustrac¸a˜o vamos
tratar o caso de duas varia´veis x e t. Considere uma malha no plano x − t como sendo o
conjunto de pontos (xi, tj) = (x0+ ih, t0+ jk), portanto com espac¸amento h em x e k em t.
As diferenc¸as finitas ja´ obtidas em uma dimensa˜o podem agora ser utilizadas em cada
varia´vel para gerar aproximac¸o˜es para as derivadas parciais de uma func¸a˜o de va´rias varia´veis.
Assim, temos as seguintes fo´rmulas: progressiva
ut(x, t) =
u(x, t+ k)− u(x, t)
k
− k
2
utt(x, ζ)
=
1
k
∆tu(x, t)− k
2
utt(x, ζ) (t < ζ < t+ k) (1.11)
regressiva
ut(x, t) =
u(x, t)− u(x, t− k)
k
+
k
2
utt(x, ζ)
=
1
k
∇tu(x, t) + k
2
utt(x, ζ), (t− k < ζ < t) (1.12)
central
ux(x, t) =
u(x+ h, t)− u(x− h, t)
2h
− h
2
6
uxxx(ξ, t)
=
1
2h
δxu(x, t)− h
2
6
uxxx(ξ, t), (x − h < ξ < x+ h)
uxx(x, t) =
u(x+ h, t)− 2u(x, t) + u(x− h, t)
h2
− h
2
12
uxxxx(ξ, t)
=
1
h2
δ2xu(x, t)−
h2
12
uxxxx(ξ, t), (x− h < ξ < x+ h) (1.13)
1.3. PROBLEMA DE VALOR INICIAL EM EQUAC¸O˜ES ORDINA´RIAS 15
utt(x, t) =
u(x, t+ k)− 2u(x, t) + u(x, t− k)
k2
− k
2
12
utttt(x, ζ)
=
1
k2
δ2t u(x, t)−
k2
12
utttt(x, ζ), (t− h < ζ < t+ h)
uxt(x, t) =
u(x+ h, t+ k)− u(x+ h, t− k)− u(x− h, t+ k) + u(x− h, t− k)
4hk
− h
2
6
uxxxt(ξ1, ζ1)− k
2
6
uxttt(ξ2,ζ2) x− h < ξ1,ξ2 < x+ h, e
t− k < ζ1,ζ2 < t+ k,
onde denotamos por ut(x, t) a derivada parcial da func¸a˜o u com relac¸a˜o a t e por ∆tu(x, t)
a diferenc¸a progressiva de u na varia´vel t.
1.3 Problema de Valor Inicial em Equac¸o˜es Ordina´rias
Uma equac¸a˜o diferencial na qual a varia´vel dependente e´ func¸a˜o de apenas uma varia´vel, e´
dita equac¸a˜o diferencial ordina´ria. Se mais de uma varia´vel independente esta´ presente enta˜o
teremos uma equac¸a˜o diferencial parcial, como ja´ haviamos dito oportunamente. A teoria
das equac¸o˜es parciais e´ consideravelmente mais complexa do que a das equac¸o˜es ordina´rias.
Essa observac¸a˜o e´ um tanto mais verdadeira no contexto do tratamento nume´rico, uma vez
que a presenc¸a de um nu´mero maior de varia´veis introduz maior chance de erro. O tema
principal deste texto e´ o estudo das equac¸o˜es parciais, no entanto, por razo˜es dida´ticas, os
conceitos fundamentais sera˜o apresentados no contexto das equac¸o˜es diferenciais ordina´rias.
Um problema de valor inicial (PVI) e´ um problema de evoluc¸a˜o, no qual a informac¸a˜o
inicial (conhecida), e´ propagada para o interior do domı´nio pela equac¸a˜o diferencial.
Matematicamente, o mais simples dos problemas de valor inicial pode ser apresentado
na forma:
y′ = f(x, y) (1.14)
y(a) = α,
onde f : IR2 → IR e´ uma func¸a˜o cont´ınua. A func¸a˜o y = y(x) (x ≥ a) e´ a soluc¸a˜o e α e´ o
valor inicial em a. Equac¸o˜es ordina´rias de ordem mais alta podem sempre ser transformadas
nessa forma, desde que estejamos preparados para admitir que y seja um vetor. Por exemplo,
a equac¸a˜o de segunda ordem,
y′′ = f(x, y, y′) (1.15)
y(a) = α, y′(a) = β,
com soluc¸a˜o y(x) e os valores iniciais α e β dados no ponto a, pode ser reescrita na forma
vetorial fazendo:
y′ = v, y(a) = α
v′ = f(x, y, v), v(a) = β.
16 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
Considerando Y =
(
y
v
)
e γ =
(
α
β
)
o problema acima transforma-se na equac¸a˜o
vetorial de primeira ordem
Y ′ = F (x, Y ),
Y (a) = γ
onde
F (x, Y ) =
(
v
f(x, y, v)
)
.
Um ponto fundamental na resoluc¸a˜o de uma equac¸a˜o diferencial, diz respeito a` existeˆncia
e unicidade de soluc¸o˜es. Por princ´ıpio, so´ tem sentido a soluc¸a˜o nume´rica de um problema
quando este tem soluc¸a˜o. Pode-se ainda argumentar que para determinarmos a soluc¸a˜o
nume´rica de um problema esta deveria ser u´nica, pois caso contra´rio o processo nume´rico
poderia ficar indeciso sobre qual soluc¸a˜o perseguir. No caso de um problema de valor inicial
de equac¸o˜es ordina´rias esta dificuldade esta´ bem resolvida.
Vamos supor que a func¸a˜o f(x, y) satisfaz as seguintes condic¸o˜es:
1. f : E → IRn e´ cont´ınua, onde E = {(x, y), x ∈ [a, b], y ∈ IRn};
2. existe uma constante L tal que para todo x ∈ [a, b] e quaisquer dois vetores y e y∗ em
IRn
|f(x, y)− f(x, y∗)| ≤ L|y − y∗|. (1.16)
Neste caso dizemos que a func¸a˜o f satisfaz a condic¸a˜o de Lipschitz na varia´vel y.
Teorema 1.1 (Existeˆncia e Unicidade-Teorema de Picard) Seja f(x, y) satisfazendo
as condic¸o˜es acima e seja α um vetor dado. Enta˜o, existe exatamente uma func¸a˜o y(x) com
as seguintes propriedades:
i) y = y(x) e´ cont´ınua e diferencia´vel para x em [a, b];
ii) y′(x) = f(x, y(x)), x ∈ [a, b];
iii) y(a) = α.
Para a prova ver [Coddington 1955].
Exemplos de Me´todos Nume´ricos
Euler em 1768, usando diferenc¸as progressivas, desenvolveu uma aproximac¸a˜o para (1.14),
que posteriormente descobriu-se tratar-se de um me´todo nume´rico pertencente a` classe dos
me´todos de passo mu´ltiplo. Esse me´todo e´ hoje conhecido como Me´todo de Euler. Sua
deduc¸a˜o segue os passos abaixo, acompanhe: O intervalo [a, b] e´ dividido em N partes
iguais, cada uma de comprimento h, formando um conjunto discreto, com x0 = a e xN = b,
Rh = {xi = a+ ih, i = 0, 1, 2, . . . , N}. Aqui h = (b − a)/N .
Sejam yi ≈ y(xi) aproximac¸o˜es para y(xi), i = 1, 2, . . .N e y0 = α. Em cada um dos
pontos x0, x1, . . . , xN−1, aproximamos a equac¸a˜o diferencial (1.14) por
1.3. PROBLEMA DE VALOR INICIAL EM EQUAC¸O˜ES ORDINA´RIAS 17
∆yi =
yi+1 − yi
h
= f(xi, yi), i = 0, 1, . . . , N − 1
ou,
yi+1 − yi = hf(xi, yi), i = 0, 1, . . . , N − 1.
Este e´ o me´todo de Euler que e´ um me´todo expl´ıcito. Note que o valor de yi+1 e´ calculado
explicitamente em func¸a˜o de yi, mesmo quando a func¸a˜o f e´ na˜o linear. Se no lugar de ∆
fosse usado ∇ ter´ıamos a versa˜o impl´ıcita do Me´todo de Euler:
∇yi = yi − yi−1
h
= f(xi, yi), i = 1, . . . , N
ou,
yi+1 − yi = hf(xi+1, yi+1), i = 0, 1, . . . , N − 1. (1.17)
Note que no caso em que f e´ na˜o linear, a equac¸a˜o (1.17) define o vetor yi+1 implicitamente
e portanto seu ca´lculo exige a utilizac¸a˜o de me´todos para a soluc¸a˜o de sistemas de equac¸o˜es
alge´bricas na˜o lineares, como por exemplo, o me´todo das aproximac¸o˜essucessivas ou o
me´todo de Newton.
Observe tambe´m que a equac¸a˜o diferencial no processo de discretizac¸a˜o, foi substitu´ıda
por uma outra equac¸a˜o chamada equac¸a˜o de diferenc¸as que tem sua pro´pria soluc¸a˜o inde-
pendente daquela da equac¸a˜o diferencial. Assim o me´todo nume´rico sera´ considerado eficaz
se a soluc¸a˜o da equac¸a˜o de diferenc¸as tiver comportamento similar ao da soluc¸a˜o da equac¸a˜o
diferencial. Para mais informac¸o˜es sobre equac¸o˜es de diferenc¸as, ver [18].
Exemplo 1.3.1 Considere o PVI,
y′ = −y + x
y(0) = 1 (1.18)
cuja soluc¸a˜o exata e´ y(x) = x− 1 + 2e−x, de onde obtemos:
y(0.1) = 0.9097 y(0.3) = 0.7816 y(0.5) = 0.7131 y(0.7) = 0.6932 y(0.9) = 0.7131
y(0.2) = 0.8375 y(0.4) = 0.7406 y(0.6) = 0.6976 y(0.8) = 0.6986 y(1.0) = 0.7357.
A func¸a˜o f e´ dada por: f(x, y) = −y+ x e, portanto o me´todo de Euler expl´ıcito toma
a forma:
yi+1 = yi + h(−yi + xi).
Tomando h = 0.2, encontraremos uma soluc¸a˜o aproximada em [0, 1] resolvendo a equac¸a˜o
de diferenc¸as:
yi+1 = yi − 0.2yi + 0.2xi ou seja yi+1 = 0.8yi + 0.2xi,
cuja soluc¸a˜o para n = 1, 2, . . . , 5 e´:
y0 = 1.0 = y(0.0) (condic¸a˜o inicial)
y1 = 0.8y0 + 0.2x0 = 0.8000 ≃ y(0.2)
y2 = 0.8y1 + 0.2x1 = 0.6800 ≃ y(0.4)
y3 = 0.8y2 + 0.2x2 = 0.6240 ≃ y(0.6)
y4 = 0.8y3 + 0.2x3 = 0.6192 ≃ y(0.8)
y5 = 0.8y4 + 0.2x4 = 0.6554 ≃ y(1.0).
18 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
Tomando h = 0.1, obtemos um problema semelhante, mas com 10 valores de y para serem
calculados. Um procedimento ana´logo produz as aproximac¸o˜es:
y1 = 0.9000 y3 = 0.7580 y5 = 0.6810 y7 = 0.6570 y9 = 0.6748
y2 = 0.8200 y4 = 0.7122 y6 = 0.6629 y8 = 0.6609 y10 = 0.6974.
Note que quando h = 0.1 calculamos a aproximac¸a˜o em um nu´mero maior de pontos e
portanto devemos tomar o cuidado de comparar as aproximac¸o˜es para h = 0.2 e h = 0.1 nos
mesmos pontos, por exemplo a aproximac¸a˜o y1 obtida com h = 0.2 deve ser comparada com
a aproximac¸a˜o y2 do passo h = 0.1 e assim por diante. Vemos assim que a soluc¸a˜o obtida
com h = 0.1 e´ melhor que aquela para h = 0.2, isto ocorre porque a equac¸a˜o de diferenc¸as
para h = 0.1 representa a equac¸a˜o diferencial de maneira mais fiel do que para h = 0.2.
O me´todo de Euler impl´ıcito para o problema (1.18) e´ dado pela equac¸a˜o de diferenc¸as:
yi+1 = yi + h(−yi+1 + xi+1),
donde yi+1 =
1
1+h (yi + hxi+1) fornece os seguintes valores para h = 0.1,
y1 = 0.9182 y3 = 0.8026 y5 = 0.7418 y7 = 0.7263 y9 = 0.7482
y2 = 0.8529 y4 = 0.7660 y6 = 0.7289 y8 = 0.7330 y10 = 0.7711.
A figura 1.1 ilustra o comportamento dos me´todos de Euler expl´ıcito e impl´ıcito compara-
dos com a soluc¸a˜o exata. Observamos que o me´todo expl´ıcito erra por falta, ja´ o impl´ıcito o
faz por excesso.
E´ natural, portanto, considerarmos a combinac¸a˜o linear desses dois me´todos, com o ob-
jetivo de que os erros se cancelem e possamos obter uma melhor aproximac¸a˜o. Por exemplo,
uma combinac¸a˜o linear e´ a me´dia aritme´tica, que fornece o Me´todo dos Trape´zios, tambe´m
conhecido como Regra dos Trape´zios:
yi+1 − yi = h
2
(f(xi, yi) + f(xi+1, yi+1)), i = 0, 1, . . . , N − 1.
tambe´m impl´ıcito e que, para o PVI em questa˜o, produz:
y1 = 0.9095 y3 = 0.7816 y5 = 0.7130 y7 = 0.6926 y9 = 0.7125
y2 = 0.8372 y4 = 0.7402 y6 = 0.6970 y8 = 0.6981 y10 = 0.7351.
O me´todo dos trape´zios e´ claramente mais preciso que ambos os me´todos anteriores,
conforme pode ser observado na figura 1.2.
1.3.1 Convergeˆncia
Seja y(xn) a soluc¸a˜o exata do PVI (1.14) no ponto xn e yn uma aproximac¸a˜o obtida
por um me´todo nume´rico para essa soluc¸a˜o.
Definic¸a˜o 1.1 O erro global no ponto xn e´ definido por:
en = y(xn)− yn.
1.3. PROBLEMA DE VALOR INICIAL EM EQUAC¸O˜ES ORDINA´RIAS 19
0 0.2 0.4 0.6 0.8 1
0.4
0.6
0.8
1
1.2
* Euler Implicito
x Solucao Exata
+ Euler Explicito
 x
 
y
Figura 1.1: Soluc¸a˜o do problema (1.18) pelos me´todos de Euler Expl´ıcito e
Impl´ıcito com h = 0.1
A ana´lise do erro global nos permite estabelecer a convergeˆncia de um me´todo. Esta e´
uma propriedade muito deseja´vel para um me´todo nume´rico, pois ela ira´ nos garantir que
ao refinarmos a malha o resultado nume´rico se torna mais pro´ximo do resultado exato. A
definic¸a˜o formal de convergeˆncia e´ a seguinte:
Definic¸a˜o 1.2 Dizemos que um me´todo nume´rico e´ convergente se o erro global en converge
para zero quando n tende para infinito de maneira que o ponto xn permanec¸a fixo.
Isto quer dizer que a convergeˆncia esta´ sendo analisada no ponto xn = a+nh e para que
este ponto permanec¸a fixo, a quantidade nh deve permanecer fixa, portanto se n tende para
infinito necessariamente h deve tender a zero, ou seja, a malha esta´ sendo refinada para cada
novo valor de n da sequeˆncia. E´ poss´ıvel verificar a convergeˆncia de um me´todo nume´rico
diretamente a partir de sua definic¸a˜o, como faremos a seguir para o me´todo dos trape´zios.
No entanto esta e´ uma maneira muito trabalhosa e e´ enta˜o prefer´ıvel fazeˆ-lo indiretamente
utilizando os conceitos de consisteˆncia e zero-estabilidade.
Ana´lise da convergeˆncia do me´todo dos trape´zios
A equac¸a˜o de diferenc¸as para o me´todo dos trape´zios e´:
yn+1 = yn +
h
2
[f(xn, yn) + f(xn+1, yn+1)]
20 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
0 0.2 0.4 0.6 0.8 1
0.4
0.6
0.8
1
1.2
x Solucao Exata
+ Regra dos Trapezios
x
y
Figura 1.2: Soluc¸a˜o do problema (1.18) pelo me´todo dos trape´zios com h = 0.1
e a equac¸a˜o de diferenc¸as para a soluc¸a˜o exata, utilizando a se´rie de Taylor e´:
y(xn+1) = y(xn) +
h
2
[f(xn, y(xn)) + f(xn+1, y(xn+1))]− h
3
12
y′′′(ξ).
Subtraindo a primeira equac¸a˜o da segunda obtemos a seguinte relac¸a˜o para o erro global
en:
en+1 = en +
h
2
[f(xn, y(xn))− f(xn, yn) + f(xn+1, y(xn+1))− f(xn+1, yn+1)]− h
3
12
y′′′(ξ).
Tomando o valor absoluto e utilizando a condic¸a˜o de Lipschitz de f , equac¸a˜o (1.16), obtemos
a desigualdade:
|en+1| ≤ |en|+ h
2
(L|en|+ L|en+1|) + h
3
12
|y′′′(ξ)|.
Supondo agora que |y′′′(x)| e´ limitada por 12M no intervalo de interesse, a desigualdade
acima transforma-se em:
(1− hL
2
)|en+1| ≤ (1 + hL
2
)|en|+ h3M.
Admitindo que h seja suficientemente pequeno de forma que (1−hL2 ) > 0, podemos reescrever
1.3. PROBLEMA DE VALOR INICIAL EM EQUAC¸O˜ES ORDINA´RIAS 21
a expressa˜o acima como:
|en+1| ≤
(
1 + hL2
1− hL2
)
|en|+ h
3M
1− hL2
.
E´ poss´ıvel mostrar por induc¸a˜o sobre n (ver exerc´ıcio (1.8)) que a soluc¸a˜o ξn da equac¸a˜o de
diferenc¸as
ξn+1 = Aξn +B, com A =
1 + hL2
1− hL2
e B =
h3M
1− hL2
(1.19)
domina o erro |en|, isto e´ |en| ≤ ξn.
A soluc¸a˜o geral de (1.19) e´:
ξn =
h2M
L
[(
1 + hL2
1− hL2
)n
− 1
]
. (1.20)
Reescrevendo o termo 1+hL2 na forma 1−hL2 +hL, substituindo no quociente entre pareˆnteses
da expressa˜o (1.20) e simplificando os fatores comuns, esse quociente transforma-se em:
1 + h
(
L
1− hL2
)
.
Utilizando agora o fato que 1 + x ≤ ex para x ≥ 0, podemos escrever:(
1 + hL2
1− hL2
)n
=
(
1 + h
(
L
1− hL2
))n
≤ enh
(
L
1−hL
2
)
e portanto a expressa˜o (1.20) transforma-se na desigualdade:
ξn ≤ h
2M
L
(
e
nh
(
L
1−hL
2
)
− 1
)
.
Lembrando que xn = a+ nh e portanto nh = xn − a a expressa˜o acima torna-se:
ξn ≤ h
2M
L
(
e
(xn−a)
(
L
1−hL
2
)
− 1
)
.
Como |en| ≤ ξn obtemos o resultado final na forma:
|en| ≤ h
2M
L
(
e
(xn−a)
(
L
1−hL
2
)
− 1
)
. (1.21)
22 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
Tomando n tendendo para infinito, ou equivalentemente h tendendo para zero, vemos que
|en| tende para zero com ordem O(h2), pois o termo entre pareˆnteses tende para eL(x−a)−1
quando h → 0, onde x (fixo) e´ o ponto onde estamos analisando a convergeˆncia. Agora, o
termo eL(x−a) − 1 e´ limitadoindependentemente de h, o que demonstra a convergeˆncia da
regra dos trape´zios.
Exemplo 1.3.2 Apresentamos a seguir um exemplo ilustrando o fato de que o limitante
para o erro global obtido em (1.21) e´ em geral bastante conservador. Consideremos o PVI:
y′ = −2y, y(0) = 1 , 0 ≤ x ≤ 1.
cuja soluc¸a˜o e´: y(x) = e−2x. Resolvendo esse problema pela regra dos trape´zios com
h = 0.1 e com h = 0.01 obtemos os seguintes valores para o erro global ma´ximo: emax =
0.00123160912 no primeiro caso e emax = 0.0000122631795 no segundo. Por outro lado, cal-
culando o limitante do erro global atrave´s da fo´rmula (1.21), obtemos os valores 0.005727606
e 0.00005727606 respectivamente para h = 0.1 e h = 0.01. Comparando esses valores, vemos
que o limitante e´ 5 vezes o erro verdadeiro.
Note que o conceito de convergeˆncia aqui definido na˜o e´ de aplicac¸a˜o pra´tica direta, pois
para calcular o erro global precisamos da soluc¸a˜o exata que na˜o esta´ dispon´ıvel. Portanto
para testar na pra´tica a convergeˆncia de um me´todo nume´rico, resolve-se o problema para
va´rios valores decrescentes de h e observa-se se a sequeˆncia assim obtida esta´ se aproxi-
mando de algum valor fixado, ou seja testa-se, de maneira emp´ırica, se essa sequeˆncia e´ de
Cauchy e portanto convergente. No entanto, o fato da sequeˆncia ser convergente na˜o implica
que ela converge para um limite que e´ a soluc¸a˜o da equac¸a˜o diferencial. Assim, algumas
questo˜es pra´ticas relacionadas ao conceito de convergeˆncia merecem ser salientadas, como
por exemplo:
• Havendo convergeˆncia no sentido pra´tico, a soluc¸a˜o nume´rica obtida representa a
soluc¸a˜o do problema?
• A discretizac¸a˜o nume´rica introduz alguma caracter´ıstica (“ru´ıdo”) que na˜o esta´ pre-
sente na equac¸a˜o original? Neste caso e´ este ru´ıdo controla´vel ou na˜o?
• Qua˜o ra´pida e´ a convergeˆncia?
As respostas a essas perguntas esta˜o associadas aos conceitos de Consisteˆncia, Zero
Estabilidade e Ordem de Convergeˆncia que estudaremos logo a seguir.
Consisteˆncia
Parece natural admitir-se que a soluc¸a˜o nume´rica aproxime a soluc¸a˜o teo´rica do problema
em estudo. No entanto, este na˜o e´ sempre o caso pois a equac¸a˜o de diferenc¸as tem vida
pro´pria e portanto sua pro´pria soluc¸a˜o, que na˜o necessariamente guarda relac¸a˜o alguma com
a soluc¸a˜o do problema original, a menos que assim seja explicitamente exigido desta equac¸a˜o
de diferenc¸as. Essa propriedade, chamada consisteˆncia, e´ imposta a` equac¸a˜o de diferenc¸as e
“amarra-a” a` equac¸a˜o diferencial. Inversamente, a soluc¸a˜o do problema cont´ınuo na˜o e´, em
geral, soluc¸a˜o da equac¸a˜o de diferenc¸as e o erro cometido ao substituirmos a soluc¸a˜o exata
na equac¸a˜o de diferenc¸as e´ chamado de erro de truncamento local. Se o erro de truncamento
1.3. PROBLEMA DE VALOR INICIAL EM EQUAC¸O˜ES ORDINA´RIAS 23
local e´ de O(hp), com p ≥ 1 dizemos que esse me´todo tem ordem de consisteˆncia p, ou
simplesmente que esse me´todo e´ de ordem p.
O erro de truncamento local pode tambe´m ser interpretado como o erro cometido ao
calcularmos a soluc¸a˜o num ponto xi+1 sem levar em considerac¸a˜o os erros dos passos ante-
riores, isto e´, levando em conta somente o erro cometido no passo atual. Para o me´todo de
Euler, suponha que ate´ o ponto xi todos os ca´lculos foram executados exatamente e que o
valor em xi+1 devera´ ser aproximado atrave´s do me´todo nume´rico. Essa aproximac¸a˜o tera´
enta˜o erro de truncamento local τi+1, dado por:
y(xi+1)− y(xi)
h
= f(xi, y(xi)) + τi+1
m
y′(xi)− h
2
y′′(ξ) = f(xi, y(xi)) + τi+1
Podemos interpretar ainda o erro de truncamento local como aquele erro introduzido ao
substituirmos a equac¸a˜o diferencial por uma equac¸a˜o de diferenc¸as. Portanto, para o me´todo
de Euler com h → 0 a equac¸a˜o de diferenc¸as aproxima-se da equac¸a˜o diferencial com erro
de truncamento local τi+1 = −h2 y′′(ξ), ou seja, com ordem O(h). Dizemos que o me´todo
de Euler e´ consistente de ordem 1. Analogamente, o me´todo dos trape´zios e´ consistente de
ordem 2, porque o seu erro e´ O(h2).
O erro global e´ formado pela acumulac¸a˜o dos erros de truncamento local a cada passo
juntamente com os erros de arredondamento. Portanto a ordem de consisteˆncia de um
me´todo esta´ relacionada com sua ordem de convergeˆncia. Essa ordem e´ um indicativo da
velocidade com que se da´ a convergeˆncia. Assim um me´todo com erro O(h2) deve convergir
mais rapidamente que outro com erro O(h). Isso significa, pelo menos teoricamente, que
para o mesmo tamanho de passo h um me´todo com ordem mais alta produz aproximac¸o˜es
mais precisas.
Zero Estabilidade
E´ poss´ıvel melhorar a ordem de consisteˆncia de um me´todo simplesmente aumentando-se
a quantidade de informac¸o˜es nas quais esse me´todo esta´ baseado. Por exemplo a fo´rmula
do ponto me´dio que utiliza 3 pontos,
yi+1 − yi−1 = 2hf(xi, yi), i = 1, 2, . . . , N − 1
ou equivalentemente,
yi+2 − yi = 2hf(xi+1, yi+1), i = 0, 1, . . . , N − 2
tem ordem de consisteˆncia 2, e a fo´rmula impl´ıcita
3yi+2 − 4yi+1 + yi = 2hf(xi+2, yi+2), i = 0, 1, . . . , N − 2. (1.22)
tambe´m tem ordem de consisteˆncia 2.
No entanto, como veremos, ao aumentarmos a quantidade de informac¸a˜o utilizada au-
mentamos tambe´m a possibilidade de algo dar errado. Na argumentac¸a˜o abaixo, tentare-
mos convencer o leitor da necessidade de uma ana´lise mais detalhada das propriedades das
24 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
equac¸o˜es de diferenc¸as com o objetivo de esclarecer onde e como esse “algo” pode dar er-
rado. Acompanhe a ana´lise dos exemplos acima, onde estamos aproximando uma equac¸a˜o
diferencial de primeira ordem por uma equac¸a˜o de diferenc¸as de segunda ordem, ou seja,
uma equac¸a˜o que relaciona 3 valores consecutivos. Similarmente a`s equac¸o˜es diferenciais,
uma soluc¸a˜o da equac¸a˜o de diferenc¸as de segunda ordem e´ uma combinac¸a˜o linear de duas
soluc¸o˜es ba´sicas. De forma que, enquanto a equac¸a˜o diferencial tem apenas uma soluc¸a˜o
ba´sica (no caso em que y ∈ IR) a equac¸a˜o de diferenc¸as tem duas distintas. Assim, essas duas
soluc¸o˜es na˜o podem convergir para a soluc¸a˜o do problema, ao mesmo tempo. A soluc¸a˜o que
na˜o converge para a soluc¸a˜o exata e´ denominada soluc¸a˜o espu´ria a` qual na˜o e´ permitido
crescer de maneira ilimitada, sob pena de inviabiliziar a aplicac¸a˜o do me´todo. O efeito das
soluc¸o˜es espu´rias na soluc¸a˜o nume´rica e´, comumente, o aparecimento de oscilac¸o˜es. Quando
o me´todo e´ razoa´vel essas oscilac¸o˜es decrescem com a evoluc¸a˜o, o que significa que as soluc¸o˜es
espu´rias decaem para zero a` medida que i cresce. Um primeiro teste para esta condic¸a˜o de
decaimento e´ feito para o caso em que f(x, y) = 0, impondo que a soluc¸a˜o do me´todo deva
oscilar em torno de uma constante. Se o me´todo e´ razoa´vel, as oscilac¸o˜es se dissipam e o
me´todo e´ dito ser zero-esta´vel; caso contra´rio as oscilac¸o˜es sa˜o amplificadas prejudicando a
convergeˆncia.
Definic¸a˜o 1.3 A condic¸a˜o necessa´ria para um me´todo ser zero-esta´vel e´ que as soluc¸o˜es
ba´sicas da equac¸a˜o de diferenc¸as associada, considerando f(x, y) = 0, sejam limitadas.
As soluc¸o˜es ba´sicas da equac¸a˜o de diferenc¸as sa˜o dadas pelas ra´ızes do polinoˆmio caracte-
r´ıstico a ela associado. Por exemplo, a fo´rmula do ponto me´dio, yi+2−yi = 0 tem polinoˆmio
caracter´ıstico ρ(z) = z2 − 1 = 0, cujas ra´ızes sa˜o 1 e −1, e as soluc¸o˜es ba´sicas, ϕ e ψ, dadas
por ϕn = (1)
n e ψn = (−1)n, n = 1, 2, . . . Ja´ a fo´rmula 3yi+2−4yi+1+yi = 0 tem polinoˆmio
caracter´ıstico ρ(z) = 3z2 − 4z + 1 = 0, cujas ra´ızes sa˜o 1 e 13 , e as soluc¸o˜es ba´sicas, ϕ e ψ,
dadas por ϕn = (1)
n e ψn = (
1
3 )
n, n = 1, 2, . . ..
Note que o comportamento dos dois me´todos e´ distinto; enquanto a fo´rmula do ponto
me´dio conserva a amplitude das oscilac¸o˜es, a fo´rmula (1.22) amortece-as. Como em nenhum
dos dois casos,as oscilac¸o˜es sa˜o amplificadas, isto e´, nenhuma das soluc¸o˜es ba´sicas cresce,
ambos os me´todos sa˜o zero-esta´veis.
As condic¸o˜es: consisteˆncia e zero-estabilidade sa˜o necessa´rias e suficientes para a con-
vergeˆncia de um me´todo [20].
Teorema 1.2 (Equivaleˆncia de Lax) No contexto de problema de valor inicial para
equac¸o˜es ordina´rias, um me´todo baseado em diferenc¸as finitas e´ convergente se e somente
se ele e´ consistente e zero-esta´vel.
E´ importante notar entretanto, que no caso de equac¸o˜es parciais, apesar de existir de-
terminada analogia, os resultados na˜o sa˜o suficientemente abrangentes, ou seja, sa˜o va´lidos
somente em alguns casos particulares. O teorema acima garante a convergeˆncia, mas na˜o
especifica a ordem de convergeˆncia. Como regra, pode-se dizer que a ordem de convergeˆncia
e´ a mesma da consisteˆncia.
Estabilidade
Os conceitos estudados ate´ aqui dizem respeito a h → 0. No entanto, quanto menor o
valor de h, tanto maior sera´ o nu´mero de ca´lculos que devera˜o ser efetuados para produzir
1.3. PROBLEMA DE VALOR INICIAL EM EQUAC¸O˜ES ORDINA´RIAS 25
uma soluc¸a˜o aproximada em determinado ponto x, fixado. Isto pode levar a um acu´mulo,
na˜o controlado, de erros e nesse caso dizemos que esse me´todo e´ na˜o esta´vel ou insta´vel.
A estabilidade (h fixo) e´ muito importante para problemas que exigem um nu´mero muito
grande de aplicac¸o˜es repetidas de um me´todo para obter a soluc¸a˜o em um dado intervalo.
Como em geral na˜o e´ poss´ıvel predizer o nu´mero de passos necessa´rios para obtenc¸a˜o da
soluc¸a˜o, e´ aconselha´vel que, por precauc¸a˜o, tenhamos um controle do comportamento do
erro, isto e´, um controle da propagac¸a˜o do erro. Na verdade, a escolha do tamanho h da
malha esta´ intimamente ligada a` estabilidade do me´todo, e e´ o paraˆmetro atrave´s do qual e´
poss´ıvel exercer esse controle. Podemos interpretar um me´todo nume´rico como uma ma´quina
de produzir nu´meros a partir de dados iniciais. Como esses dados iniciais conte´m sempre
erro (de arredondamento do computador por exemplo), se essa ma´quina tiver o defeito
de amplifica´-los, em pouco tempo o crescimento do erro dominara´ a soluc¸a˜o produzida pela
ma´quina e esta perdera´ seu significado. Vamos ilustrar essas ide´ias pela ana´lise, em paralelo,
de dois me´todos nume´ricos.
Considere os me´todos do ponto me´dio e do trape´zio
yi+2 − yi = 2hf(xi+1, yi+1), i = 0, 1, . . . , N − 2 (1.23)
yi+1 − yi = h
2
(f(xi, yi) + f(xi+1, yi+1)), i = 0, 1, . . . , N − 1. (1.24)
Substituindo a soluc¸a˜o exata nas expresso˜es (1.23-1.24) e considerando a definic¸a˜o de erro
de truncamento local obtemos:
y(xi+2)− y(xi) = 2hf(xi+1, y(xi+1)) + τi+2, i = 0, 1, . . . , N − 2 (1.25)
y(xi+1)− y(xi) = h
2
(f(xi, y(xi)) + f(xi+1, y(xi+1))) + τi+1, i = 0, 1, . . . , N − 1.
(1.26)
Considerando o erro global ei = y(xi)−yi podemos analisar seu comportamento derivando
uma equac¸a˜o que o represente. Na verdade, o erro global esta´ intimamente relacionado com
o erro de truncamento local, e em geral satisfaz uma equac¸a˜o de diferenc¸as obtida a partir
da expressa˜o do me´todo nume´rico. Por exemplo, subtraindo (1.25) de (1.23), tem-se, para
o me´todo do ponto me´dio:
ei+2 − ei = 2h(f(xi+1, y(xi+1))− f(xi+1, yi+1)) + τi+2
= 2hfy(xi+1, ζ)ei+1 + τi+2, i = 0, 1, . . . , N − 2. (1.27)
Aqui foi usado o teorema do valor me´dio e a notac¸a˜o ∂f/∂y = fy.
Analogamente, no caso do me´todo do Trape´zio, subtraindo-se (1.26) de (1.24) obte´m-se:
ei+1 − ei = h
2
(fy(xi, ζ1)ei + fy(xi+1, ζ2)ei+1) + τi+1, i = 0, 1, . . . , N − 1. (1.28)
As equac¸o˜es (1.27)-(1.28) sa˜o aquelas que determinam o erro global como func¸a˜o do erro
de trucamento local para os exemplos da regra do ponto me´dio e trape´zios.
Dessas considerac¸o˜es derivamos a motivac¸a˜o para o uso de uma equac¸a˜o teste para a
verificac¸a˜o da estabilidade. As equac¸o˜es (1.27) e (1.28) sugerem que o erro deve satisfazer
26 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
uma equac¸a˜o do tipo e′ = λe+ g, onde λ = 2hfy(xi+1, ζ) em (1.27) e λ tem uma expressa˜o
um pouco mais complicada no caso (1.28). Obviamente λ na˜o sera´ sempre constante, mas
essas expresso˜es nos motivam a definir a estabilidade, baseada no problema teste y′ = λy,
pois nesse caso, assumindo que o erro de truncamento local satisfaz τi = T , as equac¸o˜es
(1.27) e (1.28) se transformam, respectivamente, em:
ei+2 = 2hλei+1 + ei + T (1.29)
ei+1 = ei +
hλ
2
(ei+1 + ei) + T
que podem ser solucionadas, permitindo uma ana´lise detalhada do seu comportamento.
Essa ana´lise e´ feita para cada me´todo em separado, assim, por exemplo, no caso da regra
do ponto me´dio, como (1.29) e´ uma equac¸a˜o para o erro global, deseja-se que suas soluc¸o˜es
sejam decrescentes para que na˜o haja acu´mulo desse erro. Como ja´ sabemos as soluc¸o˜es
sera˜o combinac¸o˜es lineares de (ξ1)
n e (ξ2)
n onde ξ1 e ξ2 sa˜o as duas ra´ızes do polinoˆmio
π(r, hλ) = r2−2hλr−1 = 0, que e´ chamado de polinoˆmio de estabilidade da regra do ponto
me´dio. Portanto para que o erro seja decrescente precisamos impor a condic¸a˜o de que as
ra´ızes ξ1 e ξ2 tenham mo´dulo menor que 1. Observe que essas ra´ızes sa˜o func¸o˜es de hλ e
portanto a regia˜o do plano complexo h¯ = hλ para a qual essas ra´ızes esta˜o no disco unita´rio
e´ chamada de regia˜o de estabilidade absoluta do me´todo.
O estudo da estabilidade esta´ relacionado com o crescimento ou na˜o do erro global
quando o nu´mero de aplicac¸o˜es do me´todo nume´rico aumenta, mantendo-se o passo h fixo.
Para simplificar a ana´lise, a estabilidade e´ verificada por meio da aplicac¸a˜o do me´todo a`
equac¸a˜o teste y′ = λy. Note que isso significa consideramos T = 0 em (1.29, termo que na˜o
influencia no crescimento de ei. A seguir apresentamos uma breve ana´lise da estabilidade
para os me´todos mais populares.
• Euler expl´ıcito ei+1 − ei = hλei ⇒ ei+1 = (1 + hλ)ei ⇒ estabilidade para
−2 ≤ hλ ≤ 0. Este me´todo e´ enta˜o condicionalmente esta´vel.
• Euler impl´ıcito ei+1 − ei = hλei+1 ⇒ ei+1 = 1(1−hλ)ei ⇒ estabilidade para hλ ≤ 0.
Este me´todo e´ considerado absolutamente esta´vel.
• Trape´zio ei+1 − ei = hλ2 (ei+1 + ei) ⇒ ei+1 =
(1+hλ
2
)
(1−hλ
2
)
ei ⇒ estabilidade para hλ ≤ 0.
Este me´todo e´ considerado absolutamente esta´vel.
• Ponto me´dio ei+2 − 2hλei+1 − ei = 0 ⇒ as ra´ızes sa˜o: ζ1,2 = hλ ±
√
1 + (hλ)2 e
teremos estabilidade para | ζ1,2 |< 1.
• Fo´rmula 3(1 + 2λh)ei+2 − 4ei+1 + ei = 0⇒ as ra´ızes sa˜o: ζ1,2 = 2±
√
1−6hλ
3+6hλ e teremos
estabilidade para | ζ1,2 |< 1.
Observac¸a˜o 1.3.1 Um me´todo absolutamente esta´vel permite qualquer tamanho de passo
h, enquanto um condicionalmente esta´vel impo˜e restric¸a˜o ao passo h.
Observac¸a˜o 1.3.2 Para -hλ crescente, o me´todo de Euler impl´ıcito apresenta amorteci-
mento e o do Trape´zio apresenta oscilac¸a˜o, pois o fator de amplificac¸a˜o do me´todo de Euler
1.3. PROBLEMA DE VALOR INICIAL EM EQUAC¸O˜ES ORDINA´RIAS 27
impl´ıcito e´
(1 − hλ)−1 → 0 e o do trape´zio e´ (1+hλ2 )
(1−hλ
2
)
→ −1 quando −hλ → ∞. Portanto, esses
me´todos sa˜o mais indicados, respectivamente, para problemas que apresentam soluc¸a˜o com
decaimento e para problemas com soluc¸a˜o perio´dica. No contexto das equac¸o˜es parciais, um
me´todo similar ao de Euler impl´ıcito seria adequado para equac¸o˜es parabo´licas e um similar
ao Trape´zio, mais adequado para equac¸o˜es hiperbo´licas (na forma de conservac¸a˜o).
Em resumo, nesta sec¸a˜o apresentamos os principais conceitos e delineamos as maiores
dificuldades dos me´todos nume´ricos e quais as ferramentas para diagnostica´-las. Na pra´tica
pore´m, podem acontecer dificuldades extras que sa˜o inerentes ao problema. Nesse caso,
me´todos esta´veis gerais podem na˜o ser robustos o suficiente para fazer frente ao problema
e e´ recomenda´vel a construc¸a˜o de me´todos especiais, produzidos sob encomenda, que levemem conta a natureza do problema.
Exemplos
Nesta sec¸a˜o apresentamos exemplos nume´ricos ilustrando o efeito pra´tico de cada um dos
conceitos de consisteˆncia, zero-estabilidade e estabilidade.
Seja o PVI:
y′ = 4x
√
y, y(0) = 1, 0 ≤ x ≤ 2.
cuja soluc¸a˜o exata e´ y(x) = (1 + x2)2. Consideremos primeiramente a soluc¸a˜o desse PVI
pelo seguinte me´todo expl´ıcito de 2 passos:
yn+2 = yn+1 +
h
3
(3fn+1 − 2fn) (1.30)
que e´ zero-esta´vel, mas na˜o e´ consistente, como pode ser facilmente testado, ver exerc´ıcio
(1.10).
Ilustrac¸a˜o da Consisteˆncia
x |Erro| |Erro| |Erro|
h = 0.1 h = 0.05 h = 0.025
0 0 0 0
0.1 0 5.07E-03 8.82E-03
0.2 2.11E-02 3.61E-02 4.50E-02
0.3 7.21E-02 9.76E-02 1.12E-01
0.4 1.57E-01 1.95E-01 2.15E-01
0.5 2.85E-01 3.34E-01 3.61E-01
...
...
...
...
2.0 1.82E+01 1.89E+01 1.92E+01
Como o me´todo e´ de dois passos necessitamos duas condic¸o˜es iniciais para comec¸ar o processo
de ca´lculo, tomamos y0 = 1 e y1 = y(h) = (1 + h
2)2. Note que quando o passo diminui o
erro em y(2) aumenta mas na˜o de maneira explosiva. A sequeˆncia de soluc¸o˜es nume´ricas
converge para uma func¸a˜o y(x) a` medida que h→ 0, mas como o me´todo na˜o e´ consistente,
essa func¸a˜o na˜o e´ soluc¸a˜o do PVI dado, e sim de um outro problema.
28 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
Nosso pro´ximo exemplo ilustra as consequeˆncias desastrosas da falta da zero-estabilidade.
Seja o me´todo nume´rico expl´ıcito de dois passos:
yn+2 = (1 + a)yn+1 − ayn + h
2
((3− a)fn+1 − (1 + a)fn) (1.31)
Veja o exerc´ıcio (1.10) para concluir que esse me´todo e´ zero-esta´vel para a = 0, mas na˜o o
e´ quando a = −5. Assim:
Ilustrac¸a˜o da Zero-Estabilidade
x |Erro| |Erro| |Erro| |Erro| |Erro| |Erro|
a = 0 a = −5 a = 0 a = −5 a = 0 a = −5
h = 0.1 h = 0.1 h = 0.05 h = 0.05 h = 0.025 h = 0.025
0 0 0 0 0 0 0
0.1 0 0 5.63E-05 2.50E-05 2.23E-05 2.8E-05
0.2 9.00E-04 4.00E-04 3.61E-04 4.42E-04 1.09E-04 1.65E-02
0.3 2.85E-03 1.14E-03 9.38E-04 1.03E-02 2.64E-04 1.01E+02
0.4 5.97E-0.3 6.73E-3 1.82E-03 2.50E-01 4.95E-04 6.29E+3
0.5 1.04E-02 3.04E-02 3.05E-03 5.99E+00 8.17E-04 3.93E+6
...
...
...
...
...
...
...
2.0 3.67E-01 6.96E+08 9.63E-02 5.46E+21 4.26E-02 3.41E+48
Observe que no caso a = 0 ao refinarmos a malha as aproximac¸o˜es produzidas ficam mel-
hores, mas o mesmo na˜o acontece com a = −5. No primeiro caso o me´todo e´ convergente e
no segundo na˜o.
Finalmente vamos exemplificar o efeito da estabilidade absoluta no comportamento da
soluc¸a˜o nume´rica de longo prazo. Para isso consideremos o PVI y′ = −10(y− 1)2 y(0) = 2
cuja soluc¸a˜o exata e´ y(x) = 1+1/(1+ 10x). Resolvemos esse problema utilizando o me´todo
de Simpson
yn+2 − yn = h
3
(fn+2 + 4fn+1 + fn)
com h = 0.1 que e´ impl´ıcito de dois passos, mas na˜o e´ absolutamente esta´vel, veja o exerc´ıcio
(1.10).
Observe o crescimento do erro ao longo do processo.
Ilustrac¸a˜o dos efeitos da estabilidade
x 0.5 1.0 1.5 2.0 2.5
|Erro| 2.96E-02 5.64E-02 4.58E-02 9.59E-02 2.68E-01
x 3.0 3.5 4.0 4.8 4.9
|Erro| 1.28E-01 3.02E-01 1.73E-01 9.79E-01 Complexo
Estamos utilizando a palavra complexo na u´ltima coluna significando que o resultado pro-
duzido pelo me´todo de Simpson e´ complexo. Isto ocorre porque, sendo o me´todo de Simpson
impl´ıcito, devemos resolver uma equac¸a˜o do segundo grau para determinar a soluc¸a˜o. Essa
equac¸a˜o e´ resolvida exatamente, resultando em ra´ızes complexas. Na pra´tica, a resoluc¸a˜o
exata, na˜o seria poss´ıvel e ter´ıamos que resolveˆ-la por um me´todo iterativo, que teria se´rios
problemas de convergeˆncia, dado que a equac¸a˜o na˜o possui soluc¸a˜o real.
1.4. PROBLEMA DE VALOR DE FRONTEIRA EM EQUAC¸O˜ES ORDINA´RIAS 29
1.4 Problema de Valor de Fronteira em Equac¸o˜es Or-
dina´rias
No problema de valor inicial (PVI), como visto, uma condic¸a˜o inicial se propaga passo
a passo. No caso do problema de segunda ordem, foi preciso fixarmos uma condic¸a˜o inicial
e tambe´m uma inclinac¸a˜o (a derivada da soluc¸a˜o), ambas no mesmo ponto veja (1.15).
Para um problema de valor de fronteira (PVF) em equac¸o˜es ordina´rias e´ preciso fixar
valores da soluc¸a˜o em 2 pontos distintos, o que significa que a soluc¸a˜o propaga a condic¸a˜o
inicial em x = a, na direc¸a˜o do valor fixado em x = b, e vice-versa. E´ intuitivo que isso nem
sempre e´ poss´ıvel, ou seja, a condic¸a˜o de existeˆncia e unicidade de soluc¸a˜o e´ mais complicada
nestes casos. Por exemplo,
y′′ + π2y = 0, y(0) = 0, y(1) = 1 na˜o tem soluc¸a˜o,
y′′ + π2y = 0, y(0) = 0, y(1) = 0 tem infinitas soluc¸o˜es.
Um dos PVFs mais simples e de alguma importaˆncia pra´tica pode ser colocado como:
Sejam quatro constantes a, b, α e β, com a < b, e uma equac¸a˜o diferencial de segunda
ordem
y′′ = f(x, y, y′), (1.32)
o problema consiste em encontrar uma func¸a˜o y(x) ∈ C2[a, b], satisfazendo (1.32) em
a < x < b, e as duas condic¸o˜es de fronteira:
y(a) = α y(b) = β. (1.33)
Teorema 1.3 (Existeˆncia e Unicidade) Se o PVF (1.32) com as condic¸o˜es de fronteira
(1.33) e´ tal que, para −∞ < a < b < ∞, com α e β constantes arbitra´rias; a func¸a˜o f e´
Lipschitziana com relac¸a˜o a`s varia´veis y e y′, fy e´ cont´ınua e satisfaz fy > 0 e ainda fy′ e´
limitada para a ≤ x ≤ b e −∞ < y <∞, enta˜o ele tem soluc¸a˜o u´nica. [4]
O problema a ser resolvido pode apresentar dificuldades inerentes se ele for mal condi-
cionado. Por exemplo, considere o problema:
y′′ + k2y = 2k2x, y(0) = 0, y(1) = 1,
cuja soluc¸a˜o geral e´: y(x) = c1 sen kx+ c2cos kx + 2x. Das condic¸o˜es de contorno tem-se
c2 = 0 e c1 = −2/ sen k. Para k = π, o problema na˜o tem soluc¸a˜o. Para k = π + ε, ε
pequeno, c1 = −2/ sen ε ≃ 2ε−1, ou seja, uma pequena variac¸a˜o nos dados do problema
implica em grande variac¸a˜o na soluc¸a˜o que passou a existir. O problema acima constitui um
exemplo de um problema mal posto, ou mal condicionado.
Note que instabilidade e´ uma dificuldade do me´todo nume´rico, enquanto mal condiciona-
mento diz respeito ao problema.
A soluc¸a˜o do problema (1.32) pode ser obtida via me´todos para problema de valor inicial,
(Me´todo Shooting), que usa a condic¸a˜o dada em x = a e escolhe um valor para a derivada
(inclinac¸a˜o) tambe´m em x = a, calcula a soluc¸a˜o ate´ x = b e, de acordo com o resultado,
ajusta o valor da derivada em x = a e resolve novamente, repetindo estes passos ate´ que a
30 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
soluc¸a˜o propagando-se a partir de a encontre o valor de fronteira fixado em b. Na˜o vamos
detalhar este me´todo aqui, vide [Ames 1992]. Outra abordagem procura resolver o problema,
como problema de contorno; o me´todo, didaticamente mais recomendado e´ aquele que utiliza
diferenc¸as finitas, o qual detalharemos a seguir.
1.4.1 Me´todos de Diferenc¸as Finitas
Consideraremos aqui, Me´todos de Diferenc¸as Finitas para soluc¸a˜o de PVF que na˜o sa˜o
baseados na resoluc¸a˜o de problemas de valor inicial. O objetivo e´ fornecer um apanhado
geral, enfatizando os conceitos e dificuldades mais frequentes, como e´ o caso da restric¸a˜o do
tamanho h da malha e do tratamento da condic¸a˜o de fronteira com derivadas.
Seja o PVF
y′′(x) − p(x)y′(x)− q(x)y(x) = r(x)
y(a) = α, y(b) = β
que sera´ rescrito usando-se o operador L, como
L[y(x)] ≡ y′′(x) − p(x)y′(x)− q(x)y(x) = r(x), (1.34)
y(a) = α, y(b) = β.
Impomos a restric¸a˜o q(x) ≥ Q∗ > 0, a ≤ x ≤ b que e´ mais restritiva do que Q∗ ≥ 0
requerida para a prova de existeˆncia e unicidade de soluc¸a˜o .
Supomos no que segue, que o PVF objeto de estudo possui uma u´nica soluc¸a˜o com quatro
derivadas cont´ınuas em a ≤ x ≤ b. Usaremos uma malha uniforme com h = b−aN+1 e diferenc¸as
centradas para aproximar L[y(x)] para obter:
Lh[yj ] ≡ yj+1 − 2yj + yj−1
h2
−p(xj)yj+1 − yj−1
2h
−q(xj)yj = r(xj), j = 1, 2, . . . , N. (1.35)
As condic¸o˜es de fronteira sa˜o substitu´ıdas por:
y0 = α, yN+1 = β. (1.36)Multiplicando (1.35) por −h22 , obtemos:
−h
2
2
Lh[yj ] ≡ −1
2
{yj+1 − 2yj + yj−1}+ p(xj)h
4
{yj+1 − yj−1}+
q(xj)
h2
2
yj = −h
2
2
r(xj), j = 1, 2, . . . , N.
Rearranjando os termos de mesmo ı´ndice:
−h
2
2
Lh[yj ] ≡ −1
2
{
1 +
h
2
p(xj)
}
yj−1 +
{
1 +
h2
2
q(xj)
}
yj
1.4. PROBLEMA DE VALOR DE FRONTEIRA EM EQUAC¸O˜ES ORDINA´RIAS 31
−1
2
{
1− h
2
p(xj)
}
yj+1 = −h
2
2
r(xj), j = 1, 2, . . .N.
Fazendo:
aj ≡ 1 + h
2
2
q(xj),
bj ≡ 1
2
{
1 +
h
2
p(xj)
}
,
cj ≡ 1
2
{
1− h
2
p(xj)
}
,
chegamos a` equac¸a˜o
−h
2
2
Lh[yj ] ≡ −bjyj−1 + ajyj − cjyj+1 = −h
2
2
r(xj), j = 1, 2, . . .N.
que em notac¸a˜o vetorial, juntamente com as condic¸o˜es de fronteira (1.36), pode ser escrita
na forma:
Ay = r (1.37)
onde
A ≡


a1 −c1 0 · · · 0
−b2 a2 −c2 · · · 0
...
...
...
...
...
0 · · · −bN−1 aN−1 −cN−1
0 · · · 0 −bN aN

 , (1.38)
y ≡


y1
y2
...
yN−1
yN

 e r ≡ −
h2
2


r(x1)− 2b1αh2
r(x2)
...
r(xN−1)
r(xN )− 2cNβh2

 .
Assim, resolvendo o sistema linear de ordem N (1.37), cuja matriz dos coeficientes,
A, e´ tridiagonal, determinamos os valores de yj, j = 1, 2, . . . , N , que desejamos sejam
aproximac¸o˜es da soluc¸a˜o exata, y(xj), onde xj = x0 + jh.
Tendo constru´ıdo o me´todo, deve-se procurar conhecer a consisteˆncia e convergeˆncia do
mesmo. Um ponto novo aqui e´ que deve-se tambe´m cuidar da existeˆncia e unicidade da
soluc¸a˜o nume´rica. Se A e´ invers´ıvel tem-se soluc¸a˜o u´nica. Se A for quase singular, o me´todo
sera´ na˜o esta´vel. Note que isso pode ocorrer quando h → 0 e a dimensa˜o de A cresce.
Quando A e´ quase singular, os elementos de A−1 crescem com o aumento da ordem de
A, deteriorando o erro e consequentemente a convergeˆncia. Mais uma vez, estabilidade e
consisteˆncia equivalem a convergeˆncia.
32 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
Um crite´rio a ser considerado, e´ criar condic¸o˜es para que A seja diagonalmente domi-
nante. Vamos exigir que o espac¸amento h seja ta˜o pequeno que
h
2
|p(xj)| ≤ 1, j = 1, 2, . . . , N. (1.39)
Enta˜o, segue que
|bj|+ |cj | = bj + cj = 1 e aj > 1.
Logo,
|a1| > |c1|;
|aj | > |bj|+ |cj |, j = 2, 3, . . . , N − 1;
|aN | > |bN |;
e, portanto, a matriz A e´ diagonalmente dominante, por conseguinte invers´ıvel e o problema
discreto tem soluc¸a˜o u´nica. Na linguagem dos me´todos nume´ricos para problemas de valor
inicial a condic¸a˜o (1.39) representa uma condic¸a˜o de estabilidade. Note que quanto mais
acentuada for a dominaˆncia da diagonal de A, mais esta´vel sera´ o me´todo.
A soluc¸a˜o de (1.37) pode, assim, ser obtida por simples fatorac¸a˜o de A (Gauss, por
exemplo). Sendo A diagonalmente dominante, a fatorac¸a˜o e´ esta´vel e na˜o e´ necessa´rio
pivotamento. E´ claro que isto fornece uma prova da existeˆncia e unicidade da soluc¸a˜o
nume´rica de (1.37), desde que (1.39) seja satisfeita.
Vamos, agora, estimar o erro na aproximac¸a˜o nume´rica definida em (1.35). O erro
de truncamento local e´ obtido substituindo-se a soluc¸a˜o exata na equac¸a˜o (1.35), assim
definimos τj , o erro de truncamento local, por:
Lh[y(xj)] − r(xj) = τj , j = 1, 2, . . . , N. (1.40)
Como y(x) e´ soluc¸a˜o de (1.34) temos,
τj = Lh[y(xj)]− L[y(xj)]
=
y(xj − h)− 2y(xj) + y(xj + h)
h2
− p(xj)y(xj + h)− y(xj − h)
2h
−q(xj)y(xj)− [y′′(xj)− p(xj)y′(xj)− q(xj)y(xj)]
=
[
y(xj − h)− 2y(xj) + y(xj + h)
h2
− y′′(xj)
]
−p(xj)
[
y(xj + h)− y(xj − h)
2h
− y′(xj)
]
, j = 1, 2, . . . , N. (1.41)
Supondo que y(4)(x) e´ cont´ınua e usando a expressa˜o do erro nas fo´rmulas centradas
obtemos:
1.4. PROBLEMA DE VALOR DE FRONTEIRA EM EQUAC¸O˜ES ORDINA´RIAS 33
τj =
h2
12
[y(4)(ξj)− 2p(xj)y′′′(ηj)], j = 1, 2, . . . , N, (1.42)
onde ξj ,ηj ∈ [xj−1, xj+1]. Assim o erro de truncamento local e´ O(h2) e quando h → 0
o me´todo tende para a equac¸a˜o diferencial com ordem O(h2). Dizemos que o me´todo e´
consistente de ordem 2.
A estimativa de erro global e´ obtida do seguinte resultado:
Teorema 1.4 Se o espac¸amento h da malha satisfaz (1.39) enta˜o,
|y(xj)− yj | ≤ h2
(
M4 + 2P
∗M3
12Q∗
)
, j = 0, 1, . . . , N + 1,
onde y(x) e´ a soluc¸a˜o exata de (1.34), yi e´ a soluc¸a˜o calculada de (1.35) e
P ∗ ≡ maxx∈[a,b] |p(x)|, M3 ≡ maxx∈[a,b] |y′′′(x)|, M4 ≡ maxx∈[a,b]
∣∣y(4)(x)∣∣.
Prova: Sejam ej = y(xj)− yj , j = 0, 1, . . . , N + 1, ve: vetor dos erros ej e vτ : vetor
dos erros τj . Subtraindo (1.35) de (1.40) obtemos:
Lh[y(xj)]− Lh[yj ] = r(xj)− r(xj) + τj = τj , j = 1, 2, . . . , N.
Logo,
Lh[y(xj)− yj] = τj , ou Lh[ej ] = τj .
Esta u´ltima expressa˜o nos diz que o erro global satisfaz a mesma equac¸a˜o de diferenc¸as que
define o me´todo, tendo como lado direito o erro de truncamento local.
Multiplicando ambos os membros da expressa˜o acima por −h22 obtemos:
−h
2
2
Lh[ej] = −h
2
2
τj , j = 1, 2, . . . , N, (1.43)
que utilizando a mesma notac¸a˜o matricial de (1.37) produz, Ave = −h22 vτ , donde segue
que ve = −h22 A−1vτ , e mais uma vez observa-se a necessidade dos elementos de A−1 na˜o
crescerem (estabilidade). Os vetores ve e vτ sa˜o formados respectivamente pelas compo-
nentes do erro global e do erro de truncamento local.
Escrevendo individualmente cada equac¸a˜o de (1.43) obtemos
−bjej−1 + ajej − cjej+1 = −h
2
2
τj , j = 1, 2, . . . , N.
Da´ı,
ajej = bjej−1 + cjej+1 − h
2
2
τj , j = 1, 2, . . . , N
|ajej | ≤ |bjej−1|+ |cjej+1|+ h
2
2
|τj |, j = 1, 2, . . . , N
|ajej | ≤ (|bj |+ |cj |)e+ h
2
2
|τj |, j = 1, 2, . . . , N
34 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
e usando (1.39)
|ajej | ≤ e+ h
2
2
τ,
onde:
e ≡ max
0≤j≤N+1
|ej | e τ ≡ max
1≤j≤N
|τj |.
Entretanto, tambe´m pela condic¸a˜o (1.39) segue que,
|aj | = aj ≥ 1 + h
2
2
Q∗.
Logo, (
1 +
h2
2
Q∗
)
|ej| ≤ |ajej | ≤ e+ h
2
2
τ, j = 1, 2, . . . , N.
De (1.36) temos e0 = eN+1 = 0. Portanto, a desigualdade acima vale para todo j em
0 ≤ j ≤ N + 1.
Assim, conclu´ımos que (
1 +
h2
2
Q∗
)
e ≤ e+ h
2
2
τ
ou seja,
h2
2
Q∗e ≤ h
2
2
τ.
Da´ı,
e ≤ 1
Q∗
τ.
Finalmente, usando as hipo´teses do teorema em (1.42) encontramos que,
|τj | ≤ h
2
12
(M4 + 2P
∗M3) , j = 1, 2, . . . , N.
Logo,
τ ≤ h
2
12
(M4 + 2P
∗M3) ,
e portanto,
e ≤ h2
(
M4 + 2P
∗M3
12Q∗
)
.2
Do teorema 1.4 vemos que, sob a condic¸a˜o (1.39), a soluc¸a˜o discreta converge para a
soluc¸a˜o exata quando h→ 0, e ainda, que o erro e´ O(h2).
Muitas vezes, na pra´tica, esta ana´lise e´ dif´ıcil ou imposs´ıvel, calcula-se enta˜o a soluc¸a˜o
com va´rios tamanhos de h, observando-se o comportamento das soluc¸o˜es obtidas. Isto
permite-nos ter uma ide´ia da convergeˆncia e precisa˜o dos resultados (nu´mero de d´ıgitos
fixados) e tambe´m da taxa de convergeˆncia. Por exemplo, se o erro e´ O(h2), recalcular o
problema com nova malha de espac¸amento h2 deve resultar na duplicac¸a˜o da precisa˜o.
Alguns aspectos interessantes podem ainda ser considerados. Por exemplo:
1.4. PROBLEMA DE VALOR DE FRONTEIRA EM EQUAC¸O˜ES ORDINA´RIAS 35
i. outro tipo de condic¸a˜o de fronteira impo˜e alguma dificuldade?
ii. e´ poss´ıvel relaxar a condic¸a˜o (1.39) permitindo o uso de h maior?
iii. e´ possivel obter me´todos mais precisos?
iv. ocorre alguma dificuldade no caso da equac¸a˜o ser na˜o linear?
Vamos tratar de alguns desses aspectos a seguir. Considere que no lugar de y(a) = α em
(1.36), a condic¸a˜o dada seja y′(a) + γy(a) = 0. Essa nova condic¸a˜o tem consequeˆncias no
sistema (1.37) que passamos a analisar.
Uma primeira abordagem e´ tomar uma aproximac¸a˜o progressiva (de ordem O(h)) para
a derivada, ou seja,
y′(a) ≃ y(a+ h)− y(a)
h=
y1 − y0
h
donde
y1 − y0 + hγy0 = 0⇒ y0 = y1
1− hγ
que incorporada ao sistema, modifica apenas a primeira equac¸a˜o para
(a1 − b1
1− hγ )y1 − c1y2 = −
h2
2
r(x0)
deixando o restante inalterado. O inconveniente aqui e´ o fato de introduzirmos um erro
local de O(h) que deteriora um pouco o resultado final.
Uma outra abordagem, que recupera a ordem 2, se da´ pelo uso de uma aproximac¸a˜o
central
y′(a) ≃ y(a+ h)− y(a− h)
2h
=
y1 − y−1
2h
onde y−1 e´ introduzido artificialmente.
Assim,
y1 − y−1 + 2hγy0 = 0.
Como uma nova inco´gnita foi introduzida, nova equac¸a˜o e´ necessa´ria. Esta pode ser
providenciada considerando-se como primeira equac¸a˜o em (1.37), a equac¸a˜o obtida para
j = 0
−b0y−1 + a0y0 − c0y1 = −h
2
2
r(x0),
que por substituic¸a˜o de y−1 fica
(a0 − b02hγ)y0 − (b0 + c0)y1 = −h
2
2
r(x0).
36 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
O sistema (1.37) fica agora com uma equac¸a˜o extra, entrando y0 como inco´gnita. Apesar
de pequeno acre´scimo no esforc¸o computacional, o resultado final e´ O(h2), conservando-se a
ordem anterior.
Com relac¸a˜o a` utilizac¸a˜o de tamanhos de passo maiores, e´ poss´ıvel a construc¸a˜o de
um me´todo menos restritivo. Esse me´todo denominado upwind, tem bastante aplicac¸a˜o na
pra´tica, principalmente em equac¸o˜es hiperbo´licas e consiste em combinarmos as diferenc¸as
progressiva e regressiva conforme o sinal de p(x). Isto reduz a ordem do me´todo, mas a
restric¸a˜o sobre h e´ relaxada.
Assim, o me´todo upwind e´ baseado na seguinte aproximac¸a˜o para a primeira derivada:
y′(xj) ≃ η yj+1 − yj
h
+ (1− η)yj − yj−1
h
. (1.44)
Portanto para η = 1, 12 e −1 tem-se as aproximac¸o˜es progressiva, central e regressiva, re-
spectivamente.
Usando (1.44) em (1.35) temos
−h
2
2
Lh[yj ] ≡ −1
2
{yj+1 − 2yj + yj−1}+
p(xj)h {η(yj+1 − yj) + (1− η)(yj − yj−1)}+
q(xj)
h2
2
yj = −h
2
2
r(xj), j = 1, 2, . . . , N
donde
−bjyj−1 + ajyj − cjyj+1 = −h
2
2
r(xj), j = 1, 2, . . .N
com
aj ≡ 1− (1 − 2η)hp(xj)− h2q(xj),
bj ≡ 1
2
{1− (1− η)hp(xj)} ,
cj ≡ 1
2
{1 + hηp(xj)} .
No me´todo upwind, η vai assumir valores 0 ou 1 , conforme p(xj) < 0 ou p(xj) > 0,
respectivamente. Isso mante´m a dominaˆncia da diagonal de A sem a restric¸a˜o (1.39). A
ordem do me´todo diminui um pouco, mas o esforc¸o computacional pode diminuir significa-
tivamente, porque a malha pode ser menos fina, ou seja, h pode ser maior.
Foi visto anteriormente que e´ poss´ıvel a melhoria da ordem de aproximac¸a˜o, aumentando
os pontos na aproximac¸a˜o da derivada. Esta abordagem acarreta um aumento na banda da
matriz do sistema e o esforc¸o computacional cresce. Para ilustrar uma outra maneira de se
melhorar a ordem, vamos considerar o problema
y′′ = f(x, y), y(0) = α y(1) = β.
1.5. EXERCI´CIOS 37
A discretizac¸a˜o convencional fornece:
yj+1 − 2yj + yj−1 = h2f(xj , yj)
que e´ consistente de ordem O(h2). Como esta˜o envolvidos yj+1, yj e yj−1, pode-se usa´-los
para se obter uma aproximac¸a˜o melhor para f(xj , yj) da seguinte maneira:
yj+1 − 2yj + yj−1 = h2{c1f(xj−1, yj−1) + c2f(xj , yj) + c3f(xj+1, yj+1)}
onde c1, c2 e c3 sa˜o calculados para aumentar a precisa˜o. O me´todo mais conhecido, assim
obtido, e´ o me´todo de Numerov
yj+1 − 2yj + yj−1 = h
2
12
{f(xj−1, yj−1) + 10f(xj, yj) + f(xj+1, yj+1)}
que e´ consistente de ordem O(h4).
Quando f e´ na˜o linear e c3 6= 0, esse me´todo requer a soluc¸a˜o de uma equac¸a˜o (ou
um sistema) na˜o-linear. Isto pode ser resolvido por iterac¸o˜es sucessivas ou pelo me´todo de
Newton, mas claramente o esforc¸o computacional e´ bem maior.
Em geral, um problema na˜o-linear de segunda ordem pode ser colocado como
y′′ = f(x, y, y′), y(0) = α y(1) = β
e a construc¸a˜o do me´todo segue o mesmo procedimento
yj+1 − 2yj + yj−1 = h2f(xj , yj , yj+1 − yj−1
2h
), j = 1, 2, . . . , N
com a necessidade de soluc¸a˜o de um sistema na˜o linear.
1.5 Exerc´ıcios
1.1 Utilizando propriedades dos operadores e expansa˜o formal em se´rie de poteˆncias das
func¸o˜es envolvidas prove as relac¸o˜es abaixo:
∇ = E−1∆, ∆∇ = ∇∆ = ∆−∇ = δ2, µδ = 1
2
(∆ +∇), µ2 = 1 + 1
4
δ2
∇ = 1− E−1, µ = 1
2
(
E
1
2 + E−
1
2
)
, δ = E
1
2 − E− 12 ,
1 = µµ−1 = µ
(
1 +
1
4
δ2
)− 1
2
= µ
(
1− 1
8
δ2 +
3
128
δ4 + · · ·
)
hD = ln(1 + ∆) = − ln(1−∇) = 2 sinh−1
(
δ
2
)
= 2 ln
(√
1 +
1
4
δ2 +
1
2
δ
)
38 CAPI´TULO 1. CONCEITOS FUNDAMENTAIS
(hD)k = [ln(1 + ∆)]
k
= ∆k − k
2
∆k+1 +
k(3k + 5)
24
∆k+2
− k(k + 2)(k + 3)
48
∆k+3 + · · ·
(hD)k = (−1)k [ln(1 −∇)]k = ∇k + k
2
∇k+1 + k(3k + 5)
24
∇k+2
+
k(k + 2)(k + 3)
48
∇k+3 + · · ·
hD =
2µ sinh−1
(
δ
2
)
(
1 + δ
2
4
) 1
2
= µ
[
δ − 1
2
3!
δ3 +
12.22
5!
δ5 − · · ·
]
(hD)2n =
[
2 sinh−1
(
δ
2
)]2n
=
[
δ − 1
2
22.3!
δ3 +
12.32
24.5!
δ5 − 1
2.32.52
26.7!
δ7 + · · ·
]2n
(hD)2n =
µ√
1 + δ
2
4
[
2 sinh−1
(
δ
2
)]2n
= µ
[
δ − 1
2
3!
δ3 +
12.22
5!
δ5 − · · ·
]2n
(hD)2n+1 =
[
2 sinh−1
(
δ
2
)]2n+1
=
[
δ − 1
24
δ3 +
3
640
δ5 − 5
7168
δ7 + · · ·
]2n+1
.
Sugesta˜o: Utilize as expanso˜es em se´rie:
sinh−1 (x) =
∞∑
k=0
(−1)k (2k)!
22k(k!)2(2k + 1)
x2k+1
(1 + x)
q
= 1 + qx+
q(q − 1)
2!
x2 + · · ·+ q(q − 1)(q − 2) · · · (q − k + 1)
k!
xk + · · ·
1.2 Utilizando expansa˜o em se´rie de Taylor de y′(xn+1) e comparando com a expansa˜o de
y(xn+1) deduza que:
y(xn+1) = y(xn) +
h
2
[f(xn, y(xn)) + f(xn+1, y(xn+1))]− h
3
12
y′′′(ξ).
1.3 As relac¸o˜es (1.6-1.10) sa˜o casos especiais das expanso˜es do exerc´ıcio (1.1), prove-as.
1.4 Deduzir a partir da expansa˜o em se´rie de Taylor de y(x), fo´rmulas de ordem 1 e de
ordem 2 para aproximar as derivadas y(3)(x) e y(4)(x).
1.5 Mostrar, usando a equac¸a˜o teste y′ = λy, que a regra do ponto me´dio e´ insta´vel.
1.6 a. Construa o me´todo de maior ordem poss´ıvel da forma:
a1yi+1 + a2yi + a3yi−1 + a4yi−2 = f(xi, yi)
para soluc¸a˜o de um PVI.
1.5. EXERCI´CIOS 39
b. Mostre que esse me´todo na˜o e´ convergente.
1.7 Utilizando MATLAB resolva o PVI abaixo pelo me´todo do ponto me´dio, Euler expl´ıcito,
impl´ıcito e regra dos trape´zios, no intervalo [0, 2], usando h = 0.1.
y′ = 5y − 1, y(0) = 1.2.
Explique o comportamento de cada um desses me´todos a` luz das propriedades discutidas
neste cap´ıtulo. Soluc¸a˜o exata: y(x) = exp(5x) + 0.2.
1.8 Mostre por induc¸a˜o sobre n que |en| ≤ ξn onde ξn e´ a soluc¸a˜o da equac¸a˜o (1.19), com
ξ0 = |e0|. Prove que (1.20) e´ realmente soluc¸a˜o de (1.19), com ξ0 = 0. A soluc¸a˜o (1.20) e´
va´lida se e somente se A e B forem constantes, no entanto tomamos n → ∞ em (1.21) o
que certamente faz com que h → 0 e portanto A e B na˜o sa˜o constantes. Como explicar
enta˜o a passagem ao limite em (1.21)?
1.9 Deduza a fo´rmula (1.28).
1.10 Mostre que o me´todo (1.30) e´ zero-esta´vel mas na˜o e´ consistente e que o me´todo
(1.31) e´ consistente para todo valor de a; zero-esta´vel e de ordem 2 se a = 0 e zero-insta´vel
e de ordem 3 se a = −5. (Veja [20] pa´gina 41). Mostre que o me´todo de Simpson na˜o e´
esta´vel.
1.11 a. Seja o problema de valor de fronteira de segunda ordem definido no intervalo
[a, b]:
y′′ − q(x)y = r(x) y(a) = y(b) = 0
Determinar os paraˆmetros α2, α1, α0 de modo que o me´todo nume´rico:
yi+1 − 2yi + yi−1
h2
− [α2q(xi+1)yi+1 + α1q(xi)yi + α0q(xi−1)yi−1] =
α2r(xi+1) + α1r(xi) + α0r(xi−1), i = 1, 2, . . . , N
com y0 = yN+1 = 0, e h =
b−a
N+1 , seja de ordem 4.
b. Se q(x) > Q∗ > 0 mostre que a expressa˜o
h4(2M + 5N + 5R)
720Q∗
e´ um limitante superior para o erro global do me´todo em (a) onde M e R sa˜o limitantes
superiores para a derivada de