Método de Runge-Kutta

Método de Runge-Kutta

PMR 2420 – Mecânica Computacional

CAPÍTULO III MÉTODOS DE RUNGE-KUTTA
São métodos de passo simples requerem apenas derivadas de primeira ordem e pode fornecer aproximações precisas com erros de truncamento da ordem de h2, h3, h4, etc. Todas os métodos de Runge-Kutta têm a seguinte forma geral: yi +1 = yi + hφ( xi , yi , h ) onde φ, chamado de função incremento, é uma aproximação conveniente para f((x,y) no intervalo xi ≤ x ≤ xi +1 Runge-Kutta de segunda ordem Seja φ uma média ponderada de duas aproximações da derivadas k1 e k2 no intervalo xi ≤ x ≤ xi +1

φ = ak1 + bk2
Então o algoritmo de Runge-Kutta de segunda ordem fica:

yi +1 = yi + h( ak1 + bk2 )
y

(I)

k y i+1 k

2 1

y h x i x

i+1

 f ( xi , y i ) + f ( xi +1 , yi +1 )  yi +1 = yi + h   2   Assumindo: k1 = f ( xi , yi ) (a=b=1/2)

i+1

x

k 2 = f ( xi + ph, y i + qhf ( xi , y i )) = f ( xi +1 , y i +1 ) Expandindo k2 em série de Taylor:
33

PMR 2420 – Mecânica Computacional

k 2 = f ( xi + ph, y i + qhf ( xi , y i )) = f ( xi , y i ) + phf x ( xi , y i ) + qhf ( xi , y i ) f y ( xi , y i ) + O (h 2 ) (II) Substituindo (II) em (I) yi +1 = yi + h[ af ( xi , yi ) + bf ( xi , yi )] + h 2 [bpf x ( xi , yi ) + bqf ( xi , yi ) f y ( xi , yi )] + O ( h 3 ) (III) Expandindo y(x) em xi usando série de Taylor: h2 h3 f ' ( xi , y ( xi )) + f " (ξ , y (ξ )) = 2! 3! dy h 2 d 2 y h 3 d 3 y (ξ ) = y ( xi ) + h + + 3! dx 3 dx 2! dx 2 y ( xi + h ) = y ( xi ) + hf ( xi , y ( xi )) + (IV) Lembrando que: f '( xi , y ( xi )) = f x ( xi , y ( xi )) + f y ( xi , y ( xi )) f ( xi , y ( xi )) (regra da cadeia)

Igualando os termos de mesma potência em h das equações (III) e (IV), tem-se: potência de h expansão de y(x) em Taylor algoritmo de Runge-Kutta

0

y(xi)

yi

1

f ( xi , y ( xi ))
1 f x ( xi , y ( xi )) + f y ( xi , y ( xi )) f ( xi , y ( xi )) 2

( a + b ) f ( xi , y ( xi )) bpf x ( xi , yi ) + bqf ( xi , yi ) f y ( xi , yi )

2

assumindo: yi = y ( xi ) a + b = 1⇒ a = 1− b 1 bp =  2 ⇒ p=q= 1 1 2b bq =  2 família de métodos de segunda ordem

34

PMR 2420 – Mecânica Computacional

Escolhendo b =

1 1 ⇒ p = q =1 e a = 2 2

(Método de Heun)

k1 = f ( xi , yi ) k 2 = f ( xi + h, yi + hk1 )
então: 1  1 y i +1 = y i + h k 1 + k 2  2 2  fazendo b = 1 ⇒ a = 0 e p = q = 1/2

yi +1 = yi + hk2
1 1   onde: k1 = f ( xi , yi ) e k 2 = f  x i + h, y i + hk 1    2 2 isto pode ser interpretado da seguinte forma: h   h y i +1 = y i + hf  x i + , y i +1 2  onde y1+1 2 = yi + f ( xi , yi )   2 2
y y y i+1/2 h/2 h/2 x i x i+1 x k i+1 2 1

k

Este método é chamado de Método de Euler Modificado Runge-Kutta de 2a ordem tem erro local de truncamento O(h3) e erro global O(h2) Os métodos de Runge-Kutta de ordens superiores são desenvolvidos de modo análogo.

35

PMR 2420 – Mecânica Computacional

Runge-Kutta de 3a ordem y i +1 = y i + h(ak 1 + bk 2 + ck 3 ) onde k1, k2 e k3 são aproximações das derivadas em vários pontos no intervalo de integração xi , xi +1 . Nesse caso,

k1 = f ( xi , yi )
k 2 = f ( xi + ph , yi + phk1 )

k 3 = f ( xi + rh, yi + shk2 + ( r − s) hk1 )
Para determinar as constantes a, b, c, p, r e s, expandem-se k2 e k3 em Taylor, em torno de (xi,yi) e obtêm-se as seguintes equações: a +b+c =1 1 bp + cr = 2 bp 2 + cr 2 = cps = 1 6 1 2

Duas das contantes são escolhidas arbitrariamente. Para um determinado conjunto de constantes escolhido por Kutta, tem-se o seguinte método de 3a ordem h (k + 4k 2 + k 3 ) 6 1

y i +1 = y i +

k1 = f ( xi , yi ) h h   k 2 = f  xi + , yi + k 1   2 2 

k 3 = f (x i + h, yi + 2hk 2 − hk 1 ) dy = f (x ) dx

Equivale à regra de Simpson se

36

PMR 2420 – Mecânica Computacional

Runge-Kutta de 4a ordem (vários métodos dependendo da escolha das constantes). Segundo Kutta: h (k + 2 k 2 + 2 k 3 + k 4 ) 6 1 k1 = f ( xi , yi ) h h   k 2 = f  xi + , yi + k 1   2 2  h h   k 3 = f  xi + , yi + k 2   2 2  yi +1 = yi + k 4 = f ( xi + h, y i + hk 3 )

Controle do passo nos Algoritmos de Runge-Kutta Algoritmo de Runge-Kutta de ordem m → expansão em Taylor com termos até hm erro local de truncamento: onde: m = ordem do método K = função complicada dependendo de f(x,y) e suas derivadas Para escolher h dado um certo et considerando: a) erro local de truncamento et ≈ Kh m+1 (supor K constante) b) erro local de truncamento é a contribuição mais importante para o erro global
* Assumindo yn+1 como a solução exata, uma estimativa do erro local de truncamento pode ser obtida integrando-se entre xn e xn+1 com dois passos diferentes, h1 e h2, obtendo-se duas estimativas de yn+1, yn+1,1 e yn+1,2

