5.6 Integração Numérica — Quadratura de Gauss

1. Introdução

Consideremos o problema de encontrar uma aproximação numérica $I(f)$ para a integral da função linear $f(x)$

$$ I(f) = \int_a^b f(x) dx $$

onde $-\infty \leq a \leq b \infty$ . Como $I(f)$ é linear. É coerente aproximarmos isso por um funcional também linear,

$$ Q(f) = \sum_{i=0}^m \sum_{j=1}^n A_{ij} f^{(i)}(a_{ij}) $$

onde $Q(f)$ é chamada de problema da quadratura. Isto é,

$$ I(f) = Q(f) + E $$

Note que quando $E=0$ , temos a aproximação da integral como uma combinação linear dos valores de $f(x)$ e suas derivadas. O problema da quadratura numérica consiste em especificar $A\_{ij}$ e $a\_{ij}$ tais que esta aproximação obedeçam a certas propriedades, tais como convergência e acurácia.

Nossa abordagem consiste em definir uma aproximação polinomial de $f(x)$ de forma a escolher $A\_{ij}$ e $a\_{ij}$ tais que minimizem $E$ para um aproximante de grau suficientemente baixo. Assim como fizemos com a regra do trapézio, representaremos a integral como uma combinação linear de valores funcionais

únicos. Ou seja, reescreveremos a integral como

$$ \int_a^b f(x) dx = \sum_{j=1}^{n} H_j f(a_j) + E $$

Vamos assumir que os limites de integração $a$ e $b$ da equação acima são finitos. Se a equação acima é exata para polinômios de grau $2n-1$ ou menores, podemos utilizar um conjunto de $2n$ equações para descobrir $2n$ constantes ao substituir $f(x) = x^k$ para $k=0,1,\ldots,2n-1$ na equação acima. Setando $E=0$ , temos que

$$ \alpha_k = \sum_{j=1}^{n} H_j a_j^k $$

para $k=0,1,\ldots,2n-1$ , onde

