Equações Diferenciais
Métodos Adaptativos e Rigidez
Computação – 2º Semestre 2011/2012
Métodos Adaptativos

Por vezes as soluções das EDOs têm
diferentes escalas temporais:



Os métodos com o tamanho do
passo constante teriam que aplicar
sempre um pequeno passo:


Em alguns intervalos de tempo a
solução varia gradualmente;
Noutros as variações são abrutas.
Ineficiente em regiões de variação
gradual.
Os métodos adaptativos podem
variar o tamanho do passo de acordo
com as variações da solução.
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
2
Métodos Adaptativos

A implementação requer uma estimativa do erro local de
truncatura cometido em cada passo:


Com base nesta estimativa o tamanho do passo pode ser
diminuído ou aumentado.
Existem 2 abordagens alternativas para incorporar o
controlo adaptativo do tamanho do passo:

Divisão do passo – usar o método de passo simples duas vezes:



Uma vez com o tamanho do passo completo h
Outra vez com 2 meios passos de tamanho h/2
Estimar o erro local de truncatura cometido no passo n como a
diferença entre os resultados obtidos:
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
3
Método RK Adaptativo de 4ª Ordem

Um passo completo h:

Dois meios passos h/2:

Estimativa do erro local de truncatura:

Correcção de 5ª ordem:
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
4
Métodos Adaptativos

A implementação requer uma estimativa do erro local de
truncatura cometido em cada passo:


Com base nesta estimativa o tamanho do passo pode ser
diminuído ou aumentado.
Existem 2 abordagens alternativas para incorporar o
controlo adaptativo do tamanho do passo:

Métodos RK Embebidos (RK Fehlberg) – calcular no mesmo passo
duas estimativas com métodos de RK de ordens diferentes:


Estimar o erro local de truncatura como a diferença entre os resultados
obtidos
Estes métodos costumam ser mais eficientes que os da divisão
do tamanho do passo.
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
5
Método RK Embebido de 2ª/3ª Ordem

Usa os métodos RK de 2ª e 3ª ordem em simultâneo
para resolver a EDO e controlar o tamanho do passo de
acordo com a estimativa do erro:

Estimativa do erro local de truncatura:
(apenas avalia a função 3 vezes: k1 é o mesmo que o k4 do passo anterior)
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
6
Método RK Embebido de 4ª/5ª Ordem

Usa os métodos RK de 4ª e 5ª ordem em simultâneo:

4ª:

5ª:
(apenas avalia a função 6 vezes: k1,…, k6)
( 5)
( 4)
E

y

y
 Estimativa do erro local de truncatura:
i 1
i 1
i 1
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
7
Métodos RK em MATLAB



2º/3º Ordem: Implementado pela função ode23
4º/5º Ordem: Implementado pela função ode45
Após cada passo, é verificado se o erro está dentro da
tolerância desejada:


Se estiver, yi+1é aceite;
Caso contrário, o cálculo é repetido com o tamanho do passo
reduzido até ser obtida uma estimativa de erro aceitável
(tolerância relativa: por omissão103)
5 Junho 2012
(tolerância absoluta: por omissão106)
Equações Diferenciais – Métodos Adaptativos e Rigidez
8
Métodos RK em MATLAB


[t,y] = ode23(odefun,tspan,y0)
[t,y] = ode45(odefun,tspan,y0)
 y: solução, em que cada coluna representa uma das variáveis e
as linhas os seus valores nos tempos definidos no vector t
 odefun: função que calcula as derivadas do sistema de EDOs
 tspan: intervalo de tempo para resolver o sistema de EDOs:
 Se tspan=[ti tf], os resultados são mostrados para esses


tempos assim como para os tempos intermédios adicionais utilizados.
Se tspan=[t0 t1 … tn], os resultados são mostrados apenas para
os tempos especificados.
y0: vector de valores iniciais para cada uma das variáveis.
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
9
Métodos RK em MATLAB
dy1
1.2y1  0.6y1 y2
dt
dy2
 0.8y2  0.3y1 y2
dt

Exemplo:

com y1(0)=2 e y2(0)=1 para t entre 0 e 20
Solução:
>>
>>
>>
>>
function yp = predprey(t, y)
 yp = [1.2*y(1)-0.6*y(1)*y(2); -0.8*y(2)+0.3*y(1)*y(2)];
