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).