Use o botão abaixo para reportar erros ou dar sugestões.

Cálculo Numérico - Versão Python

10.7 Métodos de Runge-Kutta explícitos


Nas seções anteriores, exploramos os métodos de Euler e Euler modificado para resolver problemas de valor inicial. Neste momento, deve estar claro ao leitor que o método de Euler melhorado produz soluções de melhor qualidade que o método de Euler para a maior parte dos problemas estudados. Isso se deve, conforme Seção 10.6, à ordem de precisão do método.

Os métodos de Runge-Kutta generalizam o esquema do método de Euler melhorado, inserindo mais estágios de cálculo e buscando ordens de precisão mais altas. Nesta seção, trataremos dos métodos de Runge-Kutta explícitos2

Para tal, considere o problema de valor inicial: u(t) = f(t,u(t)), (10.135) u(t0) = a. (10.136)

Integrando a EDO em [t(n),t(n+1)] obtemos u(n+1) = u(n) +t(n)t(n+1) f(t,u(t))dt (10.137)

O método de Euler aproxima a integral no lado direito da expressão acima utilizando apenas o valor de f(t,u) em t = t(n), o método de Euler melhorado aproxima esta integral utilizando os valores de f(t,u) t = t(n) e t = t(n+1). Os métodos de Runge-Kutta explícitos admitem estágios intermediários, utilizando outros valores de f(t,u) nos pontos {τ1,τ2,,τν} dentro do intervalo [tn,t(n+1)]. Veja esquema abaixo:

un un+1
| |





t(n) t(n+1)
τ1 τ2 τν
t(n)      t(n) + c 2h     t(n) + c νh  

Observe que τj = t(n) + c jh com 0 = c1 c2 cν 1, isto é, o primeira ponto sempre coincide com o extremo esquerdo do intervalo, mas o último ponto não precisa ser o extremo direito. Ademais, um mesmo ponto pode ser usado mais de uma vez com aproximações diferentes para u(τj).

Desta forma, aproximamos a integral por um esquema de quadratura com ν pontos: t(n)t(n+1) f(t,u(t))dt un + h j=1νb jf τj,u(τj) (10.138)

onde bj são os pesos da quadratura numérica que aproxima a integral. Assim como na construção do método de Euler melhorado, não dispomos dos valores u(τj) antes de calculá-los e, por isso, precisamos estimá-los com base nos estágios anteriores:

ũ1 = u(n) ũ2 = un + ha21k1 ũ3 = un + h a31k1 + a32k2 ũ4 = un + h a41k1 + a42k2 + a43k3 ũν = un + h aν1k1 + aν2k2 + + aννkν u(n+1) = u n + h[b1k1 + b2k2 + + bνkν], (10.139)

onde kj = f τj,ũj, A := aij é a matriz Runge-Kutta (triangular inferior com diagonal zero), bj são os pesos Runge-Kutta e cj são os nós Runge-Kutta. Estes coeficientes podem ser escritos de forma compacta em uma tabela conforme a seguir:

c
 
A


bT
    =   
c1 0 0 0
c2 a21 0 0
c3 a31 a22 0




b1 b2 b3.
    =   
c1
c2 a21
c3 a31 a22




b1 b2 b3.

Na tabela mais à direita, omitimos os termos obrigatoriamente nulos.

Exemplo 10.7.1. O método de Euler modificado pode ser escrito conforme: ũ1 = u(n) (10.140) ũ2 = u(n) + hk 1 (10.141) u(n+1) = u(n) + h 1 2k1 + 1 2k2 (10.142) (10.143)

Identificando os coeficientes, obtemos ν = 2, c1 = 0, c2 = 1, a21 = 1, b1 = 1 2 e b2 = 1 2. Escrevendo na forma tabular, temos:

c
A



b
   =   
0
1 1



1 2 1 2

Observação 10.7.1. Nos métodos chamados explícitos, os elementos da diagonal principal (e acima dela) da matriz A devem ser nulos, pois a aproximação em cada estágio é calculada com base nos valores dos estágios anteriores. Nos métodos implícitos, essa restrição é removida, neste caso, o cálculo da aproximação em um estágio envolve a solução de uma equação algébrica.

