Derivada numérica

$\newcommand{\fl}{\mbox{fl}}$

Este exemplo é bem instrutivo pois mostra como podemos nos enganar ao tentar reproduzir numericamente um resultado teórico.

Do ponto de vista teórico, é verdade que a diferenças centradas, dada por $$ D_f(x) = \frac{f(x+h)-f(x-h)}{2h}, $$ realmente é uma aproximação de segunda ordem para $f'(x)$, quando $h$ tende a zero, isto é, \begin{equation} \label{limitanteexato} |D_f(x) - f'(x)| \le \frac{M_3h^2}{6}, \end{equation} sempre que $h \le \bar{h}$, para algum $\bar{h}$ e supondo que $f$ tem pelo menos até terceira derivada contínua. Portanto, é correto dizer que a aproximação de diferenças centradas converge para $f'(x)$, quando $h$ tende a zero.

Este resultado coloca até mesmo a pergunta se precisamos conhecer de fato a derivada de uma função. Não seria suficiente avaliar a função duas vezes e usar a fórmula de diferença centrada para computar o valor da derivada? Se quisermos um valor muito preciso, bastaria tomar $h$ muito pequeno.

Mesmo confiando no resultado acima, nossa tarefa como cientistas é desconfiar sempre. Nada mais justo que um experimento numérico para aumentar nossa segurança, não?

