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

e
\( \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:

  1. Obter o polinômio de Legendre de ordem \(n \) , onde \(n \) é o número de pontos utilizados.
  2. Encontrar as \(n \) raízes do polinômio obtido.
  3. Encontrar a primeira derivada do n-ésimo polinômio, pois o peso é dependente desta função. Com isso podemos determinar os pesos \(H_j \) .
  4. 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).