Física Computacional - FSC-5705

só um divisor

Diferenciação numérica

A ideia por trás da derivação numérica consiste na discretização da derivada. Dessa forma devemos pensar em quais pontos queremos calcular nossa derivada. Está técnica se conhece como método das diferencias finitas já que diferentemente do que se faz no cálculo da derivada de uma função não consideramos incrementos infinitamente pequenos e sim finitos ou dito de outro modo o limite indo para zero é desconsiderado.

Do anterior vemos que é necessário que a função para a qual pretendemos calcular a derivada, $f(x)$ esteja definida nos pontos onde será feito o calculo, suponhamos que esses pontos sejam pontos definidos dentro do intervalo $[a,b]$, isto é, ${x_0=a,\,x_1,\,x_2,\,\ldots, x_i,\,\ldots, x_{N+1}=b}$. Por comodidade vamos supor que para todos os pares contíguos de pontos se verifica a seguinte condição $x_{i+1}-x_i=\Delta x$, ou dito de outro modo, qualquer ponto desse conjunto pode ser obtido de \[ x_i = a + i\Delta x\;\;\;\;\;\;\; i=0,\,1,\,2,\,\ldots,\,N+1 \nonumber \] onde \[ \Delta x = \dfrac{b - a}{N+1} \nonumber \] com essa escolha de pontos definimos a derivada da numérica da função $f(x)$ como \[ f'(x_i) = \dfrac{f\left(x_i+\Delta x\right) - f\left(x_i\right)}{\Delta x} \label{doispontos} \] ou, utilizando a notada da serie de pontos previamente escolhidos: \[ f'(x_i) = \dfrac{f\left(x_{i+1}\right) - f\left(x_i\right)}{\Delta x} + O\left(\Delta x \right) \label{adiantada} \]

derivada numerica

Na figura acima se mostra a diferencia existente entre a derivada numérica e a deriva exata da função $f(x)$ calculada no ponto $x_0$. Pode se mostrar, utilizando os polinomios de Lagrange, que o erro, $O\left(\Delta x \right)$, de se utilizar essa equação está dado por \[ O\left(\Delta x \right) = \dfrac{\Delta x}{2} f''(\xi) \nonumber \] onde $\xi$ é algum número entre $x_{i+1}$ e $x_i$.