tspan = [0 20];
y0 = [2, 1];
[t, y] = ode45(@predprey, tspan, y0);
plot(t,y);
plot(y(:,1),y(:,2));
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
10
Métodos RK em MATLAB
Podem ser definidos parâmetros do algoritmo (options) e
passados parâmetros para a função odefun
[t,y]=ode23(odefun,tspan,y0,options,p1,p2,…)
[t,y]=ode45(odefun,tspan,y0,options,p1,p2,…)
 p1,p2,…: são os parâmetros da função odefun
 options: são os parâmetros do algoritmo


Para definir os parâmetros do algoritmo:
options=odeset('par1',val1,'par2',val2,…)
 'RelTol': tolerância relativa.
 'AbsTol': tolerância absoluta.
 'InitialStep': passo inicial.
 'MaxStep': passo máximo (1/10 do intervalo inicial)
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
11
Métodos RK em MATLAB

Exemplo:

com y (0)=0.5 para t entre 0 e 4
Solução:
function yp = dydt(t, y)
yp = 10*exp(-(t-2)*(t-2)/(2*.075^2))-0.6*y;
>> options=odeset('RelTol',1e-4);
>> ode23(@dydt, [0 4], 0.5);
>> ode23(@dydt, [0 4], 0.5, options);
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
12
Métodos Multipasso



Os métodos de passo simples usam a informação de um único
ponto ti para calcular uma estimativa da variável dependente yi+1
num ponto futuro ti+1.
Os métodos multipasso baseiam-se na utilização da informação
calculada em pontos anteriores.
A função MATLAB ode113 implementa um método multipasso
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
13
Método de Heun

Melhoria do método de Euler baseada na média das derivadas calculadas no
início e na estimativa do fim do intervalo:

Primeiro estima o novo valor de y, depois corrige-o baseado na derivada
calculada nesse novo valor.

predictor O(h2):

corrector O(h3):
29 Maio 2012
Equações Diferenciais – Problemas de Valor Inicial
14
Método Multipasso Heun

Melhorar o preditor para O(h3) usando a informação de yi1:
(a) Predictor

0
i 1
y
(b) Corrector
Versão iterativa:
y
m
i 1
 f ti , y 2h
5 Junho 2012
m
i
j
i 1
y
f ti , yim   f ti 1 , yij11 
y 
h
2
m
i
Equações Diferenciais – Métodos Adaptativos e Rigidez
15
Rigidez


Um sistema rígido (stiff) contém componentes que variam muito
rapidamente e outras que variam muito lentamente.
Em alguns casos a variação rápida ocorre numa pequena fracção
dy
de tempo. Ex:
t
dt
 1000y  3000  2000e
cuja solução com y(0)=0 é:
y  3 0.998e1000t  2.002et


5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
16
Rigidez




Estabilidade: amplificação ou decréscimo dos erros cometidos.
Um método numérico é estável se o erro cometido numa fase
do algoritmo não é amplificado para as fases seguintes.
Uma equação diferencial diz-se mal condicionada se os erros
numéricos são amplificados independentemente do método
Um sistema rígido (stiff) requer, em alguns métodos numéricos,
passos muito pequenos para evitar a instabilidade numérica.
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
17
Rigidez

Exemplo:

Com um pequeno erro no valor inicial:

logo:
O erro cresce exponencialmente: problema instável
O erro permanece constante: problema neutralmente estável
O erro decresce exponencialmente: problema estável
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
18
Rigidez

Exemplo:

Método de Euler (explícito):

Critério de estabilidade:
(factor de amplificação)


Região estável:
Se a=1000, para manter a estabilidade é necessário h2/1000
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
19
Rigidez

Exemplo:

Método de Euler Implícito:

Critério de estabilidade:
(incondicionalmente estável!)
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
20

Rigidez

Exemplo:
dy
 1000y  3000  2000et
dt

cuja solução com y(0)=0 é:
y  3 0.998e1000t  2.002et

Solução:
a)
b)
Método de Euler (explícito)
Método de Euler Implícito
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
21
Rigidez

Sistemas de ODEs rígidos também podem ser resolvidos por
métodos implícitos.
Exemplo:

Cuja solução para y1(0)=52.29 e y2(0)=83.82 é:

Método de Euler Implícito:

(resolver o sistema em cada passo)
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
22
Funções MATLAB para Sistemas Rígidos



ode15s: Método de Gear (Backward-Difference-Formulae)
ode23s: Método de Rosenbrock
Ode23ts: Método Implícito de Runge-Kutta
5 Junho 2012
Equações Diferenciais – Métodos Adaptativos e Rigidez
23