Observação 10.7.2. Além da condição fixa c1 = 0, para que um método seja de ordem pelo menos unitária, isto é, p 1.

10.7.1 Métodos de Runge-Kutta com dois estágios


Os métodos de Runge-Kutta com dois estágios ν = 2 são da seguinte forma: ũ1 = u(n) (10.144) ũ2 = u(n) + ha 21k1 (10.145) u(n+1) = u(n) + h[b 1k1 + b2k2], (10.146)

onde k1 = f(t(n),u(n)) e k2 = f(t(n) + c 2h,ũ2)

Assumindo suavidade suficiente em f, usamos o polinônio de Taylor: k2 = f(t(n) + c 2h,ũ2) (10.147) = f(t(n) + c 2h,u(n) + a 21hk1) (10.148) = f(t(n),u(n)) + h c 2f t + a21k1f u + O(h2) (10.149) = k1 + h c2f t + a21k1f u + O(h2) (10.150)

fazendo com que (10.146) se torne u(n+1) = u(n) + hb 1k1 + hb2k2 (10.151) = un + hb1k1 + hb2 k1 + h c2f t + a21k1f u + O(h2) (10.152) = un + h(b1 + b2)k1 + h2b 2 c2f t + a21k1f u + O(h3) (10.153)

Usando a equação diferencial ordinária que desejamos resolver e derivando-a em t, obtemos: u(t) = f(t,u(t)), (10.154) u(t) = f t + f uu(t) = f t + f uf(t,u(t)). (10.155)

Agora, expandimos em série de Taylor a solução exata u(t) em t = t(n), u(t(n+1)) = u(t(n) + h) = u(n) + hu(t) + h2 2 u(t) + O(h3) (10.156) = u(n) + hf(t,u(n)) + h2 2 d dtf(t,u(t)) + O(h3) (10.157) = u(n) + hf(t,u(n)) + h2 2 f t + f uu(t(n)) + O(h3) (10.158) = u(n) + hf(t,u(n)) + h2 2 f t + f uf(t,u(n)) + O(h3) (10.159) = u(n) + hf(t,u(n)) + h2 2 f t + f uk1 + O(h3) (10.160)

Finalmente comparamos os termos em (10.153) e (10.160) de forma a haver concordância na expansão de Taylor até segunda ordem, isto é, restanto apenas o erro de ordem 3 e produzindo um método de ordem de precisão p = 2:

b1 + b2 = 1,b2c2 = 1 2ea21 = c2 (10.161)

Este sistema é formada por três equações e quatro incógnitas, pelo que admite infinitas soluções. Para construir toda a família de soluções, escolha um parâmetro α (0,1] e defina a partir de (10.161):

b1 = 1 - 1 2α,b2 = 1 2α,c2 = α,a21 = α (10.162)

Portanto, obtemos o seguinte esquema genérico:

c
A



b
   =   
c 2 a21



b 1 b2
   =   
α α



1 - 1 2α 1 2α
,   0 < α 1.

Algumas escolhas comuns são α = 1 2, α = 2 3 e α = 1:

0
1 2 1 2



0 1
,
0
2 3 2 3



1 4 3 4
   e   
0
1 1



1 2 1 2
.

Note que a tabela da direita fornece o método Euler modificado α = 1. O esquema iterativo assume a seguinte forma: ũ1 = u(n) (10.163) ũ2 = u(n) + hαk 1 (10.164) u(n+1) = u(n) + h 1 - 1 2αk1 + 1 2αk2 , (10.165)

onde k1 = f(t(n),u(n)) e k2 = f(t(n) + hα,ũ 2). Ou, equivalentemente: k1 = f(t(n),u(n)) (10.166) k2 = f(t(n) + hα,u(n) + hαk 1) (10.167) u(n+1) = u(n) + h 1 - 1 2αk1 + 1 2αk2 , (10.168)

10.7.2 Métodos de Runge-Kutta com três estágios


Os métodos de Runge-Kutta com 3 estágios podem ser descritos na forma tabular como:

0
c 2 a21
c3 a31 a32




b1 b2 b3

