Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como fazer uma classe em Python para diferenciação numérica?

0 votos
31 visitas
perguntada Ago 6 em Ciência da Computação por Pedro Watuhã (1 ponto)  
Compartilhe

1 Resposta

0 votos
respondida Ago 6 por Pedro Watuhã (1 ponto)  
editado Ago 12 por Pedro Watuhã
 
Melhor resposta

Boa Tarde a todos,
A seguinte questão foi enunciada no exercício 7.15 do livro A primer on Scientific Programming with Python First Edition.

Uma formula utilizada amplamente para a diferenciação numérica de uma função \(f(x)\) é dada por: \[f\'(x) \approx \frac{f(x+h) - f(x-h)}{2h} \]
Esta formula normalmente fornece-nos derivadas mais aproximadas, pois aplica uma diferença central em fez de unidirecional.
Inicialmente, definiremos uma classe Central responsável pelo processo de diferenciação e utilizaremos \(10^{-9}\) como um valor para \(h\) em uma primeira estimativa. Em seguida, aplicaremos essa classe sob a função \(0.25x^4\) e compararemos os resultados àqueles obtidos pela derivação da função, \(x^3\), nos pontos \(x = \{1,5,10\}\).

class Central():
def __init__(self, f):
    self.f = f
def __call__(self,x):
    self.x=x
    return (f(x+10**(-9))-f(x-10**(-9)))/(2*10**(-9))


def f(x):
    return 0.25*x**4
if __name__=="__main__":
    df = Central(f)
    for x in (1, 5, 10):
        df_value = df(x)
        exact = x**3
        print ("f’(%d)=%g (error=%.2E)" % (x, df_value, exact-df_value))

\[ f’(1)=1 (error=-2.72E-08) \]
\[ f’(5)=125 (error=-1.39E-05) \]
\[ f’(10)=1000 (error=1.16E-04) \]

Com esse resultado, observamos o quão precisa pode ser a estimativa do modelo utilizando um \(h\) pequeno. Então, testemos o quão preciso o modelo pode ser em estimar \(log(x)\) utilizando diferentes valores de \(h\) no ponto \(x=10\). Para isso, alteraremos a classe Central para incluir as possíveis alterações em \(h\).

import numpy as np
class Central():
    def __init__(self, f):
        self.f = f
    def __call__(self,x,h):
        self.x = x
        self.h = h
        return (f(x+h)-f(x-h))/(2*h)
def f(x):
    return np.log(x)

if __name__=="__main__":
    df = Central(f)
    x=10
    for h in (0.5, 0.1, 10**(-3),10**(-5),10**(-7),10**(-9),10**(-11)):
        df_value = df(x,h)
        exact = 1/x
        print ("h = "+ str(h), "f’(%d)=%g (error=%.2E)" % (x, df_value, exact-df_value))

\[ h = 0.5: f’(10)=0.100083 (error=-8.35E-05) \]
\[ h = 0.1: f’(10)=0.100003 (error=-3.33E-06) \]
\[ h = 0.001: f’(10)=0.1 (error=-3.33E-10) \]
\[ h = 1e-05: f’(10)=0.1 (error=8.23E-12) \]
\[ h = 1e-07: f’(10)=0.1 (error=6.08E-10) \]
\[ h = 1e-09: f’(10)=0.1 (error=-8.27E-09) \]
\[ h = 1e-11: f’(10)=0.0999867 (error=1.33E-05) \]

Nesse caso, percebe-se que, dentre os \(h\) utilizados, \(10^{-5}\) fornece-nos o menor erro e a utilização de valores à direita ou à esquerda tendem a aumentá-lo. Podemos observar essa tendência nos seguintes gráficos em que o intervalo é seccionado de formas diferentes :

import numpy as np
import matplotlib.pyplot as plt
def f(x):
    return np.log(x)

if __name__=="__main__":
    error=[]
    df = Central(f)
    x=10
    for h in np.arange(0,0.5,10**(-6)):
        df_value = df(x,h)
        exact = 1/x
        error.append(exact-df_value)
    x=np.arange(0,0.5,10**(-6))
    y=error
    plt.xlabel('h')
    plt.ylabel('Erro')
    plt.plot(x,y)

A imagem será apresentada aqui.

import numpy as np
import matplotlib.pyplot as plt
def f(x):
    return np.log(x)

if __name__=="__main__":
    error=[]
    df = Central(f)
    x=10
    for h in np.arange(0,0.001,10**(-10)):
        df_value = df(x,h)
        exact = 1/x
        error.append(exact-df_value)
    x=np.arange(0,0.001,10**(-10))
    y=error
    plt.plot(x,y)
    plt.xlabel('h')
    plt.ylabel('Erro')

A imagem será apresentada aqui.

Por conseguinte, observamos que a tendência global é diminuir o erro conforme o h diminui, mas pela forma com que o Python ou o Numpy aparentam lidar com números cada vez menores, o erro tende a aumentar.

...