Runge-Kutta 2ª ordem
No método explícito de euler tínhamos:
![{\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+f\left(t_{n},y_{n}\right)\Delta t\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/77f98ba37bde64a1366f737c1dfdf0218bc35cb5)
Sendo
. Podemos reescrever como:
![{\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+ak_{1}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/24e70dca45042acafa3c3c4405f8ef440aa209e2)
Onde
e
. Agor se supormos uma solução:
![{\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+ak_{1}+bk_{2}\qquad \left(1\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/56257a353eafaedc11c14b0e3195e5154dbfbd31)
Com o termo adicional dependendo de uma posição genérica
![{\textstyle y_{n}+cf\left(t_{n},y_{n}\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f640e472e699c92c7f3eee414b9bc4cd6f54c68c)
em um tempo genérico
![{\textstyle t+d\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/257ac0f513eae076d1e55f415e0cb116d5f70cfc)
, 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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b374e2dcd5862d9589dd0558089748d3f6a141e6)
. Usando o fato de que
![{\textstyle k_{1}=f\left(t_{n},y_{n}\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/07d9159df32ead44ca7d1cb11a3c10ccd743f574)
, podemos escrever então que:
![{\textstyle k_{1}=f\left(t_{n},y_{n}\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/07d9159df32ead44ca7d1cb11a3c10ccd743f574)
![{\textstyle k_{2}=f\left(t_{n}+d\Delta t,y_{n}+ck_{1}\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9e56577613dc521be16ff5f1d1f093551e241076)
Agora lembrando a expansão em série de taylor que também vimos no método explícito e Euler:
![{\displaystyle 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!}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/67bc92aa31b8ed764dbe178fb33947a5feaef431)
Abrindo a segunda derivada, temos:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8df1fea6269c53ea2cccb61cbca019a2eee077b3)
Substituindo então, e escrevendo apenas
, temos a seguinte expansão em série de Taylor:
![{\displaystyle 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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3888dd0bbf18c7e31cdf066ee7c3d0d5319e74a7)
Vamos expandir
. Uma expansão de Taylor de primeira ordem para uma função de 2 variáveis em torno de
é dado por [1]:
![{\displaystyle 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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9b7f7b7bbe4b662d57a3f3608075dd7a05e58d20)
Onde
denota a derivada da função
na variável
. Para o nosso caso, temos então para uma expansão em torno de
:
![{\displaystyle 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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/08bfa5f7f7107445fed8dfc6361a549c568d715d)
Expandindo então
em torno de
temos:
![{\displaystyle 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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/10bb8199985f29badd44796ca7e63524a9566294)
Aqui podemos notar que
multiplica a expansão da função, então quando desprezamos os termos de segunda ordem da expansão de
, deprezamos os termos de terceira ordem de
. Substituindo então o
aproximado e
na equação 1, temos:
![{\displaystyle 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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/45b93c9c39adcf56eb02baa3be698e0b5b9e71f8)
Manipulando:
![{\displaystyle 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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5f54c94ed040c30dcf26e8b2c364bde41a533715)
Comparando a aproximação 3 com a expansão 2 temos a seguinte relação:
![{\displaystyle {\begin{aligned}a+b&=1\\bd&={\frac {1}{2}}\\bc&={\frac {1}{2}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e1c935e0f5ebfbeed7df635f539f28241238d9c1)
Diferentes conjuntos de valore satisfazem este sistema. O método do ponto médio é obtido se ecolhermos:
![{\textstyle d=c={\frac {1}{2}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/98f4f67c4e258565a66cb7d04691adbbf2891806)
,
![{\textstyle b=1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2049de5a23808e51215bbf8e0ec13fef73bdbb27)
e
![{\textstyle a=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/52fab3c0f8991d8efe516667862df39e89785f63)
:
![{\textstyle k_{2}=f\left(t_{n}+{\frac {\Delta t}{2}},y_{n}+{\frac {k_{1}}{2}}\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e15e1b0a0cd5f083c76136a2539cb1171bb283f9)
Então:
![{\displaystyle {\begin{aligned}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{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0847716a1d883b499def970a6dff21427ad6e12d)
O método de Heun é obtido se for escolhido
e
:
![{\textstyle k_{1}=f\left(t_{n},y_{n}\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/07d9159df32ead44ca7d1cb11a3c10ccd743f574)
![{\textstyle k_{2}=f\left(t_{n}+\Delta t,y_{n}+k_{1}\Delta t\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/98beedbddf656ae0fd0e50a63deeff652c611529)
![{\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+\left(k_{1}+k_{2}\right){\frac {1}{2}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/599e9019279e06225f14403a2f75ea5273c619e2)
Uma observação, é que o erro global no algoritmo de Runge-Kutta de segunda ordem é
e o local é
.
Exemplo
Aplicando o algoritmo para o sistema massa-mola visto no método de Euler-Cromer:
![{\displaystyle {\frac {d^{2}x}{dt^{2}}}=-{\frac {k}{m}}x=-\omega ^{2}x}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e1d1fabd87d796e63c17f0a8d1c31ac438b82ef4)
Podemos ressaltar ainda que
e
.
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
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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/010a1dffdd1a9cd7c1aff2fd847dc397211ef16f)
![{\textstyle k_{2}=f\left(y_{n}+{\frac {k_{1}}{2}},t_{n}+{\frac {\Delta t}{2}}\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e7231f05262df80d27800b4241181552d717b95a)
![{\textstyle k_{3}=f\left(y_{n}+{\frac {k_{2}}{2}},t_{n}+{\frac {\Delta t}{2}}\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/373b7d7d49260ae8280d3fb9bdc59fc19cf2c39d)
![{\textstyle k_{4}=f\left(y_{n}+k_{3},t_{n}+\Delta t\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/486e7870a94c38c11aa2ed40bd458a282fadaac9)
E por fim, temos então que o novo valor será dado por:
![{\displaystyle y_{n+1}=y_{n}+{\frac {1}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/baea1e2c7da18d3948e96e4f1230c138655a61ad)
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
em ambos os métodos, uma vez que coeficientes seguintes dependem dos valores anteriores.
Principais materiais utilizados
- Métodos de Runge-Kutta explícitos (REAMAT, UFRGS)
- Runge-Kutta Methods (Michael Zeltkevic, Instituto de Tecnologia de Massachusetts)
- Second Order Runge-Kutta (Erik Cheever, Swarthmore)
Citações