4.2 Derivação Numérica — Extrapolação de Richardson

1. Introdução

No post anterior, mostramos que as estimativas das derivadas podem ser melhoradas ao diminuir o passo $h$ ou ao usar mais pontos através de sucessivas substituições na fórmula de Taylor. Além disso, na seção de implementação, comentamos que a escolha de $h$ pode acarretar em problemas, uma vez que o computador possui precisão limita. Uma abordagem para se melhorar a estimativa da derivada consiste na utilização da extrapolação de Richardson.

No caso geral da extrapolação de Richardson, temos um procedimento para aproximar um funcional $\phi = \phi(f)$ de uma dada função $f$ por uma sequência, dependente do tamanho de $h$ , em que, como $h \rightarrow 0$ , o erro possui a seguinte forma assintótica

$$ E = \sum_{j=1}^{\infty} a_j h^{\gamma_j} $$

onde $\gamma\_1 \lt \gamma\_2 \lt \cdots$ , os $a\_j$ são constantes que podem depender da função $f$ mas não dependem de $h$ . Tais constantes $a\_j$ nos são desconhecidas, mas assumimos conhecer $\gamma\_j$ . Se a computação é gerada para dois valores distintos de $h$ , $h\_1$ e $h\_2$ , tal que $h\_1 \gt h\_2$ , então há duas aproximações distintas $\phi\_1$ e $\phi_2$ , que se relacionam com o valor verdadeiro de $\phi$ pela fórmula

$$ \phi = \phi_i + \sum_{j=1}^{\infty} a_j h_i^{\gamma_j} $$

para $i=1,2$ . Multiplicando essa equação para quando $i=1$ por $h\_2^{\gamma\_1}$ e por $i=2$ por $h\_1^{\gamma\_1}$ , e subtraindo, temos que $\phi$ é dado por

