Classes para derivadas

Relembrando a classe para a derivação numérica

$$f'(x) \approx {f(x+h)-f(x) \over h}$$

class Derivada:
    def __init__(self, f, h=1E-5):
        self.f = f
        self.h = float(h)

    def __call__(self, x):
        f, h = self.f, self.h      ## formando formas curtas
        return (f(x+h) - f(x))/h

def f(x):
    return exp(-x)*cos(tanh(x))

from math import exp, cos, tanh
dfdx = Derivada(f)
print dfdx(2.0)

Existem diversas fórmulas para a derivação numérica

$$\begin{align*} f'(x) &= \frac{f(x+h)-f(x)}{h} + { \cal O}(h) \quad \text{diferencial a frente de 1 ordem} \\ f'(x) &= \frac{f(x)-f(x-h)}{h} + { \cal O}(h) \quad \text{diferencial à trás de 1 ordem} \\ f'(x) &= \frac{f(x+h)-f(x-h)}{2h} + { \cal O}(h^2) \quad \text{diferencial central de 2 ordem} \\ f'(x) &= \frac{4}{3} \frac{f(x+h)-f(x-h)}{2h} - \frac{1}{3} \frac{f(x+2h) - f(x-2h)}{4h} \\ &+ { \cal O}(h^4) \quad \text{diferencial central de 4 ordem} \\ f'(x) &= \frac{3}{2} \frac{f(x+h)-f(x-h)}{2h} - \frac{3}{5} \frac{f(x+2h) - f(x-2h)}{4h} \\ &+ \frac{1}{10} \frac{f(x+3h) - f(x-3h)}{6h} + { \cal O}(h^6) \quad \text{diferencial central de 6 ordem} \\ f'(x) &= \frac{1}{h} \left( - \frac{1}{6}f(x+2h) + f(x+h) - \frac{1}{2}f(x) - \frac{1}{3}f(x-h) \right) \\ &+ { \cal O}(h^3) \quad \text{diferencial a frente de 3 ordem} \end{align*}$$

Como podemos escrever um módulo que oferece todas essas fórmulas?

É fácil:

class Frente1:
    def __init__(self, f, h=1E-5):
        self.f = f
        self.h = float(h)

    def __call__(self, x):
        f, h = self.f, self.h
        return (f(x+h) - f(x))/h

class Tras1:
    def __init__(self, f, h=1E-5):
        self.f = f
        self.h = float(h)

    def __call__(self, x):
        f, h = self.f, self.h
        return (f(x) - f(x-h))/h

class Central2:
    ## mesmo contrutor
    ## coloque as fórmulas relevantes em __call__

Qual é o problema com esse tipo de código?

Todos os construtores são iguais então muito código é duplicado.

• Ideia geral: colocar o código em comum em uma superclasse;
• Herdar o construtor da superclasse fazendo subclasses para as diferentes fórmulas de derivação utilizando o __call__

Classe hierárquica para derivação numérica

Superclasse

class derivada:
    def __init__(self, f, h=1E-5):
        self.f = f
        self.h = float(h)

Subclasse para a equação de derivação de primeira ordem para frente

class frente1(derivada):
    def __call__(self, x):
        f, h = self.f, self.h
        return (f(x+h) - f(x))/h

Subclasse para a equação de derivação de quarta ordem central

class central4(derivada):
    def __call__(self, x):
        f, h = self.f, self.h
        return (4./3)*(f(x+h)   - f(x-h))  /(2*h) - \
               (1./3)*(f(x+2*h) - f(x-2*h))/(4*h)

Uso das classes de derivação

Exemplo: sendo $ f(x) = sen(x) $, calcule $ f'(x) $ para $ x = \pi $, usando a derivação central de ordem 4.

>>> from derivada import *
>>> from math import sin
>>> meucos = central4(sin)
>>> # calcule sin'(pi):
>>> meucos(pi)
-1.000000082740371

central4(sin) chama o construtor herdado na superclasse, enquanto meucos(pi) chama __call__ na subclasse central4

Exercício 2: avaliando o fluxo do programa

Um programa flexível para a derivação numérica

Suponha que nós queremos derivar equações da linha de comando:

Terminal> python df.py 'exp(sin(x))' Central 2 3.1
-1.04155573055

Terminal> python df.py 'f(x)' difftype difforder x
f'(x)

Usando o eval e a classe hierárquica derivada este programa principal pode ser feito em poucas linhas (muitas linhas em C# e Java):

import sys
from derivada import *
from math import *
from scitools.StringFunction import StringFunction

f = StringFunction(sys.argv[1])
derivada_tipo = sys.argv[2]
derivada_ordem = sys.argv[3]
nome_classe = derivada_tipo + derivada_ordem
df = eval(nome_classe + '(f)')
x = float(sys.argv[4])
print df(x)

Investigando erros numéricos por aproximações

• Podemos investigar, empiricamente, a precisão de nossa família de 6 equações de derivação numérica;
• Função de teste: $ f(x) = exp(-10x) $;
• Queremos o resultado na seguinte forma:

   h         Frente1        Central2        Central4
6.25E-02 -2.56418286E+00  6.63876231E-01 -5.32825724E-02
3.12E-02 -1.41170013E+00  1.63556996E-01 -3.21608292E-03
1.56E-02 -7.42100948E-01  4.07398036E-02 -1.99260429E-04
7.81E-03 -3.80648092E-01  1.01756309E-02 -1.24266603E-05
3.91E-03 -1.92794011E-01  2.54332554E-03 -7.76243120E-07
1.95E-03 -9.70235594E-02  6.35795004E-04 -4.85085874E-08

Observações:

• Reduzindo à metade $h$ de linha a linha o erro se reduz por um fator de 2, 4 e 16, isto é, os erros se reduzem em $h$, $h^2$ e $h^4$;
• central4 tem uma precisão muito superior se comparada com frente1.

Implementação

class derivada2(derivada):
    def __init__(self, f, h=1E-5, dfdx_exact=None):
        derivada.__init__(self, f, h)
        self.exact = dfdx_exact
    def error(self, x):
        if self.exact is not None:
            df_numerical = self(x)
            df_exact = self.exact(x)
            return df_exact - df_numerical

class frente1(derivada2):
    def __call__(self, x):
        f, h = self.f, self.h
        return (f(x+h) - f(x))/h

class central2(derivada2):
    def __call__(self, x):
        f, h = self.f, self.h
        return (f(x+h) - f(x-h))/(2*h)

class central4(derivada2):
    def __call__(self, x):
        f, h = self.f, self.h
        return (4./3)*(f(x+h) - f(x-h)) /(2*h) - (1./3)*(f(x+2*h) - f(x-2*h))/(4*h)

Uso

from derivada2 import *
from math import exp

def f1(x):
    return exp(-10*x)
def df1dx(x):
    return -10*exp(-10*x)
def tabela(f, x, h_values, metodos, dfdx=None):
    # imprimir o cabeçalho (h e os nomes das classes para os métodos)
    print ’    h    ’,
    for metodo in metodos:
        print ’%-15s’ % method.__name__,
    print # nova linha
    # imprimir a tabela
    for h in h_values:
        print ’%10.2E’ % h,
    for metodo in metodos:
        if dfdx is not None:
        # escrever o erro
        d = metodo(f, h, dfdx)
        output = d.error(x)
        else:
        # escrever o valor
        d = metodo(f, h)
        output = d(x)
        print ’%15.8E’ % output,
    print # nova linha


tabela(f1, 0, [2**(-k) for k in range(10)], [frente1, central2, central4], df1dx)