Introdução à dinâmica molecular clássica

só um divisor

Resolução de equações diferencial de II grau

Método de Runge Kutta

A ideia básica do método de Euler é aproximar a derivada pelo método das diferencias finitas

\[ \frac{df}{dt} = \frac{f_{i+1} - f_i}{t_{i+1} - t_i} \nonumber \]

de onde

\[ f^E_{i+1} = \frac{df}{dt}\left(t_{i+1} - t_i\right) + f_i \nonumber \]

Na figura acima comparamos a aproximação feita pelo método de Euler com a aproximação feita com um outro método conhecido como méetodo de Runge-Kutta. Observe que a solução dada pelo método de Euler, $f_{i+1}^E$, está bem mais afastada do ponto certo, $f_{i+1}$, do que aquela obtida pelo método de Runge-Kutta, $f_{i+1}^{R}$, por tanto parece que este método é bastante superior ao método de Euler.

Método de Runge-Kutta de II ordem

A técnica de Runge-Kutta permite varias formas diferentes de equações de integração já que é um algoritmo que não se converte diretamente na equação, a diferencia do Euler. Na figura acima foi utilizado o método de Runge-Kutta de II ordem que a continuação passo a descrever:

  1. calcule com o método de Euler o valor da função no ponto $t_{i+\frac{1}{2}}$, isto é, $t+\frac{\Delta t}{2}$, na figura chamamos esse ponto de $f^E_{i+\frac{1}{2}}$.
  2. para calcular o ponto $f_{i+1}^{R}$, aplique o método de de Euler com a diferencia de que a derivada não é calculada no ponto inicial e sim no ponto previamente calculado pelo passo anterior, $f^E_{i+\frac{1}{2}}$.
Agora da para entender o porque na figura acima a inclinação das curvas é diferente, a curva em vermelha é a reta calculada utilizando a derivada no ponto inicial, com essa curva podemos estimar a posição do em $t+\frac{\Delta t}{2}$; a curva em azul tem a inclinação dada pela derivada no ponto calculado previamente em $t+\frac{\Delta t}{2}$, mas a reta passa pelo ponto em $t$ e não em $t+\frac{\Delta t}{2}$

Matematicamente podemos expressar o acima dito como: \begin{eqnarray*} f_{i+\frac{1}{2}} =& f_i + \frac{\Delta t}{2} &\left. \frac{df}{dt} \right|_{\,i} \\ f_{i+1} =& f_i + \Delta t &\left.\frac{df}{dt} \right|_{\,i+\frac{1}{2}} \end{eqnarray*}

Tradicionalmente o método de Runge-Kutta de II ordem define duas contantes: \begin{align*} k_1 =& \Delta t\; Df\left(t_i, \,f_i\right)\\ k_2 =& \Delta t\; Df\left(t_i+\frac{\Delta t}{2}, \,f_i+\frac{1}{2}k_1\right) \end{align*} onde $Df$ é a derivada da função avaliada no argumento (veja que pode depender de f), e a partir da última constante calculamos o novo valor: \[ f_{i+1} = f_i + k_2 \nonumber \]

Método de Runge-Kutta de IV ordem

Como vimos a técnica de Runge-Kutta consiste ir nos aproximando do ponto onde queremos calcular a derivada de forma sucessiva, utilizando para isso o valor da derivada no ponto no meio do caminho. Na verdade podemos aplicar essa ideia quantas vezes queiramos, utilizando o ponto médio, o ponto a 3/4, o ponto 7/8, o ponto a 15/16, o ponto a 31/32, etc. Em especial o ponto a 7/8 é muito utilizado e se conhece como método de Runge-Kutta de IV ordem.

A constantes de Runge-Kutta no método de IV ordem são: \begin{align*} k_1 = & \Delta t\; Df\left(t_i, \,f_i\right)\\ k_2 = & \Delta t\; Df\left(t_i+\frac{\Delta t}{2}, \,f_i+\frac{1}{2}k_1\right)\\ k_3 = & \Delta t\; Df\left(t_i+\frac{\Delta t}{2}, \,f_i+\frac{1}{2}k_2\right)\\ k_4 = & \Delta t\; Df\left(t_i, \,f_i+k_3\right) \end{align*} de onde \[ f_{i+1} = f_i + \frac{1}{6}\left( k_1 + 2\,k_2 + 2\,k_3 + k_4 \right) \nonumber \]

Veremos no exemplo a seguir que o Runge-Kutta de IV ordem da um resultado muito melhor do que o Euler-Crome, contudo ele é pouco utilizado na dinâmica molecular clássica devido a que são necessários vários passos para realizar a integração.

