Baixe métodos numéricos para problemas diferenciais ordinários e outras Manuais, Projetos, Pesquisas em PDF para Cálculo Numérico, somente na Docsity!
Cap´ıtulo 6
M´etodos num´ericos para problemas
diferenciais ordin´arios
6.1 Breve referˆencia hist´orica
Os trabalhos preliminares na ´area dos m´etodos num´ericos para a solu¸c˜ao de equa¸c˜oes di- ferenciais s˜ao devidos, entre outros, a Isaac Newton (1643-1729) e Gottfried Wilhelm Leibniz (1643-1716), no s´eculo XVII, bem como a Lehonard Euler (1707-1783), no s´eculo XVIII. Mas foi sobretudo devido aos trabalhos deste ´ultimo que foram impulsionados os estudos que conduziram aos m´etodos que hoje conhecemos. Euler deduziu um processo iterativo que permitia determinar, de forma aproximada, a solu¸c˜ao de um problema de condi¸c˜ao inicial num determinado ponto. A demonstra¸c˜ao rigo- rosa de que, de facto, o processo por ele apresentado conduzia a solu¸c˜ao da equa¸c˜ao s´o foi apresentada mais tarde, em 1824, por Augustine Cauchy (1789-1857) e melhorada por Rudolf Lipschitz (1832-1908). No entanto, nem assim, o processo apresentado por Euler se tornou popular. A t´ıtulo de exemplo, Karl Weierstrass (1815-1897), famoso matem´atico alem˜ao do s´eculo XIX trabalhou nestes assuntos sem conhecer os trabalhos de Cauchy e Lipschitz. Os finais do s´eculo XIX e princ´ıpios do s´eculo XX foram muito prof´ıquos ao florescimento de m´etodos num´ericos para a resolu¸c˜ao de equa¸c˜oes diferenciais. Com os desenvolvimentos efectuados na teoria do calor por Fourier e na mecˆanica celeste por Adams, Bessel, Cauchy, Gauss, Laguerre, Laplace, Legendre, Leverrier, Poincar´e e outros vieram tornar imprescind´ıvel a existˆencia de esquemas para determinar a solu¸c˜ao num´erica de equa¸c˜oes diferenciais. A bal´ıstica foi outra ciˆencia que come¸cou a exigir resultados nesta ´area. Poderemos dividir os m´etodos num´ericos para determinar a solu¸c˜ao de equa¸c˜oes diferen- ciais em duas grandes classes: por um lado, os chamados m´etodos de passo ´unico aos quais pertence o m´etodo de Euler-Cauchy-Lipschitz; por outro os chamados m´etodos de passo m´ultiplo. Os sucessores do m´etodo de Euler-Cauchy-Lipschitz foram apresentados por K. Heun em 1900 e, sobretudo, por Carl Runge (1856-1927) em 1895 e 1908 e por Martin Wilhelm Kutta (1867-1944) em 1901, tendo sido considerados como generaliza¸c˜oes das regras de integra¸c˜ao. Estes m´etodos tornaram-se bastante populares devidoas suas propriedades e `a sua f´acil utiliza¸c˜ao.^1 (^1) De notar que o primeiro sistema de equa¸c˜oes diferenciais a ser resolvido pelo ENIAC foi integrado pelo m´etodo de Heun.
Dos primeiros e mais conhecidos m´etodos de passo m´ultiplo para a resolu¸c˜ao de equa¸c˜oes diferenciais s˜ao os chamados m´etodos de Adams. John Couch Adams (1819-1892), famoso astr´onomo britˆanico que descobriu, em co-autoria, o planeta Neptuno, baseou-se nos m´etodos te´oricos propostos por Cauchy para apresentar um m´etodo novo que usou na integra¸c˜ao da equa¸c˜ao de Bashforth. Ali´as, foi num trabalho de Francis Bashforth (1819-1912) de 1883 que o m´etodo proposto por Adams foi apresentado, sendo por isso tamb´em conhecido por m´etodo de Adams-Bashforth. Foi possivelmente a primeira Grande Guerra que veio dar um forte impulso ao floresci- mento dos m´etodos num´ericos. A grande quantidade de c´alculos e a complexidade dos proble- mas que a bal´ıstica exige n˜ao poderiam ser efectuadas facilmente sem a ajuda destes processos alternativos. O primeiro contributo para o melhoramento dos m´etodos existentes foi dado pelo matem´atico americano Forest Ray Moulton (1872-1952) em 1925, propondo uma classe de m´etodos conhecida por Adams-Moulton. Em 1928 apareceu um trabalho da autoria de Richard Courant (1888-1972), Kurt Friedrichs (1901-1982) e Hans Lewy que revolucionou toda a teoria dos m´etodos num´ericos para a resolu¸c˜ao de equa¸c˜oes diferenciais. No entanto foi s´o depois da segunda Grande Guerra e sobretudo depois do aparecimento do primeiro computador e dos trabalhos de Herman Goldstine (1903-) e John von Neumann (1903- 1957), em 1947, que estes m´etodos come¸caram a ser usados de forma sistem´atica. Conceitos como ’erros de arredondamento’, ’n´umero de condi¸c˜ao’ e ’instabilidade num´erica’ come¸caram a surgir e a tornar-se de capital importˆancia para o estudo da teoria subjacente. Os estudos sobre m´etodos num´ericos para a resolu¸c˜ao de equa¸c˜oes diferenciais come¸caram a merecer a aten¸c˜ao de um n´umero crescente de matem´aticos e outros cientistas e hoje ´e uma das ´areas mais importantes e prof´ıquas da Matem´atica em geral e da An´alise Num´erica em particular.
6.2 Introdu¸c˜ao
As primeiras equa¸c˜oes diferenciais s˜ao t˜ao antigas quanto o c´alculo diferencial. Newton considerou-as, em 1671, no seu tratado de c´alculo diferencial e discutiu a sua solu¸c˜ao por integra¸c˜ao e por expans˜ao em s´erie. Leibniz, o segundo inventor do c´alculo, chegou as equ¸c˜oes diferenciais por volta de 1676 considerando o problema geom´etrico do ’inverso das tangentes’: para que curva y(x) a tangente em cada ponto P tem um comprimento constante (com o eixo dos x’s), digamos a? Este problema conduziua equa¸c˜ao y′^ = −y/
√ a^2 − y^2. Em 1696, Johann Bernoulli (1667-1748) convidou os mais ilustres matem´aticos do seu tempo para resolver o problema da braquist´ocrona (curva de tempo m´ınimo), principalmente para refutar a resposta, que esperava errada, do seu irm˜ao Jacob Bernoulli (1657-1705). O problema consistia em determinar a curva y(x) que une dois pontos P 0 e P 1 de tal modo que um ponto, partindo de P 0 e ’deslizando’, nessa curva, sujeito apenas a for¸cas grav´ıticas, atinja P 1 no menor tempo poss´ıvel. A resposta a este problema foi dada dada por v´arios matem´aticos (inclusiv´e Jacob Bernoulli) e ´e, como se sabe, a cicl´ıode. Essa curva pode ser determinada como sendo a solu¸c˜ao de uma equa¸c˜ao diferencial ordin´aria. Muitos problemas da engenharia e da ciˆencia tˆem como modelo equa¸c˜oes diferenciais. Neste curso iremos efectuar uma breve introdu¸c˜ao ao estudo dos m´etodos num´ericos para a resolu¸c˜ao de problemas que envolvem equa¸c˜oes diferenciais. Os problemas que iremos con- siderar ser˜ao de dois tipos: problemas com condi¸c˜ao inicial e problemas com condi¸c˜oes de fronteira. Neste cap´ıtulo n˜ao iremos apresentar a demostra¸c˜ao de todos os teoremas; o aluno interessado poder´a ver essas demonstra¸c˜oes nas obras referenciadas no final do cap´ıtulo.
Observa¸c˜ao 6.3 A demonstra¸c˜ao deste teorema estabelece um processo iterativo de aproxi- ma¸c˜ao da solu¸c˜ao do PCI (6.2) conhecido por m´etodo de Emile Picard (1856-1941). Se f for cont´ınua em rela¸c˜ao a t, determinar a solu¸c˜ao do PCI (6.2) ´e equivalente a determinar y, continuamente diferenci´avel, que verifica
y(t) = α +
∫ (^) t
a
f (τ, y(τ ))dτ. (6.3)
O que se prova na demonstra¸c˜ao do Teorema de Picard ´e que a sucess˜ao de fun¸c˜oes (yj (t)) definida recursivamente por
y 0 (t) = α, yj+1 = α +
∫ (^) t a f^ (τ, yj^ (τ^ )dτ.^ j^ = 0,^1 ,... ,
converge para a ´unica solu¸c˜ao de (6.3).
Como corol´ario do Teorema de Picard temos o seguinte resultado que apresentamos igual- mente sem demonstra¸c˜ao.
Corol´ario 6.4 Suponhamos que f (t, y) est´a definida num conjunto convexo^2 D ⊂ IR^2. Se existir uma constante L > 0 tal que ∣∣ ∣∣^ ∂f ∂y
(t, y)
∣∣ ∣∣ ≤ L, ∀(t, y) ∈ D,
ent˜ao f satisfaz a condi¸c˜ao de Lipschitz, na vari´avel y, com L a respectiva constante e, como tal, o PCI (6.2) tem solu¸c˜ao ´unica y(t) para t ∈ [a, b].
Observa¸c˜ao 6.5 Note-se que o conjunto D = {(t, y) : a ≤ t ≤ b, y ∈ IR} ´e, obviamente, convexo.
Exerc´ıcio 6.3.2 Mostre que o problema de condi¸c˜ao inicial
y′(t) =
1 + y^2
y(a) = α,
para t ∈
[a, b], tem solu¸c˜ao ´unica.
Resolu¸c˜ao: Seja D = {(t, y) : a ≤ t ≤ b, y ∈ IR} e
f (y) =
1 + y^2
Vamos provar que a fun¸c˜ao (^) ∣ ∣∣ ∣
∂f ∂y (t, y)
∣∣ ∣∣ =
∣∣ ∣∣^ −^2 y (1 − y^2 )^2
∣∣ ∣∣
´e limitada em D. Para isso h´a que determinar
L = max y∈IR
∣∣ ∣∣^2 y (1 − y^2 )^2
∣∣ ∣∣.
(^2) Um conjunto D ⊆ R (^2) diz-se convexo se, para qualquer (t 1 , y 1 ), (t 2 , y 2 ) ∈ D, se verifica
((1 − θ)t 1 + θt 2 , (1 − θ)y 1 + θy 2 ) ∈ D, θ ∈ [0, 1].
Como a fun¸c˜ao que queremos provar limitada ´e par temos que
L = max y∈IR+ 0
2 y (1 − y^2 )^2
Consideremos g(y) = 2 y (1 − y^2 )^2
Como g′(y) = 0 ⇒ y = ±
temos que
L = max
{ g(0), g
) , lim y→+∞ g(y)
} = max{ 0 , 0. 6594 , 0 } = 0. 6594.
Est´a assim provado o pretendido.
Outra quest˜ao que se coloca ´e a de saber se o PCI (6.2) ´e sens´ıvel a pequenas perturba¸c˜oes nos dados do problema, isto ´e, se pequenas altera¸c˜oes nas condi¸c˜oes iniciais n˜ao provocam grandes altera¸c˜oes nos resultados. Jacques Hadamard, em 1923, sugeriu duas condi¸c˜oes que devem ser verificadas quando se formula um PCI: (i) a solu¸c˜ao deve existir e ser ´unica; (ii) a solu¸c˜ao deve depender de forma cont´ınua dos dados iniciais do problema. Temos ent˜ao a seguinte defini¸c˜ao.
Defini¸c˜ao 6.6 O PCI (6.2) ´e bem posto se:
- possui solu¸c˜ao ´unica y(t);
- for est´avel, isto ´e, se existirem constantes positivas ξ e k (independente de ξ) tais que o problema perturbado { z′(t) = f (t, z) + δ(t), t ∈ [a, b] z(a) = α + δ 0
possui solu¸c˜ao ´unica z(t) com
|z(t) − y(t)| < kξ, ∀t ∈ [a, b],
sempre que |δ 0 | < ξ e |δ(t)| < ξ.
As hip´oteses do Corol´ario 6.4 s˜ao suficientes para garantir que o PCI (6.2) ´e bem posto no sentido de Hadamard. De facto, o referido corol´ario j´a estabelece a existˆencia e unicidade e solu¸c˜ao; falta apenas provar que o problema ´e est´avel. Considerando
w(t) = z(t) − y(t),
temos w′(t) = f (t, z(t)) − f (t, y(t)) + δ(t) = fy(t, y(t))w(t) + δ(t),
com a condi¸c˜ao inicial w(a) = δ 0 ,
Os m´etodos num´ericos permitem determinar valores yi ≈ y(ti) por meio de rela¸c˜oes de recorrˆencia deduzidas do PCI (6.2) de modo a que o valor de yi+1 venha expresso em fun¸c˜ao de yi, yi− 1 ,... , y 0 , sendo y 0 = y(a) = α. E usual agrupar os m´´ etodos num´ericos para a resolu¸c˜ao de problemas de condi¸c˜ao inicial em duas grandes classes.
- M´etodos de passo ´unico: S˜ao m´etodos que determinam o valor de yi+1 apenas `a custa de yi.
- M´etodos de passo m´ultiplo: S˜ao m´etodos que determinam o valor de yi+1 `a custa de yi, yi− 1 ,... , yi−r+1. Neste caso diz-se que o m´etodo ´e de r passos.
Neste curso iremos apenas abordar os m´etodos de passo ´unico. Estes m´etodos, por sua vez, podem ainda ser de dois tipos.
- M´etodos expl´ıcitos: S˜ao m´etodos em que o valor de yi+1 ´e determinado directamente a partir de yi. Estes m´etodos podem ser escritos na forma
yi+1 = yi + hφ(ti, yi, ; h). (6.5)
- M´etodos impl´ıcitos: S˜ao m´etodos em que o valor de yi+1 depende implicitamente de si mesmo atrav´es de f. Estes m´etodos podem ser escritos na forma
yi+1 = yi + hφ(ti, ti+1, yi, yi+1; h). (6.6)
A fun¸c˜ao φ que define os m´etodos (6.5) e (6.6) ´e chamada fun¸c˜ao de itera¸c˜ao ou fun¸c˜ao incremento do m´etodo num´erico.
6.3.3 M´etodos baseados na s´erie de Taylor
Consideremos o PCI (6.2) com f uma fun¸c˜ao suficientemente diferenci´avel nas vari´aveis t e y. Ent˜ao
y(t) = y(t 0 ) + (t − t 0 )y′(t 0 ) +
(t − t 0 )^2 2! y′′(t 0 ) + · · · ,
com t 0 = a. As derivadas desta express˜ao n˜ao s˜ao conhecidas explicitamente visto que a solu¸c˜ao tamb´em n˜ao ´e conhecida. No entanto, podemos escrever
y′(t) = f (t, y),
y′′(t) = df dt
(t, y) = (ft + fyy′)(t, y) = (ft + fyf )(t, y),
y′′′(t) = df dt^2 (t, y) = (ftt + 2ftyf + fyyf 2 + ftfy + f (^) y^2 f )(t, y),
.. .
onde
ft(t, y) =
∂f ∂t (t, y), fy(t, y) =
∂f ∂y (t, y),....
Por raz˜oes pr´aticas temos que limitar o n´umero de termos na expans˜ao em s´erie de y(t) a um n´umero razo´avel, o que nos conduz a restri¸c˜oes nos valores de t para os quais a expans˜ao nos d´a uma boa aproxima¸c˜ao. Se tomarmos a s´erie de Taylor truncada temos, para t = t 1 ,
y(t 1 ) ≈ y 1 = y 0 + hf (t 0 , y 0 ) + h^2 2 f ′(t 0 , y 0 ) + · · · + hk k! f (k−1)(t 0 , y 0 ),
onde
f (j)(t 0 , y 0 ) = dj^ f dtj^
(t 0 , y 0 ).
Podemos definir assim, para cada k = 1, 2.. ., um m´etodo de passo ´unico expl´ıcito que permite obter solu¸c˜oes aproximadas yi ≈ y(ti) da forma (6.5) em que
φ(t, y; h) = f (t, y) +
h 2 f ′(t, y) + · · · +
hk k! f (k−1)(t, y). (6.7)
Os m´etodos assim definidos s˜ao conhecidos por m´etodos de Taylor. O m´etodo desta classe mais simples ´e quando k = 1, isto ´e, o m´etodo
yi+1 = yi + hf (ti, yi), i = 0,... , n, y 0 = α, (6.8)
designado por m´etodo de Euler (expl´ıcito).
O seguinte algoritmo permite determinar a solu¸c˜ao do PCI (6.2) em t = b, usando o m´etodo com fun¸c˜ao incremento (6.7).
Algoritmo 6.1 M´etodo de Taylor
Ler n e k; Ler a, b e α; h := b− na ; t := a; y := α; Para i de 1 at´e n fazer φ := 0; Para j de 1 at´e k fazer φ := φ + f (j)(t, y)hj^ /j!; y := y + hφ; t := t + h; Escrever y.
Exerc´ıcio 6.3.5 Considere o problema de condi¸c˜ao inicial
{ y′(t) = − 2 y, y(0) = 1. Determine, u-
sando o m´etodo de Euler, o valor aproximado de y(1), fazendo h = 1, h = 0. 5 e h = 0. 25. Compare os resultados obtidos sabendo que y(t) = e−^2 t.
se pretendermos uma precis˜ao de, digamos, 6 casas decimais, o m´etodo de Euler necessita de aproximadamente um milh˜ao de passos. Se usarmos outros m´etodos de Taylor, a precis˜ao pode ser aumentada. A grande desvan- tagem destes m´etodos reside no facto de termos necessidade de calcular muitas derivadas da fun¸c˜ao f para obter m´etodos precisos. Esse c´alculo, al´em de muito fastidioso, torna im- pratic´avel a aplica¸c˜ao de tais m´etodos na resolu¸c˜ao de (6.2) quando a fun¸c˜ao f tem uma express˜ao anal´ıtica complicada. Uma alternativa a esses m´etodos foi dada por Runge, em 1875, e que consistia em, partindo do conhecimento de y(a), considerar
y(a + h) ≈ α + hf
( a + h 2
, y
( a + h 2
)) ;
mas, que valor atribuir a y
( a + h 2
) ? A sugest˜ao de Runge foi a de considerar o m´etodo de
Euler com passo h 2. A aplica¸c˜ao sucessiva deste processo permitiu a Runge definir o seguinte m´etodo iterativo: k 1 = f (ti, yi), k 2 = f
( ti + h 2 , yi + h 2 k 1
) , yi+1 = yi + hk 2.
Como veremos este m´etodo, apesar de recorrer ao m´etodo de Euler, vai ser mais preciso e n˜ao necessita de calcular derivadas de f. A generaliza¸c˜ao desta ideia deu origem `a seguinte defini¸c˜ao.
Defini¸c˜ao 6.8 Seja s um n´umero inteiro e a 2 , 1 , a 3 , 1 , a 3 , 2 ,... , as, 1 ,... , as,s− 1 , c 2 , c 3 ,... , cs, b 1 , b 2 ,... , bs, coeficientes reais. O m´etodo
k 1 = f (ti, yi), k 2 = f (ti + c 2 h, yi + a 2 , 1 hk 1 ), k 3 = f (ti + c 3 h, yi + a 3 , 1 hk 1 + a 3 , 2 hk 2 ), .. . ks = f (ti + csh, yi + as, 1 hk 1 + as, 2 hk 2 + · · · + as,s− 1 hks− 1 ), yi+1 = yi + h[b 1 k 1 + b 2 k 2 + · · · + bsks],
´e chamado m´etodo de Runge-Kutta expl´ıcito de s etapas para o PCI (6.2).
Usualmente considera-se
ci =
∑^ i−^1 j=
ai,j , i = 2, 3 ,... , s. (6.10)
Uma nota¸c˜ao muito usada na pr´atica para os m´etodos de Runge-Kutta foi apresentada por Butcher em 1964 e ´e dada pelo seguinte quadro, designado por quadro de Butcher:
0 c 2 a 2 , 1 c 3 a 3 , 1 a 3 , 2 .. .
cs as, 1 as, 2 · · · as,s− 1 b 1 b 2 · · · bs− 1 bs
Antes de continuarmos, notemos que os m´etodos de Runge-Kutta constituem uma exce- lente ideia. A ´unica solu¸c˜ao do PCI bem posto (6.2) ´e uma curva integral em IR^2. No entanto, devido aos erros cometidos, a solu¸c˜ao num´erica vai ser afectada pelo comportamento das cur- vas integrais vizinhas. E assim importante conhecer o comportamento de toda a fam´´ ılia de curvas integrais e n˜ao apenas o de uma ´unica curva. Os m´etodo de Runge-Kutta usam, deliberadamente, informa¸c˜ao de v´arias curvas integrais em simultˆaneo. A t´ıtulo de exemplo considere-se o m´etodo de trˆes etapas
k 1 = f (ti, yi), k 2 = f (ti + c 2 h, yi + c 2 hk 1 ), k 3 = f (ti + c 3 h, yi + (c 3 − a 3 , 2 )hk 1 + a 3 , 2 hk 2 ), yi+1 = yi + h[b 1 k 1 + b 2 k 2 + +b 3 k 3 ].
Para determinar a solu¸c˜ao num´erica do PCI (6.2) por este m´etodo, come¸ca-se pelo ponto (ti, yi) e aplica-se um passo do m´etodo de Euler com passo c 2 h. Seguidamente, calcula-se o valor de k 2 como sendo o vector derivada no ponto obtido. Temos assim dois valores para a derivada: k 1 e k 2 ; iremos usar uma m´edia pesada entre estes dois valores (c 3 −a 3 , 2 )hk 1 +a 3 , 2 hk 2 numa nova aplica¸c˜ao do m´etodo de Euler, a partir do ponto (ti, yi), com passo c 3 h. Calculando a derivada novamente obt´em-se o valor de k 3. O ´ultimo passo do algoritmo ´e mais uma aplica¸c˜ao do m´etodo de Euler, a partir do ponto (ti, yi), com passo h.
Exerc´ıcio 6.3.7 Considere o problema de condi¸c˜ao inicial
{ y′(t) = ty^2 , y(1) = 2. Determine um valor aproximado para y(1.1), usando o m´etodo de Heun, dado por
k 1 = f (ti, yi); k 2 = f (ti + h, yi + hk 1 );
yi+1 = yi +
h 2 (k 1 + k 2 ),
com h = 0. 05.
Resolu¸c˜ao: Seja f (t, y) = ty^2. Temos que
y(1) = y 0 = 2 y(1.05) ≈ y 1 = y 0 + h 2 (k 1 + k 2 ) = 2 + 0.025(k 1 + k 2 ). Por outro lado k 1 = f (t 0 , y 0 ) = f (1, 2) = 4 k 2 = f (t 0 + h, y 0 + hk 1 ) = f (1. 05 , 2 .2) = 5. 082.
Assim, y(1.05) ≈ y 1 = 2. 22705. Continuando a aplica¸c˜ao do m´etodo
y(1.1) ≈ y 2 = y 1 +
h 2 (k 1 + k 2 ) = 2.22705 + 0.025(k 1 + k 2 ).
Para este segundo passo temos que voltar a calcular k 1 e k 2. Assim, k 1 = f (t 1 , y 1 ) = f (1. 05 , 2 .22705) = 5. 207739 k 2 = f (t 1 + h, y 1 + hk 1 ) = f (1. 1 , 2 .487437) = 6. 806077.
Logo, y(1.1) ≈ y 2 = 2. 5273954.
o m´etodo diz-se consistente com o PCI (6.2). O m´etodo diz-se de ordem p > 0 se existir um C > 0 tal que |Ti+1| ≤ Chp+1, i = 0,... , n − 1 ,
isto ´e, se |Ti+1| = O(hp+1), i = 0,... , n − 1.
Da defini¸c˜ao anterior conclui-se que o erro de truncatura local ´e definido com sendo
Ti+1 = y(ti+1) − y(ti) − hφ(ti, y(ti); h).
Assim, o erro local pode ser determinado atrav´es dos seguintes passos: (i) substituir na express˜ao que define o m´etodo num´erico a solu¸c˜ao aproximada no ponto ti+1, yi+1, pela solu¸c˜ao exacta y(ti+1); (ii) considerar a hip´otese y(ti) = yi; (iii) efectuar o desenvolvimento em s´erie de Taylor de y(ti+1) em torno de ti. Certos autores consideram a defini¸c˜ao de erro local de forma diferente. E muito frequente´ encontrar a defini¸c˜ao de erro de truncatura local no ponto ti+1 como sendo
Ti+1 = y(ti+1) − y(ti) h
− φ(ti, y(ti); h), i = 0,... , n − 1.
De acordo com esta defini¸c˜ao (que ´e a que foi usada em Ferreira e Patr´ıcio (1999)) o m´etodo (6.11) tem ordem p se |Ti+1| = O(hp), i = 0,... , n − 1.
Observa¸c˜ao 6.
- Um m´etodo ´e consistente se tiver, pelo menos, ordem um ou, o que ´e equivalente, se φ(t, y; 0) = f (t, y).
- O erro local para o m´etodo de Taylor de fun¸c˜ao incremento (6.7) ´e dado por
Ti+1 = hk+ (k + 1)!
y(k+1)(ξ), ξ ∈ (ti, ti+1),
ou seja, Ti+1 = O(hk+1) e, como tal, o m´etodo tem ordem k.
- Para o caso particular do m´etodo de Euler temos que Ti+1 = O(h^2 ) e assim sendo o m´etodo tem ordem um.
O estudo da ordem para os m´etodos de Runge-Kutta n˜ao ´e uma tarefa f´acil e, como tal, n˜ao o iremos efectuar neste curso introdut´orio. A t´ıtulo exemplificativo consideremos este estudo apenas para um caso muito simples.
Exemplo 6.11 O m´etodo de Heun ´e dado por
k 1 = f (ti, yi); k 2 = f (ti + h, yi + hk 1 );
yi+1 = yi + h 2 (k 1 + k 2 ).
Vamos determinar qual o seu erro local e, consequentemente, qual a sua ordem. Atendendo `a defini¸c˜ao de erro local temos que, supondo y(ti) = yi, Ti+1 = y(ti+1) − yi+1. Desenvolvendo y(ti+1) em s´erie de Taylor em torno do ponto ti temos
y(ti+1) = yi + hf (ti, yi) + h^2 2
df dt
(ti, yi) + h^3 6
d^2 f dt^2
(ti, yi) + · · ·.
Por outo lado, considerando o desenvolvimento de yi+1, recorrendo `a s´erie de Taylor para duas vari´aveis,
yi+1 = yi+
h 2
( f (ti, yi) + f (ti, yi) + h(ft + fyf )(ti, yi) +
h^2 2 (ftt + 2f fty + f 2 f yy)(ti, yi) + · · ·
) .
Subtraindo membro a membro, e uma vez que y(ti) = yi, temos
Ti+1 =
h^3 12 (ftt + 2f fty + f 2 fyy − 2 ftfy − 2 f f (^) y^2 )(ti, yi) + · · ·.
Assim sai que
Ti+1 = h^3 12
(ftt + 2f fty + f 2 fyy − 2 ftfy − 2 f f (^) y^2 )(ξ, y(ξ)), ξ ∈ (ti, ti+1).
Como Ti+1 = O(h^3 ) conclu´ımos que o m´etodo de Heun tem ordem 2.
Exerc´ıcio 6.3.9 Mostre que o m´etodo de Runge (6.9) tem ordem dois.
Vamos considerar agora a defini¸c˜ao de erro global.
Defini¸c˜ao 6.12 Considere-se o PCI (6.2) e um m´etodo num´erico de passo ´unico expl´ıcito (6.11) que determine aproxima¸c˜oes yi para a solu¸c˜ao exacta y(ti), i = 0, 1 ,... , n. A e(ti) = y(ti) − yi chama-se erro global do m´etodo no ponto ti. Se
lim h→ 0 max 1 ≤i≤n |e(ti)| = 0,
o m´etodo diz-se convergente.
Note-se que, uma vez que o n´umero de vezes que aplicamos um determinado m´etodo itera-
tivo ´e n = b − a h
, podemos ser levados a afirmar que o erro global e(ti) dever´a ser proporcional
a Ti h
. Por outras palavras, se Ti = O(hp+1) tudo nos leva a crer que e(ti) = O(hp). De facto,
o que se pode concluir ´e o que vem expresso no pr´oximo teorema.
Teorema 6.13 Seja y(t) a ´unica solu¸c˜ao do PCI, bem posto, (6.2) e (6.11) um m´etodo num´erico que supomos ser consistente com o problema e ter ordem p, isto ´e, |Ti+1| ≤ Chp+1, i = 0,... , n − 1 , p ≥ 1. Se φ verificar as hip´oteses do Teorema de Picard ent˜ao
|e(ti)| ≤
C
L
hp^
[ eL(ti−a)^ − 1
] , i = 1,... , n.
sendo L a constante de Lipschitz de φ.
- Estabilidade. Para provar que o m´etodo ´e est´avel vamos provar que φ(t, y; h) ´e lips- chitziana, na vari´avel y, em D = {(t, y; h) : a ≤ t ≤ b, y ∈ IR, 0 ≤ h ≤ h 0 }. Seja L a constante de Lipschitz de f (t, y) na vari´avel y. Ent˜ao
|φ(t, y 1 ; h) − φ(t, y 2 ; h)| =
∣∣ ∣∣^1 2
[f (t, y 1 ) + f (t + h, y 1 + hf (t, y 1 ))]
[f (t, y 2 ) + f (t + h, y 2 + hf (t, y 2 ))]
∣∣ ∣∣
[L|y 1 − y 2 | + L|y 1 + hf (t, y 1 ) − y 2 − hf (t, y 2 )|]
≤ L|y 1 − y 2 | +
hL^2 |y 1 − y 2 |
( L +
hL^2
) |y 1 − y 2 |.
Assim φ(t, y; h) satisfaz a condi¸c˜ao de Lipschitz, na vari´avel y, em D sendo a sua constante de Lipschitz dada por
K =
( L +
h 0 L^2
) .
Finalmente, tanto φ como f s˜ao cont´ınuas em D. Est´a assim provada a estabilidade do m´etodo.
Atendendo `a consistˆencia e estabilidade temos que o m´etodo converge para a solu¸c˜ao de (6.2).
6.3.6 M´etodos de passo ´unico impl´ıcitos
Nesta sec¸c˜ao vamos considerar a classe dos m´etodos de passo ´unico impl´ıcitos da forma (6.6). N˜ao havendo possibilidade de explicitar o valor de yi+1 temos necessidade de o calcular re- solvendo a equa¸c˜ao (geralmente n˜ao linear)
yi+1 − yi − hφ(ti, ti+1, yi, yi+1; h) = 0.
Usualmente considera-se um m´etodo num´erico na resolu¸c˜ao desta equa¸c˜ao. Se considerarmos o m´etodo de Newton, a primeira quest˜ao a resolver ´e a da determina¸c˜ao de uma aproxima¸c˜ao inicial y(0) i+1. Normalmente toma-se para aproxima¸c˜ao inicial o valor de yi; outra hip´otese ser´a a de considerar a aproxima¸c˜ao inicial obtida pela aplica¸c˜ao de um m´etodo expl´ıcito. Deteminado o valor de y i(0)+1 temos que
y( ik+1+1) = y( ik+1) −
F (y i(+1k)) F (y i(+1k))
, k = 0, 1 ,... ,
sendo F (yi+1) = yi+1 − hφ(ti, ti+1, yi, yi+1; h).
Os m´etodos impl´ıcitos s˜ao usados visto que, em geral, s˜ao mais precisos e menos sens´ıveis a erros que os m´etodos expl´ıcitos. Por outro lado, o esfor¸co computacional exigido no c´alculo de yi+1 ´e, para os m´etodos impl´ıcitos, muito maior. Assim, estes m´etodos s´o devem ser usados quando h´a necessidade de uma precis˜ao muito elevada em problemas sens´ıveis a erros. Exemplos comuns de m´etodos impl´ıcitos s˜ao o chamado m´etodo de Euler impl´ıcito, dado pela express˜ao
yi+1 = yi + hf (ti+1, yi+1), i = 0,... , n − 1 , y(a) = α,
e o m´etodo dos trap´ezios, dado por
yi+1 = yi + h 2
[f (ti, yi) + f (ti+1, yi+1), i = 0,... , n − 1 , y(a) = α.
Apesar do estudo da consistˆencia e convergˆencia de um m´etodo iterativo ter sido efectuado apenas para m´etodos expl´ıcitos, estes conceitos ainda s˜ao v´alidos para m´etodos impl´ıcitos. Para o m´etodo impl´ıcito (6.6) o erro de truncatura local ´e definido por
Ti+1 = y(ti+1) − y(ti) − hφ(ti, ti+1, y(ti), y(ti+1); h), i = 0,... , n − 1.
Vamos considerar o seguinte exerc´ıcio.
Exerc´ıcio 6.3.11 Considere o m´etodo dos trap´ezios na resolu¸c˜ao de um problema de condi¸c˜ao inicial.
- Determine a ordem e o erro de truncatura local do m´etodo.
- Aplique o m´etodo ao problema de condi¸c˜ao inicial
{ y′(t) = −ty^2 , y(0) = 2 , e obtenha uma aproxima¸c˜ao em t = 1 usando h < 1. (Considere a solu¸c˜ao exacta positiva em [0, 1].)
Resolu¸c˜ao: 1. Atendendo `a defini¸c˜ao de erro local temos que Ti+1 = y(ti+1) − yi+1, onde y(ti) = yi. Desenvolvendo y(ti+1) e yi+1 em s´erie de Taylor em torno do ponto ti, temos que
y(ti+1) = yi + hf (ti, yi) +
h^2 2
df dt (ti, yi) +
h^3 6
d^2 f dt^2 (ti, yi) + · · ·
e yi+1 = yi + h 2
( f (ti, yi) + f (ti, yi) + h df dt
(ti, yi) + h^2 2
d^2 f dt^2
(ti, yi) + · · ·
) .
Subtraindo membro a membro vem
Ti+1 = −
h^3 12
d^2 f dt^2 (ti, yi) + · · ·.
Assim sai que Ti+1 = −
h^3 12 y′′′(ξ), ξ ∈ (ti, ti+1).
Como Ti+1 = O(h^3 ) temos que o m´etodo dos trap´ezios tem ordem 2.
6.4 Sistemas de equa¸c˜oes diferenciais
A teoria apresentada nas sec¸c˜oes precedentes pode ser facilmente generalizadas para sistemas de equa¸c˜oes diferenciais ordin´arias de primeira ordem. Todos os m´etodos num´ericos apresen- tados podem ser adaptados ao c´alculo da solu¸c˜ao aproximada do PCI { Y ′(t) = F (t, Y ), t ∈ [a, b], Y (a) = α. (6.12)
onde
Y (t) =
Y 1 (t) Y 2 (t) .. . YN (t)
, F (x) =
F 1 (t, Y ) F 2 (t, Y ) .. . FN (t, Y )
Os m´etodos num´ericos ir˜ao, neste caso, determinar aproxima¸c˜oes Y (i)^ para Y (ti). O m´etodo de Euler, por exemplo, ´e dado por
Y (i+1)^ = Y (i)^ + hF (ti, Y (i)).
Equa¸c˜oes diferenciais de ordem superior a um. Uma situa¸c˜ao importante onde surgem sistemas de equa¸c˜oes diferenciais ´e quando pretendemos resolver uma equa¸c˜ao diferencial de ordem superior a um. Note-se que qualquer equa¸c˜ao diferencial de ordem N pode ser escrita como um sistema de N equa¸c˜oes diferenciais de primeira ordem. A forma como essa passagem se processa ´e bastante simples e pode ser facilmente compreendida com a ajuda de um exemplo.
Exemplo 6.15 Consideremos o problema de condi¸c˜ao inicial
y′′^ − 3 y′^ + 2y = 0, y(0) = y′(0) = 1.
Efectuando a mudan¸ca de vari´avel z = y′^ obtemos o problema de condi¸c˜ao inicial de primeira ordem
y′(t) = z z′(t) = 3 z − 2 y y(0) = 1 z(0) = 1
[ y z
]′ (t) =
[ z 3 z − 2 y
] ,
[ y z
] (0) =
[ 1 1
] .
Exerc´ıcio 6.4.1 Converta num sistema de equa¸c˜oes diferenciais de primeira ordem o problema
y′′′^ − 0 .1(1 − y^2 )y′^ + y = 0, y(0) = 1, y′(0) = y′′(0) = 0.
Resolu¸c˜ao: Efectuando a mudan¸ca de vari´avel z = y′^ e w′^ = y′′^ obtemos o problema de condi¸c˜ao inicial de primeira ordem
y′(t) = z z′(t) = w w′(t) = 0 .1(1 − y^2 )z − y y(0) = 1 z(0) = 0 w(0) = 0
Exerc´ıcio 6.4.2 Considere a equa¸c˜ao diferencial y′′^ + 4ty′^ + 2y^2 = 0 com condi¸c˜oes iniciais y(0) = 1 e y′(0) = 0. Com h = 0. 1 , utilize o m´etodo de Euler e de Heun para obter aproxima¸c˜oes para y(0.2) e y′(0.2).
Resolu¸c˜ao: Seja z = y′. Assim o nosso problema ´e equivalente a
y′(t) = z z′(t) = − 4 tz − 2 y^2 y(0) = 1 z(0) = 0
[ y z
]′ (t) =
[ z − 4 tz − 2 y^2
] ,
[ y z
] (0) =
[ 1 0
] .
Seja F
( t,
[ y z
])
[ z − 4 tz − 2 y^2
] .
- M´etodo de Euler [ y z
] (0) =
[ y z
](0)
[ 1 0
] ;
[ y z
] (0.1) ≈
[ y z
](1)
[ y z
](0)
t 0 ,
[ y z
](0) (^) =
[ 1 − 0. 2
] ;
[ y z
] (0.2) ≈
[ y z
](2)
[ y z
](1)
t 1 ,
[ y z
](1) (^) =
[
- 98 − 0. 392
] .
Temos assim que y(0.2) ≈ 0. 98 e y′(0.2) ≈ − 0. 392.
- M´etodo de Heun [ y z
](0)
[ y z
] (0) =
[ 1 0
] ,
[ y z
] (0.1) ≈
[ y z
](1)
[ y z
](0)
h 2
[K 1 + K 2 ],
onde K 1 = F
t 0 ,
[ y z
](0) (^) =
[ 0 − 2
]
K 2 = F
t 0 + h,
[ y z
](0)
(^) =
[ − 0. 2 − 1. 92
] .
Logo (^) [ y z
] (0.1) ≈
[ y z
](1)
[
- 99 − 0. 196
] .