Seguindo um procedimento similar ao da Seção 10.7.1, podemos obter as condições equivalentes às condições (10.161) para um método com ν = 3 e ordem p = 3, as quais são: b1 + b2 + b3 = 1, (10.169) b2c2 + b3c3 = 1 2, (10.170) b2c22 + b 3c32 = 1 3, (10.171) b3a32c2 = 1 6, (10.172) a21 = c2, (10.173) a31 + a32 = c3. (10.174)

Assim, temos 6 condições para determinar 8 incógnitas, o que implica a existência de uma enorme família de métodos de Runge-Kutta com três estágios e ordem p = 3. Se fixarmos os coeficientes c2 e c3, podemos os outros de forma única: b1 = 1 - 1 2c2 - 1 2c3 + 1 3c2c3 (10.175) b2 = 3c3 - 2 6c2(c3 - c2) (10.176) b3 = 2 - 3c2 6c3(c3 - c2) (10.177) a21 = c2 (10.178) a31 = c3 - 1 6b3c2 (10.179) a32 = 1 6b3c2 (10.180)

Alguns exemplos de métodos de Runge-Kutta de 3 estágios são o método clássico de Runge-Kutta c2 = 1 2 e c3 = 1 e o método de Nyström c2 = c3 = 2 3:

0
1 2 1 2
1 - 1 2




1 6 4 6 1 6
   e      
0
2 3 2 3
2 3 0 2 3




2 8 3 8 3 8
.

10.7.3 Métodos de Runge-Kutta com quatro estágios


As técnicas utilizadas nas seções 10.7.1 e 10.7.2 podem ser usadas para obter métodos de quarta ordem (p = 4) e quatro estágios (ν = 4). As seguintes tabelas descrevem os dois esquemas mais conhecidos de Runge-Kutta quarta ordem com quatro estágios. O primeiro é denominado método de Runge-Kutta 3/8 e o segundo é chamado de método de Runge-Kutta clássico.

0
1 3 1 3
2 3 - 1 3 1
1 1 - 1 1





1 8 3 8 3 8 1 8
   e       
0
1 2 1 2
1 2 0 1 2
1 0 0 1





1 6 2 6 2 6 1 6

O método de Runge-Kutta clássico é certamente o mais notório dos métodos de Runge-Kutta e seu esquema iterativo pode ser escrito como a seguir: k1 = hf t(n),u(n) (10.181) k2 = hf t(n) + h2,u(n) + hk 12 (10.182) k3 = hf t(n) + h2,u(n) + hk 22 (10.183) k4 = hf t(n) + h,u(n) + hk 3 (10.184) u(n+1) = u(n) + k1 + 2k2 + 2k3 + k4 6 (10.185)

A seguinte heurística, usando o método de Simpson para quadratura numérica, pode ajudar a compreender os estranhos coeficientes: u(t(n+1)) - u(t(n)) = t(n)t(n+1) f(t,u(s))ds (10.186) h 6 f t(n),u(t(n)) + 4f t(n) + h2,u(t(n) + h2) (10.187) + f t(n) + hu(t(n) + h) (10.188) hk1 + 4(k2+k3 2 ) + k4 6 (10.189)

onde k1 e k4 representam os valores de f(t,u) nos extremos; k2 e k3 são duas aproximações diferentes para a inclinação no meio do intervalo.

Exercícios resolvidos


ER 10.7.1. Construa o esquema iterativo o método clássico de Runge-Kutta três estágios cuja tabela é dada a seguir:

0
1 2 1 2
1 - 1 2




1 6 4 6 1 6

Solução. k1 = f t(n),u(n) (10.190) k2 = f t(n) + h2,u(n) + k 12 (10.191) k3 = f t(n) + h,u(n) - k 1 + 2k2 (10.192) u(n+1) = u(n) + hk1 + 4k2 + k4 6 (10.193)

ER 10.7.2. Utilize o método clássico de Runge-Kutta três estágios para calcular o valor de u(2) com passos h = 10-1 e h = 10-2 para o seguinte problema de valor inicial: u(t) = -u(t)2 + t, (10.194) u(0) = 0. (10.195)

Aplicando o processo iterativo obtido no Problema Resolvido 10.7.1, obtemos a seguinte rotina:

 def f(t,u):  