No exemplo a seguir calculamos o movimento de um oscilador harmônico utilizando o Runge-Kutta (solução x, v), Euler-Crome (solução xec, vec) e os comparamos à solução exata (solução xe, ve):

      #/usr/bin/env python
      #-*- coding: utf-8 -*-
      
      import numpy as np
      import matplotlib as mpl
      import matplotlib.pyplot as plt
      from mpl_toolkits.axes_grid1.inset_locator \
        import inset_axes, zoomed_inset_axes

      def Dv(w2, xo, x, v):
        return -w2 * (x - xo)
        
      def Dx(x, v):
        return v

      #******************************

      m  = 1.0
      k  = 10.0
      xo = 1.0
      xi = 1.1
      vi = 0.0
      ti = 0.0
      tf = 10.0
      dt = 1.0e-1

      w2 = 10.0 / 1.0
      N  = int ( (tf - ti) / dt )


      #Runge-Kutta de IV ordem
      x = np.zeros( N+1, dtype=float)
      v = np.zeros_like(x)
      t = np.zeros_like(x)

      x[0] = xi
      v[0] = vi
      t[0] = ti

      for i in range(N):
        
        kv1 = dt * Dv(w2, xo, x[i], v[i])
        kx1 = dt * Dx(x[i], v[i])
        
        kv2 = dt * Dv(w2, xo, x[i]+0.5*kx1, v[i]+0.5*kv1)
        kx2 = dt * Dx(x[i]+0.5*kx1, v[i]+0.5*kv1)
        
        kv3 = dt * Dv(w2, xo, x[i]+0.5*kx2, v[i]+0.5*kv2)
        kx3 = dt * Dx(x[i]+0.5*kx2, v[i]+0.5*kv2)
        
        kv4 = dt * Dv(w2, xo, x[i]+kx3, v[i]+kv3)
        kx4 = dt * Dx(x[i]+kx3, v[i]+kv3)

        t[i+1] = t[i] + dt
        x[i+1] = x[i] + 1.0 / 6.0 * (kx1 + 2.0 * kx2 + 2.0 * kx3 + kx4)
        v[i+1] = v[i] + 1.0 / 6.0 * (kv1 + 2.0 * kv2 + 2.0 * kv3 + kv4)
        

      #Euler-Crome
      xec = np.zeros( N+1, dtype=float)
      vec = np.zeros_like(xec)

      xec[0] = xi
      vec[0] = vi

      for i in range(N):
        vec[i+1] = vec[i] + Dv(w2, xo, xec[i], vec[i]) *dt
        xec[i+1] = xec[i] + Dx(xec[i], vec[i+1])*dt

      w = np.sqrt(w2)

      dx   = xi - xo
      xmax = dx
      vmax = -w * dx
        
      #Solução exata
      xe = np.zeros( N+1, dtype=float)
      ve = np.zeros_like(xe)

      xe = xo + xmax * np.cos(w * t)
      ve = vmax * np.sin(w * t)

      #a figura
      fig     = plt.figure(figsize=(10,10))  
      
      #cria uma grade 1x1 e coloca a figura na cela 1
      subfig1 = fig.add_subplot(111) 

      #coloca um inset que inicia em (0.88, 0.88)
      inset   = inset_axes(subfig1, width=2, height=2, \
                          bbox_to_anchor=(0.88, 0.88), \
                          bbox_transform=subfig1.figure.transFigure)

      inset.set_xlim( 0.95, 1.02 )
      inset.xaxis.set_ticks( np.linspace(0.95, 1.02, 3, endpoint=True) )
      inset.set_ylim( -0.321, -0.300 )
      inset.yaxis.set_ticks( np.linspace(-0.321, -0.300, 3, endpoint=True) )

      subfig1.plot(xe, ve, '-', label='exata')
      subfig1.plot(x, v, '--', label='Runge-Kutta')
      subfig1.plot(xec, vec, ':', label='Euler-Crome')

      subfig1.set_xlabel(u'Posição', fontsize=20)
      subfig1.set_ylabel('Velocidade', fontsize=20)

      inset.plot(xe, ve, '-')
      inset.plot(x, v, '--')
      inset.plot(xec, vec, ':')

      subfig1.legend(loc=4)

      plt.savefig( 'fig02.png', format='png' )
      plt.show()
      

Além do exemplo mostrar como utilizar o Runge-Kutta de IV ordem, se mostra a colocação de um inset (ver o exemplo do site do matplotlib para situações mais complexas) onde se faz um zoom para ver que nem mesmo o Runge-Kutta é exatamente igual à solução exata, finalmente se manda a imprimir em um arquivo