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

Runge-Kutta 2ª ordem

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

$${\begin{aligned}y_{n+1}&=y_{n}+f\left(t_{n},y_{n}\right)\Delta t\end{aligned}}$$

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

$${\begin{aligned}y_{n+1}&=y_{n}+ak_{1}\end{aligned}}$$

Onde ${\textstyle a=1}$ e ${\textstyle k_{1}=f\left(t_{n},y_{n}\right)\Delta t}$. Agor se supormos uma solução:

$${\begin{aligned}y_{n+1}&=y_{n}+ak_{1}+bk_{2}\qquad \left(1\right)\end{aligned}}$$

Com o termo adicional dependendo de uma posição genérica ${\textstyle y_{n}+cf\left(t_{n},y_{n}\right)\Delta t}$ em um tempo genérico ${\textstyle t+d\Delta t}$, isto é ${\textstyle k_{2}=f\left(t_{n}+d\Delta t,y_{n}+cf\left(t_{n},y_{n}\right)\Delta t\right)\Delta t}$. Usando o fato de que ${\textstyle k_{1}=f\left(t_{n},y_{n}\right)\Delta t}$, podemos escrever então que:

• ${\textstyle k_{1}=f\left(t_{n},y_{n}\right)\Delta t}$
• ${\textstyle k_{2}=f\left(t_{n}+d\Delta t,y_{n}+ck_{1}\right)\Delta t}$

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

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

Abrindo a segunda derivada, temos:

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

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

$$y\left(t+\Delta t\right)=y\left(t\right)+y'\left(t\right)\Delta t+\left({\frac {\partial f}{\partial t}}+f\left(t,y\right){\frac {\partial f}{\partial y}}\right){\frac {\Delta t^{2}}{2}}+{\mathcal {O}}\left(\Delta ^{3}\right)\qquad \left(2\right)$$

Vamos expandir ${\textstyle k_{2}}$. Uma expansão de Taylor de primeira ordem para uma função de 2 variáveis em torno de ${\textstyle \left(a,b\right)}$ é dado por :

$$f\left(x,y\right)\approx f\left(a,b\right)+f_{x}\left(a,b\right)\left(x-a\right)+f_{y}\left(a,b\right)\left(y-b\right)$$

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

$$f\left(x+\Delta x,y+\Delta y\right)\approx f\left(x,y\right)+f_{x}\left(x,y\right)\Delta x+f_{y}\left(x,y\right)\Delta y$$

Expandindo então ${\textstyle k_{2}}$ em torno de ${\textstyle \left(t_{n},y_{n}\right)}$ temos:

$$k_{2}\approx \left[f\left(t_{n},y_{n}\right)+d\Delta t{\frac {\partial }{\partial t}}f\left(t_{n},y_{n}\right)+ck_{1}{\frac {\partial }{\partial y}}f\left(t_{n},y_{n}\right)\right]\Delta t$$

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

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

Manipulando:

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

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

$${\begin{aligned}a+b&=1\\bd&={\frac {1}{2}}\\bc&={\frac {1}{2}}\end{aligned}}$$

Diferentes conjuntos de valore satisfazem este sistema. O método do ponto médio é obtido se ecolhermos: ${\textstyle d=c={\frac {1}{2}}}$, ${\textstyle b=1}$ e ${\textstyle a=0}$:

• ${\textstyle k_{2}=f\left(t_{n}+{\frac {\Delta t}{2}},y_{n}+{\frac {k_{1}}{2}}\right)\Delta t}$

Então:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \begin{align} y_{n+1} & =y_{n}+f\left(t_{n}+\frac{\Delta t}{2},y_{n}+f\left(t_{n},y_{n}\right)\frac{\Delta t}{2}\right)\Delta t\end{align}}

O método de Heun é obtido se for escolhido Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\textstyle a=b=\frac{1}{2}} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\textstyle c=d=1} :

• Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\textstyle k_{1}=f\left(t_{n},y_{n}\right)\Delta t}
• Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\textstyle k_{2}=f\left(t_{n}+\Delta t,y_{n}+k_{1}\Delta t\right)\Delta t}

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \begin{align} y_{n+1} & =y_{n}+\left(k_{1}+k_{2}\right)\frac{1}{2}\end{align}}

Uma observação, é que o erro global no algoritmo de Runge-Kutta de segunda ordem é Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathcal{O}\left(\Delta^{2}\right) } e o local é Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathcal{O}\left(\Delta^{3}\right) } .

Exemplo

Aplicando o algoritmo para o sistema massa-mola visto no método de Euler-Cromer: Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{d^{2}x}{dt^{2}}=-\frac{k}{m}x=-\omega^{2}x }

Podemos ressaltar ainda que Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle a =-\omega^{2}x } e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{dv}{dt}=\frac{d^{2}x}{dt^{2}}} .

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):
  #Primeira etapa
  k1x = fx(v[it])*dt
  k1v = fv(x[it],w2)*dt
  #Segunda etapa
  k2x = fx(v[it]+k1v/2)*dt
  k2v = fv(x[it]+k1x/2,w2)*dt
  #Solução
  x.append(x[it]+k2x)
  v.append(v[it]+k2v)
  #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 ${\textstyle \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:

• ${\textstyle k_{1}=f\left(y_{n},t_{n}\right)\Delta t}$
• ${\textstyle k_{2}=f\left(y_{n}+{\frac {k_{1}}{2}},t_{n}+{\frac {\Delta t}{2}}\right)\Delta t}$
• ${\textstyle k_{3}=f\left(y_{n}+{\frac {k_{2}}{2}},t_{n}+{\frac {\Delta t}{2}}\right)\Delta t}$
• ${\textstyle k_{4}=f\left(y_{n}+k_{3},t_{n}+\Delta t\right)\Delta t}$

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

$$y_{n+1}=y_{n}+{\frac {1}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)$$

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.01 ; 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):
  #Primeira etapa
  k1x = fx(v[it])*dt
  k1v = fv(x[it]     ,w2)*dt
  #Segunda etapa
  k2x = fx(v[it]+k1v/2)*dt
  k2v = fv(x[it]+k1x/2,w2)*dt
  #Terceira etapa
  k3x = fx(v[it]+k2v/2)*dt
  k3v = fv(x[it]+k2x/2,w2)*dt
  #Quarta etapa
  k4x = fx(v[it]+k3v)*dt
  k4v = fv(x[it]+k3x,w2)*dt
  #Solução:
  x.append(x[it]+(k1x+2*k2x+2*k3x+k4x)/6)
  v.append(v[it]+(k1v+2*k2v+2*k3v+k4v)/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)

Ainda podemos chamar a atenção para o fato de que devemos intercalar os coeficientes $k_{i}$ em ambos os métodos, uma vez que coeficientes seguintes dependem dos valores anteriores.

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)