$$ \alpha_k = \int_a^b x^k dx = \frac{b^{k+1} – a^{k+1}{k+1} $$

Seja a fórmula de interpolação de Hermite:

$$ f(x) = \sum_{j=1}^n h_j f(a_j) + \sum_{j=1}^n \hat{h_j(x)} f'(a_j) + \frac{pn^2(x)}{(2n)!} f^{(2n)}(\xi) $$

que é exata para os polinômios de grau $2n-1$ ou menores. Integrando essa equação, temos

$$ \int_a^b f(x) dx = \sum_{j=1}^n H_j f(a_j) + \sum_{j=1}^n \hat{H_j} f'(a_j) + E $$

onde

$$ H_j = \sum_a^b h_j(x) dx $$

$$ \hat{H_j} = \int_a^b \hat{h_j}(x) dx $$

Como $E$ é zero, se $f(x)$ é um polinômio de grau $2n-1$ ou menor e se podemos escolher as abcissas tal que $\hat{H_j} = 0$ , então teremos as propriedades desejadas.

2. Quadratura de Gauss-Legendre

A quadratura de Gauss-Legendre consiste em aproximar a integral $I(f)$ utilizando $n$ pontos

$$ I(f) = H_0 f(a_0) + H_1 f(a_1) + \cdots + H_{n-1} f(a_{n-1}) $$

onde as abcissas $a_j$ são as raízes do polinômio de Legendre e os pesos são obtidos por

$$ H_j = \dfrac{-2}{(n+1) P_{n+1}(a_j)P’_n(a_j)} = \dfrac{2}{(1-a_j^2) P’_n(a_j)} $$

3. Implementação

Como podemos notar das seções anteriores, a quadratura requer avaliar a função em pontos bem definidos utilizando pesos específicos. Como dito na última seção, os pontos a serem utilizados $a\_j$ são as raízes do polinômio de Legendre. Além disso, podemos ver que os pesos também dependem dessas raízes $a\_j$ . Sendo assim, o problema da quadratura de Gauss-Legendre é resolvido em resolver quatro grandes etapas:

Obter o polinômio de Legendre de ordem $n$ , onde $n$ é o número de pontos utilizados.
Encontrar as $n$ raízes do polinômio obtido.
Encontrar a primeira derivada do n-ésimo polinômio, pois o peso é dependente desta função. Com isso podemos determinar os pesos $H_j$ .
Computar a quadratura gaussiana, utilizando os pesos $H\_j$ e as raízes $a\_j$ .

3.1. Computando o Polinômio de Legendre

O primeiro passo é encontrar os polinômios de Legendre para que possamos encontrar suas raízes. Como sabemos que nos passos posteriores teremos que encontrar a derivada do n-ésimo polinômio, escolhemos gerar seus coeficientes e colocá-los em uma lista ordenada da maior ordem para a menor. Por exemplo, o polinômio

$$ p(x) = 3x^4 + 2x^3 – x + 7 $$

pode ser representado em forma de lista

$$ p(x) = [3, 2, 0, -1, 7] $$

onde a ordem fica implícita na posição do vetor em ordem decrescente. Assim, a função que computa os coeficientes, representando o polinômio como uma lista é

def compute_legendre_polynomials_coeficients(n):
    """
    Compute the coefficients of the nth Legendre polynomial.

    [coefficients] = compute_legendre_polynomials_coeficients(n)

    INPUT:
      * n: order of the Legendre polynomial

    return: a list containing the polynomial coefficients of nth Legendre
            polynomial (ordered from the highest to lowest)

    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
    """
    b = 2.0 ** n
    B = binomial_coefficients
    f = lambda k: b * B(n, k) * B((n + k - 1) / 2.0, n)
    a = [f(k) for k in xrange(0, n + 1)]
    a.reverse()
    return a

Note que esta função requer o cálculo do binômio $\binom{n}{k}$ . A função que faz esta computação é

def binomial_coefficients(n, k):
    """
    Calcule any binomial coefficients.

    coeff = binomial coefficients(n, k)

    INPUT:
      * n: power of binomial
      * k: coefficient of x^k term

    return: binomial coefficient

    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
    """
    fun = lambda i: (n - i + 1.0) / i
    l = [1.0] + [fun(i) for i in xrange(1, k + 1)]
    coeff = reduce(mul, l)
    return coeff

3.2. Computando as Raízes do Polinômio de Legendre

O próximo passo é obter as raízes $a_j$ do polinômio de Legendre. Para isso, precisamos de um algoritmo que calcule as $n$ raízes de um polinômio de ordem $n$ . Nesse caso, podemos utilizar qualquer método numérico que realize esta tarefa: Aberth, Bairstow, Laguerre, Jenkins-Traub, etc.

Como apresentado em outros posts, cada um desses métodos possui vantagens e desvantagens, tais como requerer ou não derivadas, complexidade, etc. Sabendo que todas as raízes do polinômio de Legendre estão definidas no intervalo $[-1,1]$ e sabendo os coeficientes do polinômio, o que nos permite derivá-lo, escolhemos o método de Aberth. Esse método possui uma convergência quadrática e requer a função (polinômio) e a sua derivada. Assim, podemos computar todas as raízes do polinômio utilizando a seguinte função:

def compute_legendre_roots(n):
    """
    Compute all roots of the nth Legendre polynomial

    [roots] = compute_legendre_roots(n)

    INPUT:
      * n: order of Legendre polynomial

    return: the roots of P_n(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
    """
    legendre_coefficients = compute_legendre_polynomials_coeficients(n)
    diff_coefs = polynomial_derivative_coefficients(legendre_coefficients)
    legendre_fun = polynomial_coeffs_to_function(legendre_coefficients)
    diff_legendre_fun = polynomial_coeffs_to_function(diff_coefs)
    roots = aberth(legendre_fun, diff_legendre_fun, polDegree=n)
    return roots

Onde polynomial_derivative_coefficients realiza a derivada de um polinômio representado como lista e retorna o outro polinômio derivado também em forma de lista:

def polynomial_derivative_coefficients(coeff):
    """
    Compute the derivative of a polynomial represented by a list.

    [diff_coefficients] = polynomial_derivative_coefficients(coeff)

    INPUT:
      * coeff: coefficients of the polynomial

    return: a list containing the polynomial coefficients of the derivated
            polynomial (ordered from the highest to lowest)

    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
    """
    order = len(coeff) - 1
    f = lambda i: (order - i) * coeff[i]
    diff = [f(i) for i in xrange(order)]
    return diff

A função polynomial_coeffs_to_function converte uma representação em forma de lista em uma função avaliável em qualquer ponto:

def polynomial_coeffs_to_function(coeff):
    """
    Convert the polynomial coefficients in a polynomial function.

    function = polynomial_coeffs_to_function(coeff)

    INPUT:
      * coeff: coefficients of the polynomial

    return: a function

    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
    """
    order = len(coeff) - 1
    f = lambda x: sum([coeff[i] * x ** (order - i) for i in xrange(order + 1)])
    return f

3.3. Computando os Pesos da Quadratura

Conforme a função apresentada, os $n$ pesos obtidos a partir das $n$ raízes são calculados pela seguinte função:

def compute_quadrature_weights(n):
    """
    Compute the weights of each coefficient used in gaussian quadrature.

    [weights] = compute_quadrature_weights(n)

    INPUT:
      * n: order of Legendre polynomial

    return: a list with the weights

    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 get_diff_Legendre_n(n):
        l1 = compute_legendre_polynomials_coeficients(n)
        diff_coefs = polynomial_derivative_coefficients(l1)
        return polynomial_coeffs_to_function(diff_coefs)

    def compute_weights(n):
        p1n = get_diff_Legendre_n(n)
        a = compute_legendre_roots(n)
        h = lambda i: 2.0 / ((1.0 - a[i] * a[i]) * p1n(a[i]) * p1n(a[i]))
        H = [h(j) for j in xrange(n)]
        return H

    return compute_weights(n)

3.4. Computando a Integral

O último passo é computar a integral da função escolhida utilizando a quadratura de Gauss-Legendre:

def integral_gaussian_quadrature(fun, xi, xe, order=10):
    """
    Numerical integration using Gaussian quadrature

    integral = integral_gaussian_quadrature(fun, xi, xe, order=10)

    INPUT:
      * fun: function to be integrated
      * xi: beginning of integration interval
      * xe: end of integration interval
      * order: degree of Legendre polynomial

    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():
        f = convert_intervals()
        roots = compute_legendre_roots(order)
        weights = compute_quadrature_weights(order)
        approx = lambda i: weights[i] * f(roots[i].real)
        quad = [approx(i) for i in xrange(order)]
        return sum(quad)

    integral = compute_quadrature()
    return integral

Onde convert_intervals realiza uma mudança de variáveis para converter um intervalo qualquer para o intervalo compatível com o da integral demonstrada no método, isto é, $[-1,1]$ .

Um exemplo de utilização dessas funções pode ser obtido em http://www.sawp.com.br/code/integral/integral_gauss_legendre_quadrature.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).