Introdução a Resolução Numérica de Equações Diferenciais Ordinárias

Introdução a Resolução Introdução a Resolução

Numérica de EquaçõesNumérica de Equações

Diferenciais Ordinárias Diferenciais Ordinárias

Definição 6.1.1: Uma equação que envolve derivadas até ordem n, é chamada de equação diferencial ordinária (EDO) de ordem n e pode ser escrita na forma:

(6.1.1)

Definição 6.1.2: A solução da equação (6.1.1) é qualquer função y = F(x) que é definida em [a,b] e tem n derivadas neste intervalo e que satisfaz (6.1.1).

y x f x y x y x y xn n( ) ( )( ) ( , ( ), ( ), , ( )) 1

bxa

Equações Diferenciais

A forma mais simples de uma EDO é

(6.1.2)

onde f é contínua para a < x < b.

A solução geral desta equação é:

(6.1.3)

com constante c determinada por

)(xfy

y x f x dx c( ) ( ) y x y( )0 0

De um modo geral temos

(6.1.4)

como por exemplo:

é conveniente reduzi-la a um sistema de EDO de primeira ordem chamando :

nn

nn

ayayaybxa

yyyxfxy

)(,,)(,)(

),,,,()(

121

)1(

1)0(,,0)0(,1)0(10

1)()()( 2

yyyx

fxxyexyxxy x

y y x1 ( )

1)0(,,0)0(,1)0(10

1)()()( 2

yyyx

fxxyexyxxy x

y x y

y x y

y x xy x e y x x f

y y y

x

1 2

2 3

3 2 12

1 2 3

1

0 1 0 0 0 1

( )

( )

( ) ( ) ( )

( ) , ( ) , ( )

Para reduzir a uma EDO de primeira ordem assumimos:

y y x1 ( )

10 x

)(2 xyy )(3 xyy

isto é, se

3

2

1

~

y

y

y

y ~

y

y

y

y

1

2

3

~( )

( )

( )

( )

y

y

y

y

0

0

0

0

1

2

3

~

( , , , )

~( , )

y

y

y

f x y y y

F x y2

3

1 2 3

de um modo geral (6.1.6)

dy

dtF x y x

y y

~ ~( , ~( ))

~( ) ~

0 0

Equações de Diferenças

Definição 6.2.1: Uma equação de diferenças de ordem n é uma sequência de equações da forma:

(6.2.1)

Os são funções de (n+1) variáveis e os valores são específicos.

Definição 6.2.2: Uma solução de (6.2.1) é uma sequência

que satisfaz (6.2.1).

g y y y

y

k n k n k

i i

( , , , )

1 0

k

i n

0 1

0 1 1

, ,...

, ,...,

gk i

y y y y y yn n n0 1 2 1 1, , ,..., , , ,...

Uma forma especial de (6.2.1) é:

(6.2.2)

n k n n k n k

i i

y y y

y

1 1 0 0...

k

i n

0 1

0 1 1

, ,...

, ,...,

Estas fórmulas nos lembram as fórmulas recursivas e pode ser mostrado que toda fórmula recursiva pode ser escrita como uma equação de diferença.

Suponha, por exemplo, que a equação (6.2.2) tenha uma solução da forma:

k=0,1,... (6.2.5)

onde k é constante a ser determinada, então de (6.2.2) temos:

Como 0 então p() = 0

(6.2.6)

Logo se é uma raiz de p() = 0 então satisfaz a equação (6.2.2).

ykk

nk n

nk n k

1

10... k

nn

nn( ... )

1

10 0

p nn

nn( ) ...

1

10

Exemplo : (6.2.7)

(6.2.8)

satisfaz (6.2.7). Note que, se e então

e a solução será (6.2.9)

y y yk k k 2 185 4.

p( ) 2 17 8 02

1 05 . 2 8

y a akk k 1 1 2 2

y0 2 11 y

a a

a a1 2

1 2

2

05 8 1

.

a

a1

2

2

0

ykk2 05( . )

Se em vez dos valores iniciais anteriores tivermos, por exemplo:

então teríamos , de modo que a solução

perturbada seria :

(6.2.9-a)

  Mesmo para valores moderados de k, a mudança é grande.

y

y0

1

2 000001

1000008

.

.

a

a1

26

2 0

10

.

ykk k 2 05 10 86( . )

Teorema 6.2.1: Teorema 6.2.1: A solução da equação de diferenças

que envolve somente as raízes , , ... , do

polinômio característico p() = 0 é estável (bem

condicionada numéricamente) se todas as outras

raízes , , ... , satisfazem <1

para .

1 2 m

1m 2m n k

k m n 1,...,

Estudaremos, agora, os métodos que aproximam uma EDO por uma equação de diferença. Para facilitar, seja inicialmente

 (6.2.10)

• Em alguns métodos para determinar , uma aproximação de levamos em consideração apenas o valor de f no ponto , isto é, temos métodos de um só passo.

• Quando estamos calculando podemos levar em conta, não só , mas também os passos anteriores , ... , , etc. Tais métodos são chamados de métodos de passos múltiplos.

yk 1

y xk( )

yk 1

yk

yk 1ky lky

y f x y x

y x y

( , ( ))

( )0 0

A forma geral dos métodos de k passos é dada por:

(6.2.11)

Se Método Explícito.

Se Método Implícito.

y y h fi j i j jj

k

i jj

k

1 10

11

0 00 0

Métodos de um PassoMétodos de um PassoMétodo de Runge Kutta (RK)Método de Runge Kutta (RK)

Os métodos de RK são obtidos pela série de Taylor em que se omitem os termos de mais alta ordem na expansão. Assim, expandindo a solução em série de Taylor obtemos:

...)(!3

)(!2

)()()(32

xyX

xyX

xyXxyXxy

Sabendo que

ydy

dt

df

dt

f

xx

f

tf

f

x

f

t

ydy

dtf

f

xf

f

xf

f

x t

f

x

f

t

f

t

22

2

2 2 2

22

obtemos a base para todos os métodos RK

(6.2.14)

t

f

x

ff

XXfxyXxy

!2)()(

2

Xf

f

xf

f

xf

f

x t

f

x

f

t

f

t

32

2

2

2 2 2

232

!...

Discretizando, temos:

logo

(6.2.17)

levando-se em conta

y x y x y k

y x X y x y k

f y x x f y x x f k

k

k

k k

( ) ( ) ( )

( ) ( ) ( )

( ( ), ) ( ( ), ) ( )

1 1

f

t0

x

fkf

XkXfkyky )(

!2)()()1(

2

Xf k

f k

xf k

f k

x

32

2

2

2

3!( )

( )( )

( )...

Método de Euler

A menor aproximação é da forma :

k = 0,1,... (6.2.18)

e é conhecido como MÉTODO DE EULER..

y k y k Xf k( ) ( ) ( ) 1

Runge Kutta de segunda ordemPara desenvolvermos o método de segunda ordem RK2, assumiremos que a solução assume uma aproximação tendo a expressão:

k = 0,1,... (6.2.19)

 onde e são constantes e

(6.2.20)

(6.2.21)

com constante.

y k y k c g c g( ) ( ) 1 1 1 2 2

c1 2c

g Xf k1 ( )

g Xf y k g2 2 1 ( ( ) ) 2

Desenvolvendo em série de Taylor a expressão (6.2.21), obtem-se:

(6.2.23)

Logo então

y k y k c c Xf k c X f kf k

x( ) ( ) ( ) ( ) ( )

( )... 1 1 2 2 2

2

c c

c

1 2

2 2

11

2

c c1 2

2

12

1

y k y k g g( ) ( ) ( ) 11

2 1 2

g Xf k1 ( ) g Xf y k g2 1 ( ( ) )

Runge Kutta de quarta ordemRunge Kutta de quarta ordem

Normalmente X = h, por universalidade : e , e seguindo o mesmo procedimento obteremos RK4 (o mais importante):

k = 0,1,...

(6.2.25)

g K1 122 Kg

y k y k g g g g( ) ( ) ( ) 11

62 21 2 3 4

g Xf k

g Xf y k g

g Xf y k g

g Xf y k g

1

2 1

3 2

4 3

0 5

0 5

( )

( ( ) . )

( ( ) . )

( ( ) )

Nos métodos RK temos as seguintes características:

São auto-inicializáveis, isto é, a partir da condição inicial temos condições de calcular os demais pontos;

Não precisam do cálculo manual de derivadas;

Permitem fácil troca de X=h;

São facilmente codificáveis;

Utilizam p avaliações de função, se a ordem da fórmula for p.

Método de Passos MúltiplosMétodo de Passos Múltiplos

Conforme vimos, os métodos de passo simples precisam de informação sobre a solução apenas em x = para achar uma aproximação para y( + h ).

A característica dos métodos de passo múltiplos é que eles usam informações sobre a solução em mais de um ponto. Inicialmente vamos supor que conhecemos aproximações para y(x) em e

, i = 0,1,... .

A seguir exporemos aqui uma classe de métodos de passo múltiplo que é baseado no principio de integração numérica (Métodos Adams-Bash Forth).

xn xn

x x xn0 1, ,...,x x hi i 1

Métodos ExplícitosMétodos Explícitos

São obtidos quando trabalhamos com para aproximar a integral :

Aproximamos pelo polinômio de grau m, que interpola em e então:

x x xn n n m, ,..., 1

y x y x f x y x dxn n

x

x

n

n

( ) ( ) ( , ( ))

1

1

f x y x( , ( )) P xm ( )f x y( , ) x x xn n n m, ,..., 1

y x y x P x dxn n

x

x

m

n

n

( ) ( ) ( )

1

1

Escolhendo m = 3 então vamos usar ; ; ; aproximando pelo polinômio de grau 3 ( ) que interpola nos pontos acima, fazendo as contas:

Se o polinômio fosse do quarto grau

),( nn yx ),( 11 nn yx),( 22 nn yx ),( 33 nn yx f x y x( , ( ))

P x3 ( ) f x y x( , ( ))

y yh

f f f fn n n n n n 1 1 2 32455 59 37 9[ ]

y yh

f f f f fn n n n n n n 1 1 2 3 47201901 2774 2616 1274 251[ ]

Métodos ImplícitosMétodos Implícitos

São obtidos quando trabalhamos com (Adams-Moulton). Se m = 2 vamos usar

Da mesma forma que fizemos anteriormente

x x xn n n m, ,..., 1

2111

31

519924

)()()(1

nnnnnn

x

xnn

ffffh

yy

dxxPxyxyn

n

),( ),,( ),,( ),,( 221111 nnnnnnnn yxyxyxyx

Notamos para métodos de mesma ordem :

  eles usam menos informações que os métodos explícitos;

os coeficientes de são menores e portanto são fórmulas menos sensíveis ao arredondamento;

os limites do erro de truncamento são menores;

o problema de cálculo de , pois precisamos de , o que torna o método inaplicável numéricamente.

f i

yi1

f i1

Métodos Previsores CorretoresMétodos Previsores Corretores

Para contornar o último quesito podemos adotar um esquema:

1.    Calculamos o valor de por um método explícito e o chamamos de valor previsto.

2.    Com o valor previsto, calculamos .

3.    Recalculamos , isto é, corrigimos o valor de já previsto usando uma fórmula implícita.

yi1

f i1

yi1yi1

MatlabMatlabA biblioteca do Matlab de EDOs os seguintes métodos de valor inicial:

ode23 método explícito de um passo, RK de ordem baixa.

ode45 método explícito de um passo, RK de ordem média. Este é geralmente, o primeiro método a se tentar em um novo problema.

ode113 método de passo múltiplo, de Adams-Bashforth-Moulton de ordens variadas.

ODE 23Adequado para problemas que apresentam bruscas

variações na solução, para os quais é aceitável uma baixa precisão, ou problemas em que f(t, y) não é suave, ou seja, descontínua.

Equação da MolaEquação da Molasujeita a uma forçasujeita a uma força )(tasin

function wd=spring(t,w);

a=2.0; m=2.0;

c=1.4;h=0.1;

wd = [w(2);-c/m*w(2)-h/m*w(1)+a/m*sin(t)];

)(122

21

tsinm

aw

m

hw

m

cw

ww

>> tspan=[0 100];>> wo=[0;0];>> [t,w] = ode23(@mola,tspan,wo);>> plot(t,w(:,1));>> grid>> xlabel('tempo, segundos');>> ylabel('Elongamento da mola, metros');

>> plot(w(:,1),w(:,2));

>> xlabel('Velocidade das oscilacoes da mola');

>> ylabel('Amplitude das oscilacoes da mola');

ODE 45Adequado para problemas não-stiff que exijam

precisão moderada.

function yprime=vdpol(t,y)

mu=2;

yprime = [y(2);mu*(1-y(1)^2)*y(2)-y(1)];

Equação de Van der PolEquação de Van der Pol

0)1( 2 xxxx onde é um parâmetro positivo. Escolhendo

e a equação de Van der Pol se torna:

xy 1

xy 2

12212

21

)1( yyyy

yy

>> tspan = [0 20];

>> yo = [2; 0];

>> [t,y] = ode45(@vdpol,tspan,yo);

>> plot(t,y(:,1),t,y(:,2),'--');

>> xlabel('tempo');

>> title('Solucao de Van der Pol');