Métodos
Numéricos para
Equações Diferenciais
Ordinárias
Começamos por relembrar o caso mais simples
de uma equação diferencial ordinária,
y'(t)=f(t)
cuja solução é dada simplesmente
pelo integral
| y(t) = |
ó
õ |
t
0 |
f(s)ds + C |
onde C é uma constante a determinar. Sabendo
o valor inicial y(0)=y0 obtemos imediatamente C=y0
,
e
trata-se de uma fórmula explícita para a solução
y (desde que f seja integrável).
De forma semelhante, quando temos um caso mais
geral em que
y'(t)=f( t, y(t) ),
definimos um problema com valores iniciais, denominado
problema
de Cauchy:
ì
í
î |
y'(t) = f( t, y(t) )
y(t0) = y0 |
cuja solução deverá verificar
| y(t) = y0 + |
ó
õ |
t
t0 |
f( s, y(s) )ds |
o que se trata de uma equação implícita
em y.
Para certas funções f particulares,
como a função linear
f(x,y) = a(x) y + b(x),
é possível obter uma fórmula
explícita
| y(t) = y0 eA(t)
+ eA(t) |
ó
õ |
t
t0 |
b(s) e-A(s) ds , |
em que A(t) = òt0t
a(s)
ds, e assim a solução
y fica condicionada ao cálculo
dos integrais.
Para um f geral, a fórmula explícita
já não é possível, e apenas podemos considerar
métodos numéricos. Mesmo no caso em que temos fórmula
explícita, o cálculo dos integrais poderá também
apenas ser efectuado numericamente.
Antes de estudar os métodos numéricos,
convém relembrar em que condições podemos obter uma
solução única.
Teorema (Picard-Lindelöf)
Seja D um aberto de R2,
e
f: D Ì R2®
R
uma
função que verifica as condições:
(i) f é contínua em D
(ii) f é Lipschitziana na segunda
variável, ou seja, existe L>0 :
|f(t,y)-f(t,w)| <
L |y-w| , para quaisquer (t,y), (t,w)
pertencentes
a D.
então, dado qualquer (t0 ,
y0) em D, existe uma
única solução do problema de Cauchy definida numa
vizinhança de t0 .
demonstração: Sai fora do âmbito do curso,
no entanto podemos referir que pode ser demonstrado pelo teorema do ponto
fixo aplicado a um espaço funcional conveniente (em dimensão
infinita).
Observação:
A condição (ii) pode ser facilmente verificada se
mostrarmos que,
num conjunto limitado D, a derivada ¶f
/ ¶y existe e é contínua.
Exemplo:
Para além do caso linear, em que apresentámos mesmo a
solução explícita (assumimos a e b funções
contínuas), damos um outro exemplo em que o facto da solução
ser única depende das condições iniciais consideradas.
Trata-se da equação diferencial y'(t) = 2Öy(t)
. Esta equação diferencial apresenta várias funções
C1
| yp(t) = |
ì
í
î |
(t - p)2, se t > p
0 , se t < p |
que são soluções do problema
de Cauchy com y(1) = 0, bastando considerar p >1.
Isto não contraria o teorema de Picard-Lindelöf, pois notamos
que a condição (ii) não é verificada em nenhum
ponto (x,0),
porque ¶f/¶y
= -1/Öy, ou seja, há
uma singularidade quando y=0.
Se alterarmos a condição inicial, escolhendo y(1)
= 4, então garantimos solução única numa
vizinhança de t=1.
Métodos de Taylor
Consideramos um primeiro esquema numérico, muito simples e que
consiste em efectuar uma truncatura da expansão da série
de Taylor e substituir a derivada pela expressão explícita
dada na equação diferencial ordinária. As derivadas
seguintes são obtidas efectuando a derivação sucessiva
do segundo membro.
Método de Euler
Assim, definindo um valor h>0 e pontos igualmente espaçados
tm=
t0+ mh, obtemos
y(tm+1)
= y(tm+h) = y( tm ) + y'( tm )h + y''(
xm
)h2/2,
com xm em
]tm
,tm+1[,
pelo resto de Lagrange da série de Taylor.
A partir deste desenvolvimento, podemos substituir
y'(
tm ) = f( tm, y(tm) )
e desprezar o erro de truncatura y''( xm
)h2/2.
Portanto, da expressão truncada, y(tm+1)
= y( tm ) + f( tm, y(tm) ) h,
definindo ym= y(tm),
e partindo do valor y0= y(t0), obtemos o
| Método de Euler :
ym+1 = ym + f( tm ,
ym ) h , |
que nos irá fornecer aproximações para os valores
y( tm+1 )
a partir dos valores ym .
Como os valores ym são
também aproximações, os erros vão acumular
em cada iteração.
Erro do método de Euler
Usando o erro de truncatura, reparamos que
y(tm+1)
- ym+1 = ( y( tm
) + y'( tm )h + y''(
xm
)h2/2 ) - ( ym + f( tm
,
ym ) h)
= ( y( tm ) - ym)
+ (f( tm , y(tm ) ) - f( tm ,
ym )) h + y''(
xm
)h2/2 ,
definindo em = y( tm ) -
ym , e como f é lipschitziana,
|f( tm , y(tm ) ) -
f( tm , ym )| < L|y(tm ) - ym|
= L |em|
obtemos, supondo que |y''(x)| <B
em [t0 ,tm+1]
|em+1|
<
|em| + L|em| h + |y''(
xm
)|h2 / 2 < (1+ Lh)|em|
+ B h2 / 2,
o que permite obter recursivamente, com K=1+
Lh,
|em|
<
Km|e0| + (1+K+...+Km-1)
B h2 / 2, ou seja,
| |em|
<
Km|e0| + |
Km-1
K-1 |
B h2 / 2 . |
Notando que eLh > Lh+1
= K, obtemos também
| |em|
<
emLh|e0| + |
emLh-1
Lh |
B h2 / 2 , |
ou seja, como mh=tm-t0
| |em|
<
exp((tm-t0)L)|e0| + |
exp((tm-t0)L) -1
2 L |
B h , |
o que torna claro que quando e0=0
(erro inicial nulo), temos em= O(h), e põe
em evidência a forma como o erro tende para zero quando h
converge para zero.
Observação:
Pela estimativa de erro, podemos também
verificar que quando o erro inicial não é nulo, esse erro
pode propagar-se de maneira significativa já que aparece multiplicado
por uma exponencial. Pode assim ocorrer um problema de estabilidade, para
valores tm-t0 grandes, e que pode ocorrer
mesmo com simples erros de arredondamento.
Com efeito, se considerarmos simplesmente y'(x) = c y(x), cuja solução
é uma exponencial, teremos
em+1
= em + c h em +
y''(
xm )h2/2 ,
e portanto obtemos um termo (1 + c h)m e0,
que irá tender para infinito se e0 não
for nulo e se c>0. Mesmo no caso em que c<0, se h>-2/c, teremos
o módulo desse termo a tender para infinito, o significa que apenas
haverá estabilidade para certos valores de h (estabilidade
condicionada).
Método de Euler implícito
É assim conveniente encontrar outros métodos para os
quais não se verifique o problema de instabilidade. Uma solução
possível são os métodos implícitos, por exemplo,
| Método de Euler implícito
: ym+1 = ym
+ f( tm+1 , ym+1
) h , |
em que o valor de ym+1
é calculado por resolução de uma equação
não linear.
Com este método, a mesma equação y'(x) = c y(x),
leva agora à estimativa de erro
em+1
= em + c h em+1+
y''(
xm )h2/2
,
logo,
em+1
= em /(1- c h) + y''(
xm
)h2/2/(1- c h) ,
e assim obtemos o termo (1 - c h)-m e0,
que apenas tenderá para infinito se |1 - c h| < 1, e isto
apenas pode acontecer se hc>2, pelo que basta considerar um h
suficiente pequeno para o evitar.
Outros métodos de Taylor
O método de Euler é o caso mais
simples de uma aproximação feita pela expansão em
série de Taylor.
Efectuando a derivada total em ordem a t
de f( t, y(t) ) = y'(t) podemos obter a segunda derivada de y.
Por exemplo, sendo y'(t) = f( t, y(t) ) =
t2y(t), reparamos que
y''(t) = d/dt (f( t, y(t)
)) = d/dt (t2y(t)) = 2t y(t)+ t2y'(t)
=(2t + t4)y(t)
e portanto neste caso considerando a fórmula
de Taylor
y(tm+1)
= y( tm ) + y'( tm )h + y''( tm
)h2/2 + y'''(
xm
)h3/6,
com xm em
]tm
,tm+1[,
podemos considerar o método de Taylor de segunda ordem,
ym+1
= ym + tm2 ym h
+ (2 tm + tm4)ym h2/2
,
com o erro de truncatura y'''(
xm
)h3/6.
Calculando as derivadas seguintes de y
e efectuando uma expansão maior da série de Taylor, podemos
obter um método melhor, com um erro de truncatura de ordem superior.
A desvantagem destes métodos é
exigir o sucessivo cálculo simbólico das derivadas, o que
pode conduzir a grandes expressões, e consequentemente a um maior
número de cálculos.
Generalização
do método de Euler a sistemas de EDO's
Da mesma forma que considerámos o problema de Cauchy definido
para uma função com valores reais, podemos considerar uma
incógnita y que é uma função com valores
vectoriais, tendo-se
ì
í
î |
y'(t) = F( t, y(t)
)
y(t0) = y(0) |
em que y: R ®
Rn , F : R´Rn
®
Rn, com o vector inicial y(0) ÎRn
.
Tal como explicitado para o caso de uma única
dimensão, no caso de várias dimensões o método
de Euler é semelhante, ficando
| y(m+1)
= y(m) + F( tm ,
y(m)
) h . |
Exemplo:
Como exemplo apresentamos o caso de uma equação diferencial
de ordem 2,
u''(t) + u'(t)2 u(t) = t2u(t)2,
com condições iniciais u(2) = 3, u'(2) = -4.
Definindo y1(t) = u(t), y2(t) = u'(t),
obtemos as 2 equações de ordem 1,
y1'(t) = y2(t),
y2'(t) = t2 y1(t)2 -
y2(t)2 y1(t),
podendo definir o problema de Cauchy vectorial,
ì
í
î |
(y1'(t), y2'(t))
= ( y2(t) , t2
y1(t)2 - y2(t)2
y1(t)
)
(y1(2), y2(2))
= ( 3, -4 ) |
e desta forma, o método de Euler corresponde
a efectuar as iterações nas duas componentes,
ì
ï
î |
y1(m+1) = y1(m) + y2(m)h
y2(m+1) = y2(m) +
(
tm2
(y1(m))2
-
(y2(m))2
y1(m) )h |
Assim, como y(0) = (3
,
-4), a primeira iterada será, neste
caso,
y1(1) = 3 +(-4)h,
y2(1) = -4 +(22´32-
(-4)2´3
) h = -4 - 12h,
ou seja, y(1) = (3 - 4h, -4
-12h).
O valor da 1ª componente, 3 - 4h, é uma aproximação
de u(2+h) = y1(t1) ,
e o valor da 2ª componente, -4 - 12h, é uma
aproximação de u'(2+h) = y2(t1).
Generalidades sobre os métodos
numéricos
Consideramos a equação diferencial y'(t)
= f( t, y(t) )
aproximada por um esquema do tipo ym+1
= ym + f( tm ,
ym ; h ) h
(no caso do método de Euler, f(
tm , ym ; h )= f( tm , ym )
)
Há três aspectos a ter em conta quando analisamos um esquema
numérico
(i) Consistência, (ii) Convergência, (iii)
Estabilidade.
(i) Consistência.
A consistência diz respeito à aproximação
efectuada pela discretização local da equação
diferencial.
Assim, definimos o erro de discretização local
E(h) = y'(tm) - (ym+1
- ym)/h = f( tm , y(tm) ) - f(
tm , ym ; h )
e dizemos que o esquema é consistente se E(h) ®
0, quando h ® 0,
e que tem ordem de consistência p se E(h) =
O(h p).
No caso do método de Euler, como E(h) = y''(
xm
)h /2 é um O(h) quando y'' é limitada,
concluímos que se trata de um esquema
com ordem de consistência 1.
(ii) Convergência.
A convergência é verificada quando temos
e(h) = maxm=0,...,M
|em| = maxm=0,...,M |y'(tm)
- ym | ® 0
e dizemos que tem convergência de ordem
p
se e(h) = O(h p).
(iii) Estabilidade.
A estabilidade está relacionada com a propagação
de erros de arredondamento, nos passos intermédios, e pode ser grave
no caso de métodos multipasso. Por vezes também é
utilizada a designação de estabilidade para o próprio
condicionamento, ou seja, no que diz respeito à influência
face aos erros iniciais.
Teorema: Um esquema estável que tiver ordem de consistência
p
então também terá ordem de convergência
p.
Demonstração: exercício... repetir os passos
da demonstração do erro do método de Euler.
Observação: Este teorema mostra que, desprezando
o erro inicial, basta verificar a ordem de consistência do método
para assegurar a mesma ordem de convergência. O recíproco
é também válido. No caso de métodos unipasso
a estabilidade é verificada mediante a hipótese da função
ser Lipschitziana.
Métodos de Runge-Kutta
Consideramos uma aproximação baseada nos métodos
de Taylor, mas em que não é necessário calcular explicitamente
as derivadas de f, sendo subsituída pelo desenvolvimento
em série de Taylor, com duas variáveis.
Por exemplo, consideramos
y(t+h) = y( t ) + y'( t )h + y''(t)h2/2
+ O(h3),
e como y'( t ) = f( t, y( t )), consideramos
o desenvolvimento de f, nas 2 variáveis,
f(x+a,
y+b) = f(x, y) + a¶1
f(x,
y) + b¶2
f(x,
y) + O(ab).
Como temos,
y'(t)h+ y''(t)h2/2
= f(t, y)h + d/dt (f( t, y(t) ))h2/2 =
= f(t, y)h + (¶1
f(t,
y) + f(t,y) ¶2
f(t,
y) )h2/2
podemos multiplicar a 1ª igualdade por h, e, subtraindo, obtemos
f(x+a,
y+b)h - (y'(t)h+ y''(t)c2/2)
=
= a¶1
f(x,
y)h + b¶2
f(x,
y)h - (¶1
f(t,
y) + f(t,y) ¶2
f(t,
y) )h2/2 + O(ab)h.
Podemos escolher a,b
convenientemente, de forma a ficarmos apenas com O(ab)h
e reparando que se a,b
= O(h), então O(ab)
= O(h3)...
é esta a ideia subjacente aos
métodos
de Runge-Kutta (introduzidos por C. Runge no final do século
XIX).
-- Escolhendo a
= h/2, b
= f(t,y)h/2, ficamos com
y(t+h) = y(t) + f(x+h/2
,
y+f(t,y)h/2) h + O(h3),
o que nos leva a um método de segunda ordem
Método do ponto-médio (RK ordem 2):
ym+1
= ym + f(tm+h/2
,
ym+fmh/2)
h ,
com fm
= f(tm,ym) |
Podemos ainda considerar outro tipo de métodos, deduzindo uma
fórmula mais geral,
já que, de forma semelhante, podemos deduzir,
Af(x, y)h+Bf(x+a,
y+b)h = (A+B)f(x, y)h + B (a¶1
f(x,
y)h + b¶2
f(x,
y))h + O(ab)h
y'(t)h+ y''(t)h2/2
= f(t, y)h + (¶1
f(t,
y) + f(t,y) ¶2
f(t,
y) )h2/2
ficando com as equações, A+B =1 , aB=h/2,
bB=f(t,y)h/2
-- Escolhendo A=B=1/2, a=h,
b=hf(t,y)
Método de Euler modificado (RK ordem 2):
ym+1
= ym + fmh/2+
f(tm+h,
ym+fmh)
h/2 , |
-- Escolhendo A=1/4, B=3/4, a=2h/3,
b=
2hf(t,y)/3
temos
| ym+1
= ym + fmh/4
+
3f(tm+2h/3,
ym+2fmh/3)
h/4 . |
Métodos de Runge-Kutta
de ordem superior
Por um processo semelhante, mas com cálculos exaustivos, é
possível obter outros métodos RK de ordem superior.
Enunciamos apenas o mais popular, o método de Runge-Kutta de
ordem 4,
Método de Runge-Kutta de ordem 4:
ym+1
= ym + (A1+2A2+2A3+A4)
/6 ,
com
A1= fmh,
A2= f(tm+h/2,
ym+ A1/2)
h ,
A3= f(tm+h/2,
ym+ A2
/2)
h, A4=
f(tm+h,
ym+ A3)
h |
Observação: Tal como no
caso do método de Euler, estes métodos podem ser imediatamente
adaptados para sistemas de equações diferenciais ordinárias,
usando a notação vectorial.
Métodos implícitos
e preditor-corrector
Relembramos o método de Euler implícito,
ym+1
= ym + f( tm+1 , ym+1
) h
cujo valor de ym+1deve
ser encontrado resolvendo a equação não linear em
z,
z = ym +
f( tm+1, z )
h.
Para esse efeito, podemos definir o processo iterativo do ponto fixo
zk+1
= ym + f( tm+1,
zk ) h = g( zk )
que converge para h pequeno, porque se
|¶2
f(t,
y)| < L, o que acontece normalmente
para funções regulares (e o que
implica que é Lipschitziana), porque
|g'(x)| = |¶2
f(t,
x)h| < Lh < 1, desde que h seja suficientemente
pequeno.
Assim, a função g é
contractiva e temos convergência local do método do ponto
fixo.
Note-se que se começarmos com z0
=ym, obtemos z1
= ym + f( tm+1,
ym ) h,
o que está próximo do valor ym
+ fm h, e
que pode ser visto como a aproximação
de ym+1 usando
o método de Euler explícito.
Método dos trapézios
Uma outra possibilidade é usar a relação
| y(tm+1)
= y(tm) + |
ó
õ |
tm+1
tm |
f( s, y(s) )ds |
aproximando o integral com a regra dos trapézios,
ficando com o método
| ym+1
= ym + (f(
tm , ym)+f( tm+1,
ym+1))h/2 |
que também é um método implícito.
Notamos, no entanto, que é habitual,
usar nos métodos implícitos a aproximação inicial
obtida
por um explícito, e depois aplicar
o método do ponto fixo, designando-se esta técnica por preditor-corrector.
Assim, o valor ym+1
pode ser aproximado pelo método de Euler explícito, ficando
a primeira iterada
| ym+1
= ym + (f(
tm , ym)+f( tm+1,
ym+hfm)) h/2 |
que é exactamente o método de Euler
modificado, um método RK de 2ª ordem.
Métodos multipasso
Reparamos que até aqui, procurámos apenas usar a informação
da última iterada, designando-se estes métodos por unipasso.
Mas é claro que uma melhor aproximação poderá
ser obtida considerando outros valores já calculados, tratam-se
dos métodos multipasso.
Os mais conhecidos são os métodos de Adams.
Métodos de Adams-Bashforth
Consideramos ainda a fórmula
| y(tm+1)
= y(tm) + |
ó
õ |
tm+1
tm |
f( s, y(s) )ds |
aproximando o integral através de interpolação
da função f nos pontos já calculados.
Por exemplo, podemos considerar apenas tm-1
e tm .
Se definirmos uma fórmula de quadratura
Q( f ) = A0 f(tm-
h)+A1 f(tm)
que seja correcta para polinómios de grau
1,
h = I( 1 ) = Q(1) = A0 +A1
h tm+ h2/2 = (tm+h)2/2
- tm2/2 = I( t ) = Q( t ) =
( A0 +A1) tm-A0 h
=
h
tm-A0 h
logo tm+ h/2 = tm-A0
, o que dá A0 = - h/2, A1
=
3h/2, e assim
Q( f ) = (-f(tm-1)+3f(tm))h
/
2, pelo que obtemos
| ym+1
= ym + ( - f(
tm-1 ,
ym-1)+3f(
tm , ym)) h/2 |
que corresponde à fórmula de Adams-Bashforth
de 2ª ordem, em que o
erro de truncatura local é 5y'''(c)h3/12.
Nota: Para obter esta estimativa de erro deverá
considerar-se um processo semelhante ao utilizado para demonstrar o erro
de integração, usando o teorema do valor intermédio
para integrais.
Usando mais pontos anteriores, tm-2
, tm-3 , etc. obtemos fórmulas de maior grau.
Métodos de Adams-Moulton
A única diferença face aos anteriores, é considerar
também tm+1
como nó de integração.
É claro que nesse caso tratar-se-à de uma fórmula
implícita.
Por exemplo, efectuando cálculos semelhantes, iremos obter como
primeira aproximação a
regra dos trapézios implícita, que já vimos e
sabemos ter erro -y'''(c)h3/12
Métodos predictor-corrector de Adams
Basta considerar uma fórmula implícita com aproximação
inicial de ym+1
dada pela fórmula explícita.
Assim, por exemplo,
| ym+1
= ym + (fm
+f(
tm+1
, ym
+ (3fm
- fm-1)
h/2)) |
Exemplo
Consideramos o caso da equação diferencial ordinária
de 2ª ordem
y''(t)=y(t), com condições iniciais y(0)=0,
y'(0)=1,
cuja solução é y(t)=sin(t).
Na figura seguinte mostramos a diferença entre a aproximação
obtida
com o método de Euler (de 1ªordem), o método de
Euler modificado
(método de Runge-Kutta de ordem 2) e com um método de
Runge-Kutta
de ordem 4, já apresentado. Considerámos 20 pontos no
intervalo [0, 9],
e podemos constatar que a aproximação obtida por Euler
não é adequada, sendo
necessária um passo h muito mais pequeno para que a aproximação
seja razoável.
No caso do método RK de ordem 2, a aproximação
é melhor, mas apenas com
o método RK de ordem 4 obtemos uma aproximação
excelente.