et = Kh m+1 + O ( h m+ 2 )

Empregando a extrapolação de Richardson:
* y n +1 − y n +1,1 = Kh1 m +1

(xn+1 − xn )

h2 onde xn+1-xn equivale ao número de passos. Dividindo (1) por (2) tem-se:

* y n +1 − y n +1,2 = Kh2 m+1

(x

h1
n +1

(1) (2)

− xn )

37

PMR 2420 – Mecânica Computacional

y

* n +1

=

y n +1,1 − y n +1,2 (h1 h2 ) 1 − (h1 h2 )
m

m

(3)

escolhendo h2 = h1 2 tem-se:

y

* n +1

yn +1,1 − 2 m yn +1,2 = 1 − 2m

(4)

Assim, uma estimativa para o erro local de truncamento para a solução yn+1,1 , assumindo

(x

n +1

− x n ) = h1 é dada por: et = Kh1
m +1

=

2 m y n +1,2 − y n +1,1 2m − 1

(

)

(5)

Para o algoritmo de Runge-Kutta de 4a ordem (m = 4) 16 y − y n +1,1 15 n +1,2 Usando esse procedimento para monitorar o erro, o número de cálculos é triplicado. et = Kh1 5 = Um outro critério para determinar o tamanho do passo h, conhecido como critério de Collatz é feito através da avaliação da relação (k 3 − k 2 ) (k 2 − k 1 ) após cada passo de integração. Se esta relação se torna maior que alguns centésimos, o passo h deve ser diminuido. Solução de E.D.O. simultâneas seja o seguinte sistema com n E.D.O. dy1 = f 1 ( x , y1 , y 2 , , y n ) dx dy 2 = f 2 ( x , y1 , y 2 , , y n ) dx dy n = f n ( x , y1 , y 2 , dx com condições iniciais: y1 ( x0 ) = y10 y2 ( x0 ) = y2 0
£ ¢     ¡

(

)

ou

 dy   dx  = [f ]  

, yn )