return t-u**2  
 
 
def RK3_classico(h,Tmax,u1):  
   itmax = Tmax/h;  
u=np.empty(itmax+1)  
u[0]=u1  
 
for i in np.arange(0,itmax):  
t=i*h  
 
k1 = f(t,     u[i])  
k2 = f(t+h/2, u[i] + h*k1/2)  
k3 = f(t+h,   u[i] + h*(2*k2-k1))  
 
u[i+1] = u[i] + h*(k1+4*k2+k3)/6  
return u  
 
 
 
Tmax=2 #tempo maximo de simulacao  
u1=0 #condicoes iniciais na forma vetorial  
h=1e-2 #passo  
 
sol=RK3_classico(h,Tmax,u1);  
itmax=Tmax/h  
print(sol[itmax])

Exercícios


E 10.7.1. Aplique o esquema de Runge-Kutta segunda ordem com dois estágios cujos coeficientes são dados na tabela a seguir

0
2 3 2 3



1 4 3 4
para resolver o problema de valor inicial dado por: x(t) = sen (x(t)), (10.196) x(0) = 2. (10.197)

para t = 2 com h = 1e - 1, h = 1e - 2 e h = 1e - 3. Expresse sua resposta com oito dígitos significativos corretos.

Resposta. 2,9677921, 2,9682284 e 2,9682325.

E 10.7.2. Resolva pelo método de Euler, Euler melhorado, Runge-Kutta clássico três estágios e Runge-Kutta clássico quatro estágios o problema de valor inicial tratados nos exercícios resolvidos 10.2.1 e 10.3.1 dado por: u(t) = -0,5u(t) + 2 + t (10.198) u(0) = 8 (10.199)

Usando os seguinte passos: h = 1, h = 10-1, h = 10-2 e h = 10-3 e compare a solução aproximada em t = 1 com as soluções obtidas com a solução exata dada por:

u(t) = 2t + 8e-t2u(1) = 2 + 8e-12 6,85224527770107 (10.200)

Resposta.






Euler 6,0000000 6,7898955 6,8461635 6,8516386





εrel 1,2e-01 9,1e-03 8,9e-04 8,9e-05





Euler mod, 7,0000000 6,8532949 6,8522554 6,8522454





εrel 2,2e-02 1,5e-04 1,5e-06 1,5e-08





RK3 6,8333333 6,8522321 6,8522453 6,8522453





εrel 2,8e-03 1,9e-06 1,9e-09 1,8e-12





RK4 6,8541667 6,8522454 6,8522453 6,8522453





εrel 2,8e-04 1,9e-08 1,9e-12 1,3e-15





E 10.7.3. Aplique o método de Euler, o método de Euler melhorado, o método clássico de Runge-Kutta três estágios e o método clássico de Runge-Kutta quatro estágios para resolver o problema de valor inicial dado por u = u + t (10.201) u(0) = 1 (10.202)

com passo h = 1, h = 10-1, h = 10-2 e h = 10-3 para obter aproximações para u(1). Compare com a solução exata dada do problema dada por u(t) = 2et - t - 1 através do erro relativo e observe a ordem de precisão do método. Expresse a sua resposta com oito dígitos significativos para a solução e 2 dígitos significativos para o erro relativo.

Resposta.






Euler 2,0000000 3,1874849 3,4096277 3,4338479





εrel 4,2e-01 7,2e-02 7,8e-03 7,9e-04





Euler mod 3,0000000 3,4281617 3,4364737 3,4365628





εrel 1,3e-01 2,4e-03 2,6e-05 2,6e-07





RK3 3,3333333 3,4363545 3,4365634 3,4365637





εrel 3,0e-02 6,1e-05 6,5e-08 6,6e-11





RK4 3,4166667 3,4365595 3,4365637 3,4365637





εrel 5,8e-03 1,2e-06 1,3e-10 1,2e-14





Creative Commons License Este texto é disponibilizado nos termos da licença Creative Commons Atribuição-CompartilhaIgual 3.0 Não Adaptada (CC-BY-SA 3.0). Página gerada em 15/5/2019 às 15:24:50.

Informe erros ou edite você mesmo!