O experimento foi montado assim. Consideramos a função $f(x)= \sin(x)$, para simplificar o cálculo de $M_3 = \max|f'''(x)|$. No caso da função seno, $M_3=1$. Escolhemos também um valor específico para $x$ que não fosse nenhum ângulo notável particular. Depois fomos reduzindo gradativamente o valor do passo de discretização $h$ usado no cálculo da diferença centrada e avaliando o erro no valor numérico da derivada, contra o valor exato da derivada, que sabemos ser $\cos(x)$. O valores para $h$ foram tomados como $10^{0},10^{-1},\ldots,10^{-14}$. Os resultados para este experimento foram produzidos com o bloco de código abaixo.

> f = @(x) sin(x);
> g = @(x) cos(x);                      # derivada de f
> dc = @(x,h) (f(x+h)-f(x-h))/(2*h);    # diferença centrada
> x = pi/3.2;                           # valor escolhido para x
> for n = 0:14;                         # loop para diversos valores para h
    h = 1e-n;                           # h = 10^(-n)
    erro(n,1) = abs( g(x) - dc(x,h) );
  endfor
> format short e
> erro
erro =

   8.8074e-02
   9.2549e-04
   9.2595e-06
   9.2595e-08
   9.2653e-10
   1.2262e-11
   2.1045e-11
   4.2304e-10
   6.7934e-09
   3.2064e-08
   1.9860e-07
   3.5293e-06
   1.4631e-05
   9.6391e-05
   4.5872e-04

> semilogy(0:14, erro, 'o-')

O gráfico produzido foi

Observando o gráfico, no início vemos um decrescimento linear. Repare que o gráfico foi feito em escala logarítmica. Nesta escala, toda função polinomial é exibida com um aspecto linear, porém com diferentes inclinações, que dependem do grau do polinômio. Será que a inclinação observada no início do gráfico acima está correta? Não custa conferir. O gráfico abaixo exibe em como referência o valor de $M_3h^2/6$ (em azul), que deveria ser o limitante superior para o erro.

De fato, até $h=10^{-5}$ os resultados numéricos parecem consistentes com a teoria. Porém para $h=10^{-6}$ o erro numérico está acima do previsto, uma vez que a curva azul seguira reduzindo indefinidamente. O que pode ter acontecido? Será que fizemos alguma hipótese que foi violada?

Essa discrepância entre teoria e experimento numérico vem do fato de que a função $f$ não é computada exatamente, uma vez que todas as contas são feitas em precisão finita. No caso das funções trigonométricas, mesmo os erros sendo pequenos, próximos à precisão da máquina, eles estão presentes. Ou seja, $$\fl(f(x)) = f(x) + \epsilon(x), \quad |\epsilon(x)| \le \epsilon_f.$$ Desta forma, temos que a diferença centrada na avaliação numérica de $f$ é dada por \begin{align*} D[ f(x) + \epsilon(x) ] & = {f(x + h) + \epsilon(x+h) - f(x-h) - \epsilon(x-h)\over 2h} \\[2mm] & = D_f(x) + D_\epsilon(x). \end{align*} Subtraindo $f'(x)$ em ambos os lados da igualdade acima e tirando o módulo, temos que o erro em precisão finita satisfaz \begin{align*} |D[ f(x) + \epsilon(x) ] -f'(x)| & \le |D_f(x) - f'(x)| + |D_\epsilon(x)| \\[2mm] & \le {M_3h^2\over 6} +{\epsilon_f\over h}, \end{align*} onde usamos que $|D_\epsilon(x)| \le \epsilon_f/h$, uma vez que não sabemos o sinal que $\epsilon(x\pm h)$ assume.

A última desigualdade nos mostra que o limitante para o erro, no caso da função ser avaliada em precisão finita, tem um termo a mais quando comparado com o limitante para o caso de avaliação exata de função, expresso em \eqref{limitanteexato}. Este termo extra, $\epsilon_f/h$, cresce a medida que $h$ diminui, e, em algum momento, passa a dominar o valor do erro. Para as funções trigonométricas $\epsilon_f \approx u \approx 10^{-16}$, a unidade de arredondamento da máquina em precisão dupla.

Agora é fácil se perguntar para que valor de $h$ o limitante para o erro, no caso em que as operações são feitas em um sistema de ponto flutuante, é o menor possível. Denotando esse limitante por \begin{equation} \label{erropf} E(h) = {M_3h^2\over 6} +{\epsilon_f\over h}, \end{equation} seu mínimo é determinado impondo a condição de que $E'(h^*) = 0$. Com um pouco de álgebra, chegamos ao valor ótimo para o passo de discretização, dado por $$h^* = \sqrt[3]{3\epsilon_f \over M_3} \approx 7\cdot 10^{-6},$$ uma vez que $M_3=1$. Isto explica o crescimento observado nos gráficos acima, quando $h$ assumia valores menores que $10^{-6}$ (pontos mais à direita no gráfico).

Para encerrar a questão, no gráfico abaixo novamente exibimos o erro experimental no aproximação da derivada pela fórmula de diferenças centradas (curva em preto), mas agora exibimos também o gráfico de $E(h)$, o limitante para o erro, quando o impacto da computação em precisão finita é levado em consideração, dado por \eqref{erropf}. Fica evidente que o erro experimental está de fato limitado por $E(h)$.

1. Tente reproduzir no Octave o experimento que foi descrito nesta aula. A escolha do ponto $x = \pi/3.2$ foi importante? Os resultados teriam sido diferentes se outro ponto fosse fixado?

2. Na discussão que fizemos, tomamos $M_3 = 1$. Porém, lembrando da aula sobre polinômios de Taylor, para o caso em que $x$ está fixo, e apenas $h$ está variando, poderíamos tomar $$ M_3 = \max_{x-h\le y \le x+h} |f'''(y)|. $$ A medida que $h$ diminui, $M_3$ tende ao valor de $|f'''(x)|$. Refaça o último gráfico do texto, usando esse valor para $M_3$, no cálculo da curva azul. O que você observa?

3. Experimente descobrir empiricamente qual o valor ótimo de $h$ para estimar $f'(20)$ para as funções abaixo, usando a aproximação de diferenças centradas. Qual o menor erro absoluto atingido em cada caso? Qual o menor erro relativo em cada caso?

1. $f(x) = x^3$
2. $f(x) = \exp(\sqrt{x})$