$$ \phi = \frac{1}{h_1^{\gamma_1} – h_2^{\gamma_1} (h_1^{\gamma_1} \phi_2 + h_2^{\gamma_1} \phi_1) + \sum_{j=2}^{\infty} a_j \frac{h_1^{\gamma_1} h_2^{\gamma_j} – h_2^{\gamma_1} h_1^{\gamma_j} {h_1^{\gamma_1} – h_2^{\gamma_1} $$

Escolhendo-se $h\_2 = \rho h\_1$ , $\rho > 1$ , na equação acima, temos

$$ \phi = \frac{\phi_2 – \rho^{\gamma_1}\phi_1}{1 – \rho^{\gamma_1} + \sum_{j=2}^{\infty} \frac{\rho^{\gamma_j} – \rho^{\gamma_1}{1 – \rho^{\gamma_1} a_j h_1^{\gamma_j} = \phi_{12} + \sum_{j=2}^{\infty} b_j h_1^{\gamma_j} $$

Se desejarmos eliminar p termo com $j=2$ nessa equação, devemos computar $\phi\_3$ com $h\_3 = \rho h\_2$ e computar $\phi\_{23}f$ de forma simular para, então, calcular $\phi\_{123}$ a partir de $\phi\_{12}$ e $\phi_{23}$.

Formalmente, denotamos $\phi\_i$ por $T\_0^i$ e geramos um vetor triangular de aproximantes $T_m^i$ pela fórmula

$$ T_m^i = \frac{T_{m-1}^{i+1} – \rho^{\gamma_m}T_{m-1}^i}{1 – \rho^{\gamma_m} = T_{m-1}^{i+1} + \frac{T_{m-1}^{i+1} – T_{m-1}^i}{\rho^{-\gamma_m} – 1} $$

onde o elemento $T\_{m-1}^1$ corresponde a $\phi\_{12 \cdots m}$.

Em alguns casos, um crescimento geométrico de $\frac{1}{h}$ pelo fator $\frac{1}{\rho}$ não é desejável ou possível, e nós estamos interessados em uma sequência geral ${h\_i}$ . Neste caso, podemos gerar uma tabela $T\_m^i$ por um simples algoritmo apenas se $\gamma\_j$ tem uma estrutura especial, $y\_j = j \gamma + \delta$ . Para o caso usual $\delta=0$ , temos que

$$ T_m^i = \frac{h_j^{\gamma} T_{m-1}^{i+1} – h_{i+m}^{\gamma} T_{m-1}^i}{h_i^{\gamma} – h_{i+m}^{\gamma}= T_{m-1}^{i+1} + \frac{T_{m-1}^{i+1} – T_{m-1}^{i} {(h_i/h_{i+m}^{\gamma}) – 1} $$

 

2. Extrapolação de Richardson

A derivada numérica $D$ de uma função pode ser descrita como função do passo $h$ das fórmulas apresentadas no post anterior na forma

$$ D = D(h) + E(h) $$

em que $D$ é o valor exato da derivada (analítica), $E(h)$ é o resíduo que aparece do truncamento da série de Taylor e $D(h)$ é a derivada aproximada. Se fizermos duas estimativas independentes usando passos distintos, $h\_1$ e $h\_2$, teremos que

$$ D = D(h_1) + E(h_1) $$

e

$$ D = D(h_2) + E(h_2) $$

portanto

$$ D(h_1) + E(h_1) = D(h_2) + E(h_2) $$

Supondo que as fórmulas usadas na diferenciação numérica sejam aquelas com erros $O(h^2)$ , temos que

$$ E(h_1) \approx E(h_2) \approx h^2 $$

logo, a razão entre os dois erros será

$$ \frac{E(h_2)}{E(h_1)} \approx \frac{h_2^2}{h_1^2} $$

Com isso, temos a seguinte relação entre os erros

$$ E(h_1) \approx E(h_2) \left( \frac{h_1}{h_2} \right)^2 $$

que podemos substituir na equação da igualdade das derivadas numéricas

$$ D(h_1) + E(h_2) \left( \frac{h_1}{h_2} \right)^2 \approx D(h_2) + E(h_2) $$

isolando o erro, temos

$$ E(h_2) \approx \frac{D(h_1) – D(h_2)}{1 – (h1/h2)^2} $$

Com isso, podemos obter a estimativa para o erro de truncamento em termos da relação dos tamanhos dos passos. Essa estimativa pode ser substituída em

$$ D = D(h_2) + E(h_2) $$

o que permite fornecer uma estimativa melhorada da derivada

$$ D \approx D(h_2) + \frac{1}{(h1/h2)^2 – 1} \left( D(h_2) – D(h_1) \right) $$

 

3. Implementação

def richard_extrapolation(diff, f, x, h1=0.1, h2=0.05):
    """
    Increase the accuracy of numerical derivative (diff).

    diff_f_xi = richard_extrapolation(diff, f, x, h1=0.001, h2=0.0005)

    INPUT:
      * diff: the numerical differentiation function (diff1, diff2, ...)
      * f: function to derivate
      * x: evaluation point
      * h: step size

    return: f^(i)(X=x)

    Author: Pedro Garcia [sawp@sawp.com.br]
    see: http://www.sawp.com.br

    License: Creative Commons
             http://creativecommons.org/licenses/by-nc-nd/2.5/br/

    Jun 2012
    """
    Dh2 = diff(f, x, h2)
    Dh1 = diff(f, x, h1)
    extrapolation = Dh2 + (1.0 / ((h1 / h2) ** 2) - 1.0) * (Dh2 - Dh1)
    return extrapolation

Um exemplo de utilização dessa função pode ser obtido em http://www.sawp.com.br/code/derivatives/derivative_richard_extrapolation.py .

 

4. Copyright

Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.

References

[1] Anthony Ralston and Philip Rabinowitz, A First Course in Numerical Analysis (2nd ed.), McGraw-Hill and Dover, (2001).
[2] N.B Franco, Cálculo Numérico, Pearson Prentice Hall (2006).
</fieldset>