Método de Runge-Kutta 2ª e 4ª ordem

Runge-Kutta 2ª ordem

No método explícito de euler tínhamos:

$$\begin{array}{rl}{y}_{n+1}& =y_n+f(t_n,y_n)\mathrm{\Delta }t\end{array}$$

Sendo $\frac{dy}{dt}=f(t,y)$. Podemos reescrever como:

$$\begin{array}{rl}{y}_{n+1}& =y_n+a{k}_1\end{array}$$

Onde $a=1$ e $k_1=f(t_n,y_n)\mathrm{\Delta }t$. Agor se supormos uma solução:

$$\begin{array}{rl}{y}_{n+1}& =y_n+a{k}_1+b{k}_2\left(1\right)\end{array}$$ Com o termo adicional dependendo de uma posição genérica $y_n+cf(t_n,y_n)\mathrm{\Delta }t$ em um tempo genérico $t+d\mathrm{\Delta }t$, isto é $k_2=f(t_n+d\mathrm{\Delta }t,y_n+cf(t_n,y_n)\mathrm{\Delta }t)\mathrm{\Delta }t$. Usando o fato de que $k_1=f(t_n,y_n)\mathrm{\Delta }t$, podemos escrever então que:

• $k_1=f(t_n,y_n)\mathrm{\Delta }t$
• $k_2=f(t_n+d\mathrm{\Delta }t,y_n+c{k}_1)\mathrm{\Delta }t$

Agora lembrando a expansão em série de taylor que também vimos no método explícito e Euler:

$$y(t+\mathrm{\Delta }t)=y\left(t\right)+y^{\prime }\left(t\right)\mathrm{\Delta }t+y^{″}\left(t\right)\frac{\mathrm{\Delta }{t}^2}{2}+\sum _{n=3}^{\mathrm{\infty }}{y}^{\left(n\right)}\left(t\right)\frac{\mathrm{\Delta }{t}^n}{n!}$$

Abrindo a segunda derivada, temos:

$$\begin{array}{rl}{y}^{″}\left(t\right)=\frac{d^2}{d{t}^2}y\left(t\right)& =\frac{d}{dt}f(t,y)=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}\frac{dy}{dt}=\frac{\partial f}{\partial t}+f(t,y)\frac{\partial f}{\partial y}\end{array}$$

Substituindo então, e escrevendo apenas ${𝒪}\left({\mathrm{\Delta }}^3\right)={\sum }_{n=3}^{\mathrm{\infty }}{y}^{\left(n\right)}\left(t\right)\frac{\mathrm{\Delta }{t}^n}{n!}$, temos a seguinte expansão em série de Taylor:

$$y(t+\mathrm{\Delta }t)=y\left(t\right)+y^{\prime }\left(t\right)\mathrm{\Delta }t+(\frac{\partial f}{\partial t}+f(t,y)\frac{\partial f}{\partial y})\frac{\mathrm{\Delta }{t}^2}{2}+{𝒪}\left({\mathrm{\Delta }}^3\right)\left(2\right)$$

Vamos expandir $k_2$. Uma expansão de Taylor de primeira ordem para uma função de 2 variáveis em torno de $(a,b)$ é dado por :

$$f(x,y)\approx f(a,b)+f_x(a,b)(x-a)+f_y(a,b)(y-b)$$

Onde $f_a$ denota a derivada da função $f$ na variável $a$. Para o nosso caso, temos então para uma expansão em torno de $(x,y)$:

$$f(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y)\approx f(x,y)+f_x(x,y)\mathrm{\Delta }x+f_y(x,y)\mathrm{\Delta }y$$

Expandindo então $k_2$ em torno de $(t_n,y_n)$ temos:

$$k_2\approx [f(t_n,y_n)+d\mathrm{\Delta }t\frac{\partial }{\partial t}f(t_n,y_n)+c{k}_1\frac{\partial }{\partial y}f(t_n,y_n)]\mathrm{\Delta }t$$

Aqui podemos notar que $\mathrm{\Delta }t$ multiplica a expansão da função, então quando desprezamos os termos de segunda ordem da expansão de $f(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y)$, deprezamos os termos de terceira ordem de $k_2$. Substituindo então o $k_2$aproximado e $k_1$ na equação 1, temos:

$$y_{n+1}=y_n+af(t_n,y_n)\mathrm{\Delta }t+b[f(t_n,y_n)+d\mathrm{\Delta }t\frac{\partial }{\partial t}f(t_n,y_n)+cf(t_n,y_n)\mathrm{\Delta }t\frac{\partial }{\partial y}f(t_n,y_n)]\mathrm{\Delta }t$$

Manipulando:

$$y_{n+1}=y_n+(a+b)f(t_n,y_n)\mathrm{\Delta }t+[2bd\frac{\partial }{\partial t}f(t_n,y_n)+2bcf(t_n,y_n)\frac{\partial }{\partial y}f(t_n,y_n)]\frac{\mathrm{\Delta }{t}^2}{2}\left(3\right)$$

Comparando a aproximação 3 com a expansão 2 temos a seguinte relação:

$$\begin{array}{rl}a+b& =1\\ bd& =\frac{1}{2}\\ bc& =\frac{1}{2}\end{array}$$ Diferentes conjuntos de valore satisfazem este sistema. O método do ponto médio é obtido se ecolhermos: $d=c=\frac{1}{2}$, $b=1$ e $a=-1$:

• $k_2=f(t_n+\frac{\mathrm{\Delta }t}{2},y_n+\frac{k_1}{2})\mathrm{\Delta }t$

Então:

$$\begin{array}{rl}{y}_{n+1}& =y_n+f(t_n+\frac{\mathrm{\Delta }t}{2},y_n+f(t_n,y_n)\frac{\mathrm{\Delta }t}{2})\mathrm{\Delta }t\end{array}$$

O método de Heun é obtido se for escolhido $a=b=\frac{1}{2}$ e $c=d=1$:

• $k_1=f(t_n,y_n)\mathrm{\Delta }t$
• $k_2=f(t_n+\mathrm{\Delta }t,y_n+k_1\mathrm{\Delta }t)\mathrm{\Delta }t$

$$\begin{array}{rl}{y}_{n+1}& =y_n+(k_1+k_2)\frac{1}{2}\end{array}$$

Uma observação, é que o erro global no algoritmo de Runge-Kutta de segunda ordem é ${𝒪}\left({\mathrm{\Delta }}^2\right)$ e o local é ${𝒪}\left({\mathrm{\Delta }}^3\right)$.

Exemplo

import matplotlib.pyplot as plt            #Biblioteca para plotar gráficos
import numpy as np                         #Biblitoeca de cálculos científicos

#Taxas de variação
def fv(x,w2):      #Velocidade
  return (-w2*x)
def fx(v):
  return (v)       #Posição

#Constantes
m=1  ; k= 1.; w2= k/m
#Parâmetros
dt  = 0.0001 ; tau = 2*np.pi; tf=4*tau ; Np= int(tf/dt)
#Valores iniciais
x=[1]; v=[0]; t=[0]; E=[k*x[0]**2/2+m*v[0]**2/2]

#Método Range-Kutta de segunda ordem, no método do ponto médio
for it  in range(Np):
  #Posição
  k1 = fx(v[it])*dt
  k2 = fx(v[it]+k1/2)*dt
  x.append(x[it]+k2)
  #Velocidade
  k1 = fv(x[it],w2)*dt
  k2 = fv(x[it]+k1/2,w2)*dt
  v.append(v[it]+k2)
  #Energia
  E.append(k*x[it+1]**2/2+m*v[it+1]**2/2)
  #Tempo
  t.append(dt+it*dt)

plt.plot(t,x)
plt.plot(t,v)
plt.plot(t,E)
#plt.plot(x,v)

Runge-Kutta 4ª ordem

O método de Runge-Kutta de quarta ordem segue uma ideia similar e pode ser obtido utilizando a mesma técnica. Porém agora vamos ignorar termos de ordem $\mathrm{\Delta }{t}^5$ ou superior, então será necessário lidar com uma enorme quantidade de termos, o que torna a tarefa exaustiva e repetitiva. Logo não será feito esta demonstração aqui, mas o algoritmo de Runge-Kutta de quarta ordem pode ser dado por:

• $k_1=f(y_n,t_n)\mathrm{\Delta }t$
• $k_2=f(y_n+\frac{k_1}{2},t_n+\frac{\mathrm{\Delta }t}{2})\mathrm{\Delta }t$
• $k_3=f(y_n+\frac{k_2}{2},t_n+\frac{\mathrm{\Delta }t}{2})\mathrm{\Delta }t$
• $k_4=f(y_n+k_3,t_n+\mathrm{\Delta }t)\mathrm{\Delta }t$

E por fim, temos então que o novo valor será dado por:

$$y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$$

Exemplo

Vamos resolver o mesmo exemplo anterior, porém agora utilizando o Range-Kutta de quarta ordem.

import matplotlib.pyplot as plt            #Biblioteca para plotar gráficos
import numpy as np                         #Biblitoeca de cálculos científicos

#Taxas de variação
def fv(x,w2):      #Velocidade
  return (-w2*x)
def fx(v):
  return (v)       #Posição

#Constantes
m=1  ; k= 1.; w2= k/m
#Parâmetros
dt  = 0.00001 ; tau = 2*np.pi; tf=4*tau ; Np= int(tf/dt)
#Valores iniciais
x=[1]; v=[0]; t=[0]; E=[k*x[0]**2/2+m*v[0]**2/2]

#Método Range-Kutta de quarta ordem
for it  in range(Np):
  #Posição
  k1 = fx(v[it])*dt
  k2 = fx(v[it]+k1/2)*dt
  k3 = fx(v[it]+k2/2)*dt
  k4 = fx(v[it]+k3)*dt
  x.append(x[it]+(k1+2*k2+2*k3+k4)/6)
  #Velocidade
  k1 = fv(x[it]     ,w2)*dt
  k2 = fv(x[it]+k1/2,w2)*dt
  k3 = fv(x[it]+k2/2,w2)*dt
  k4 = fv(x[it]+k3  ,w2)*dt
  v.append(v[it]+(k1+2*k2+2*k3+k4)/6)
  #Energia
  E.append(k*x[it+1]**2/2+m*v[it+1]**2/2)
  #Tempo
  t.append(dt+it*dt)

plt.plot(t,x)
plt.plot(t,v)
plt.plot(t,E)
#plt.plot(x,v)

Principais materiais utilizados

1. Métodos de Runge-Kutta explícitos (REAMAT, UFRGS)
2. Runge-Kutta Methods (Michael Zeltkevic, Instituto de Tecnologia de Massachusetts)
3. Second Order Runge-Kutta (Erik Cheever, Swarthmore)