A definição $\ref{doispontos}$ se conhece como derivada em 2 pontos e ela nos permite definir a deriva numérica utilizando 2 pontos contíguos da serie. Observe que podemos escolher uma outra definição equivalente à definição de derivada adiantada dada pela equação $\ref{adiantada}$, se colocamos $\Delta x \gets -\Delta x$ obtemos a derivada atrasada dada por \[ f'(x_i) = \dfrac{f\left(x_{i}\right) - f\left(x_{i-1}\right)}{\Delta x} + O\left(\Delta x \right) \label{atrasada} \]

Podemos aumentar a precisão se fizermos uma derivada centrada: \[ f'(x_i) = \dfrac{f\left(x_{i+1}\right) - f\left(x_{i-1}\right)}{2\Delta x} + O\left(\Delta t^2\right) \label{centrado} \] com \[ O\left(\Delta t^2\right) = -\dfrac{\Delta t^2}{6} f'''(\xi) \nonumber \] onde $\xi$ é um número entre $x_{i+1}$ e $x_{i-1}$

derivada numerica 2

Derivadas avaliadas com mais de 2 pontos

Derivada de 3 pontos à direita \[ f'(x_i) = \dfrac{-3f\left(x_{i}\right) + 4f\left(x_{i+1}\right) - f\left(x_{i+2}\right)}{2\Delta x} + \dfrac{\Delta t^2}{3} f'''(\xi) \nonumber \]

Derivada de 5 pontos à direita \[ f'(x_i) = \dfrac{-25f\left(x_{i}\right) + 48f\left(x_{i+1}\right) - 36f\left(x_{i+2}\right) + 16f\left(x_{i+3}\right)-3f\left(x_{i+4}\right)}{12\Delta x} \;\;\;\;+ \dfrac{\Delta t^4}{5} f^{(4)}(\xi) \nonumber \]

Derivada de 5 pontos centrada \[ f'(x_i) = \dfrac{f\left(x_{i-2}\right) - 8f\left(x_{i-1}\right) + 8f\left(x_{i+1}\right) - f\left(x_{i+2}\right)}{12\Delta x} \;\;\;\;+ \dfrac{\Delta t^4}{30} f^{(4)}(\xi) \nonumber \]

Derivada de II ordem

Combinando duas definições de derivada de I ordem obtemos a definição de deriva de II ordem \[ f''(x_i) = \dfrac{f\left(x_{i-1}\right) - 2 f\left(x_{i}\right) + f\left(x_{i+1}\right)}{\Delta x^2} - \dfrac{\Delta t^2}{12} f^{(4)}(\xi) \]

Exercícios

  1. Utilize as equações de derivada adianta e atrasada para calcular as derivadas das funções (a) $\sin(x)$, (b) $e^x - 2x + 3x - 1$. Calcule o erro e compare com o valor exato os resultados obtidos. Utilize $\Delta x = 0.1$
  2. Dada as funções: (a) $e^{2x}$, (b) $x\ln x$, (c) $x\cos(x) - x^2\sin(x)$ e (d) $2\left( \ln x\right)^2 + 3 \sin x$ calcule a derivada centrada. Calcule o erro e compare com o valor exato. Utilize $\Delta x = 0.1$
  3. Determine a derivada numérica tão precisa quanto puder das funções (a) $f(x) = \tan (x)$ e (b) $f(x) = e^{x/3}+ x^2$. Calcule o erro e compare com o valor exato. Utilize $\Delta x = 0.1$
  4. Seja $f(x)=3xe^x - \cos x$. Calcule $f''(1.3)$ utilizando $\Delta x = 0.1$ e $\Delta x = 0.01$. Compare os resultados com o valor exato.
  5. Num circuito $LR$ a tensão é função do tempo e está dada pela equação \[ \epsilon (t) = L \dfrac{di}{dt} + Ri \nonumber \] onde $R$ é a resistência e $i$ a corrente que passa pelo circuito. Suponha que a corrente foi medida em varios valore s de $t$ e os resultados obtidos sejam \[ \begin{array}{c|c|c|c|c|c} t & 1.00 & 1.01 & 1.02 & 1.03 & 1.04\\\hline i & 3.1 & 3.12 & 3.14 & 3.18 & 3.42 \nonumber \end{array} \] onde $t$ está em segundo e $i$ em ampères, a indutância $L$ é um valor constante de $0,98$ henry e a resistência é $0,142\Omega$. Aproxime o valor de $\epsilon (t)$ para os tempos dado na tabela.
  6. Considere a função \[ e(h) = \dfrac{\epsilon}{h} + \dfrac{h^2}{6}M \nonumber \] onde $\epsilon$ e $M$ são constantes. Mostre que $e(h)$ tem um mínimo em $\sqrt[3]{3\epsilon / M}$

Solução da questão 6

A ideia era solucionar a questão de forma totalmente numérica. O único calculo analítico era na avaliação de quais pontos seriam utilizados.

      #!/usr/bin/env python
      # -*- coding: utf-8 -*-

      import numpy as np
      import matplotlib.pyplot as plt

      #---------- Função externa ----------
      def f(e, m, h):
        return (e/h + h**2/6.0*m)

      def df(e, m, h):
        return(-e/h**2 + h/3.0 * m)
        
      #---------- interpolação de Lagrange ----------

      def Lnk(k, x, xi):
        l_nk = 1.0
        n   = len(xi)
        for i in xrange(n):
          if (i != k):
            l_nk = l_nk * ( (x - xi[i]) / (xi[k] - xi[i]) )
        
        return l_nk

      def lagrange(yi, xi, x):
        n  = len(xi)
        P = 0.0
        for k in range(n):
          P = P + yi[k]*Lnk(k, x, xi)
        return P

      #---------- calculo de raiz - bisseção ----------
      def bissecao(xa, xb, ya, yb, x, y, tol):
        a = xa
        b = xb
        iteracao = 0
        while True:
          c = a + (b - a) / 2.0
          fa = lagrange(y, x, a)
          fc = lagrange(y, x, c)
          if (a == xa):
            fa = ya
          elif (a == xb):
            fa = yb

          if (fa * fc < 0.0):
            a = a
            b = c
          else:
            a = c
            b = b
          
          if (a == 0.0):
            a = 1.0e-5
          elif (b == 0.0):
            b = 1.0e-5
            
          if (np.abs((b - a) / a) < tol):
            break
          
          iteracao = iteracao + 1
        return (a)

      #---------- calculo da derivada numérica ----------

      print
      print 'Com este programa vamos tentar determinar o mínimo'
      print 'da função dada no problema 6 de forma numérica'
      print 'utilizando a derivada centrada de 5 pontos'
      print '---------------------------------------------------'
      print
      print 'vamos definir os parâmetros da função a ser derivada'

      e = input('valor de epsilon: ')
      m = input('valor de H: ')

      print 'primeiro devemos determinar o intervalo adequado'
      print 'a função derivada será calculada com 100 pontos'

      while True:
        a = input('valor inferior do intervalo: ')
        b = input('valor superior do intervalo: ')
        x = np.linspace(a, b, 104)
        y = f(e, m, x)
        plt.plot(x, y, '-')
        plt.show()
        ok = input('é adequado o intervalo 1-sim, 0-não: ')
        if (ok == 1):
          break
          
      dy = np.zeros(100)
      xn = np.zeros(100)
      j  = 0
      dx = x[2] - x[1]
      for i in range(2, len(x)-2):
        dy[j] = (f(e, m, x[i-2]) - 8.0*f(e, m, x[i-1]) +
                8.0*f(e, m, x[i+2]) - f(e, m, x[i+2])) / (12.0 * dx)
        xn[j] = x[i]
        j     = j + 1
        
      plt.plot(xn, dy, '-')
      plt.show()

      tol = 1.0e-14

      print 'a interpolação será feita entre', x[0], 'e ', x[-1]
      print 'o erro estimado será de: ', tol

      minimo = bissecao(xn[10], xn[-10], dy[10], dy[10], xn, dy, tol)
      print
      print 'O minimo aproximado é: ', minimo
      print 'O minimo calculado é:  ', (3.0*e/m)**(1.0/3.0)
      
[usuario@pclabfis: ]# python deriva_p6.py

Com este programa vamos tentar determinar o mínimo
da função dada no problema 6 de forma numérica
utilizando a derivada centrada de 5 pontos
---------------------------------------------------

vamos definir os parâmetros da função a ser derivada
valor de epsilon: 25.0
valor de H: 3.0
primeiro devemos determinar o intervalo adequado
a função derivada será calculada com 100 pontos
valor inferior do intervalo: -10.0
valor superior do intervalo: 10.0
Integração
é adequado o intervalo 1-sim, 0-não: 0
valor inferior do intervalo: 0.1
valor superior do intervalo: 10.0
Integração
é adequado o intervalo 1-sim, 0-não: 0
valor inferior do intervalo: 0.5
valor superior do intervalo: 10.0
Integração
é adequado o intervalo 1-sim, 0-não: 1
a interpolação será feita entre 0.5 e 10.0
Integração
o erro estimado será de: 1e-14

O minimo aproximado é: 2.87031214612
O minimo calculado é: 2.92401773821