|
|
Linha 79: |
Linha 79: |
|
| |
|
| <pre> | | <pre> |
| | 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) |
| </pre> | | </pre> |
|
| |
|
Edição das 17h48min de 15 de março de 2022
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=-1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b9e85672bf966e24b4c7eb1a3f251b21716f1d6c)
:
![{\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
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
Exemplo
Principais materiais utilizados
- Runge-Kutta Methods (Michael Zeltkevic, Instituto de Tecnologia de Massachusetts)
- Second Order Runge-Kutta (Erik Cheever, Swarthmore)
Citações