ou

[y ] = [y 0 ]
38

PMR 2420 – Mecânica Computacional

yn ( x0 ) = yn0 Deve-se aplicar, qualquer um dos métodos apresentados, em paralelo em cada passo. Por exemplo, o método de Euler modificado: y1i +1 = y1i + hk 21 y2i+1 = y2 i + hk 2 2 y ni +1 = y ni + hk 2n onde:
 

ou

[y i +1 ] = [y i ] + h[k 2 ]

Exemplo: Aplicar R.K 4 à seguinte equação diferencial (problema resolvido no livro do Carnahan, página 367): d 2V dV + RC + V = 0 (circuito elétrico RLC) 2 dt dt sujeita às seguintes condições iniciais: LC V ( 0) = V0 (tensão elétrica no capacitor)

dV dV (corrente elétrica no circuito: i = − C ) =0 dt t = 0 dt Como primeiro passo, deve-se transformar a E.D.O. acima (segunda ordem) num sistema de E.D.O de primeira ordem. Fazendo:

¢

h h h  k 2n = f n  xi + , y1i + k 11 , y 2i + k 12 ,  2 2 2

¢

h h h  k 21 = f 1  xi + , y1i + k 11 , y 2i + k 12 ,  2 2 2 h h h  k 22 = f 2  xi + , y1i + k 11 , y 2i + k 12 ,  2 2 2
¢ ¤

¢

k 1n = f n xi , y1i , y 2i ,

(

¢

k 12
£

2

i

1i

, y 2i ,

¡

k11 = f 1 xi , y1i , y2i ,

( = f (x , y

y ni

) y )
ni

ou

[k 1 ] = [f ( xi , [y i ])]

y ni

)
y ni + h  k  2 1n  h  h h   y ni + k 1n  ou [k 2 ] = f ( xi + , [y i ] + [k 1 ])  2 2 2   y ni + h  k  2 1n 

39

PMR 2420 – Mecânica Computacional

dV =I dt podemos escrever a equação do circuito como:

dI + RCI + V = 0 dt Assim teremos o seguinte sistema de E.D.O. LC
1 dI R =− I− V dt L LC dV =I dt com condições iniciais: V ( 0) = V0 I ( 0) = 0 Aplicando Runge-Kutta de 4a ordem: I i +1 = I i + Vi +1 onde: h (k + 2k 2 + 2k 3 + k 4 ) 6 1 h = Vi + (q1 + 2q 2 + 2q 3 + q 4 ) 6

k 1 = f 1 (t i ,Vi , I i ) h h h   k 2 = f 1  t i + ,Vi + q1 , I i + k 1   2 2 2  h h h   k 3 = f 1  t i + ,Vi + q 2 , I i + k 2   2 2 2  k 4 = f 1 (t i + h,Vi + hq 3 , I i + hk 3 ) q1 = f 2 (t i ,Vi , I i ) h h h   q 2 = f 2  t i + ,Vi + q1 , I i + k 1   2 2 2  h h h   q 3 = f 2  t i + ,Vi + q 2 , I i + k 2   2 2 2  q 4 = f 2 (t i + h,Vi + hq 3 , I i + hk 3 )

40

PMR 2420 – Mecânica Computacional

sendo: f 1 (t ,V , I ) = − f 2 (t ,V , I ) = I 1 R I− V L LC

41

Comentários