5.7 Integração Numérica — Quadratura tanh-sinh

A quadratura tanh-sinh é um método de integração numérica que utiliza a seguinte mudança de variáveis para transformar uma integral definida no intervalo \([-1,1] \) em um somatório infinito:

\( x = tanh \left(\frac{1}{2}~\pi~sinh~t \right) \)

Após essa transformação, os integrandos decaem com uma proporção exponencial. Para um dado passo de tamanho \(h \) , a integral pode ser aproximada por
\( \int_{-1}^1 f(x) dx \approx \sum_{k=-\infty}^{\infty} H_k f(a_k) \)

com as abcissas
\( a_k = tanh\left(\frac{1}{2} \pi ~ sinh(k h) \right) \)

e pesos
\( H_k = \dfrac{\frac{1}{2} h \pi ~ cosh(k h)}{cosh^2 \left(\frac{1}{2} \pi ~ sinh(k h) \right)} \)

Assim como na quadratura gaussiana, a quadratura tanh-sinh é adequada para integração de funções quaisquer. A convergência dessa fórmula é exponencial em base dois (em relação ao número de divisões). Ou seja, duplicando-se \(n \) temos que o número de dígitos corretos dobra também, sendo o método de quadratura que possui a maior velocidade de convergência conhecido.

 

Escolha do tamanho do passo de discretização (h)

Há uma relação entre o número de pontos usados na soma ( \(n \) ) e o tamanho do passo que deve ser preservada. O \(h \) deve ser escolhido da seguinte forma

\( h = \frac{4}{2^n} \)

 

Implementação

def integral_tanhsinh_quadrature(fun, xi, xe, n=12):
    """
    Numerical integration using Tanh-sinh quadrature
 
    integral = integral_tanhsinh_quadrature(fun, xi, xe, n=12)
 
    INPUT:
      * fun: function to be integrated
      * xi: beginning of integration interval
      * xe: end of integration interval
      * n: number of points used in quadrature
 
    return: \int_{xi}^{xe} f(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/
 
    Jul 2012
    """
    def convert_intervals():
        frac = (xe - xi) / 2.0
        xd = lambda x: (xe + xi + (xe - xi) * x) / 2.0
        f = lambda x: fun(xd(x)) * frac
        return f
 
    def compute_quadrature_weights(n, h):
        w_num = lambda k: 0.5 * h * pi * cosh(k * h)
        w_den = lambda k: cosh(0.5 * pi * sinh(k * h)) ** 2
        w = lambda k: w_num(k) / w_den(k)
        return [w(k) for k in xrange(-n, n)]
 
    def compute_abcissas(n, h):
        x = lambda k: tanh(0.5 * pi * sinh(k * h))
        return [x(k) for k in xrange(-n, n)]
 
    def compute_quadrature(n):
        s = 2 ** n
        h = 4.0 / (2 ** n)
        f = convert_intervals()
        w = compute_quadrature_weights(s, h)
        x = compute_abcissas(s, h)
        quad = [w[k] * f(x[k]) for k in xrange(len(x))]
        return sum(quad)
 
    integral = compute_quadrature(n)
    return integral

Um exemplo de utilização dessas funções pode ser obtido em http://www.sawp.com.br/code/integral/integral_tan_sinh_quadrature.py.

 

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