17 Integração Numérica

Esse capítulo descreve rotinas para executar integração numérica (quadratura) de uma função em uma dimensão. Exitem rotinas para integração adaptativa e não adaptativa de funções genéricas, com rotinas especializadas para casos específicos. essas rotinas especializadas incluem integração sobre intervalos infinitos e semi-infinitos, integrais singulares, incluindo singularidades logaritmicas, cálculos dos valores principais de Cauchy e integrais oscilatórias. A biblioteca reimplementa os algoritmos usados em QUADPACK, um pacote de integração numérica escrito por Piessens, de Doncker-Kapenga, Ueberhuber e Kahaner. Código Fortran para QUADPACK está disponível na Netlib. Também incluídos estão rotinas de integração não adaptativas, de Gauss-Legendre de ordem fixa com coeficientes de alta precisão por Pavel Holoborodko.

As funções descritas nesse capítulo são declaradas no arquivo de cabeçalho ‘gsl_integration.h’.

17.1 Introdução

Cada algoritmo calcula uma aproximação para uma integral definida da forma,

Note que isso é um ou um ou outro restrito, não simultâneo. Para calcular com um erro absoluto especificado, ajuste epsrel para zero. Para calcular com um erro relativo especificado, ajuste epsabs para zero. As rotinas irão falhar na convergência se o erro associado for muito pequeno, mas sempre retorna a melhor aproximação obtida acima daquele estágio.

Os algoritmos em QUADPACK usam uma convenção de nome baseada nas seguintes letras,

Q - rotina de quadratura

N - integrador não adaptativo
A - integrador adaptativo

G - integrando genérico (definido pelo usuário)
W - função peso com integrando

S - singularidades podem ser mais prontamente integradas
P - pontos de dificuldade especial podem ser fornecidos
I - intervalo infinito de integração
O - função de peso oscilatório, cos ou sin
F - integral de Fourier
C - valor principal de Cauchy

Os algoritmos são construídos sobre pares de regras de quadratura, uma regra de alta ordem e uma regra de baixa ordem. Uma regra de alta ordem é usada para calcular a melhor aproximação para uma integral sobre um pequeno intervalo. A diferença entre os resultados de uma regra de alta ordem e a correspondente regra de baixa ordem fornece uma estimativa do erro na aproximação.

17.1.1 Integrandos sem Funções Ponderadas

Os algoritmos para funções gerais (sem uma função peso) são baseados nas regras de Gauss-Kronrod.

Uma regra de Gauss-Kronrod inicia-se com uma regra de quadratura de Gauss clássica de ordem m. Essa regra de quadratura de Gauss clássica é extendida com pontos adicionais entre cada uma das abscissas para fornecer uma regra de alta ordem de Kronrod de ordem 2m+1. A regra de Kronrod é eficiente pelo fato de reutilizar avaliações de funções existentes a partir da regra de Gauss.

A regra de Kronrod de alta ordem é usada como a melhor aproximação para a integral, e a diferença entre as duas regras é usada como uma estimativa do erro na aproximação.

17.1.2 Integrandos com funções ponderadas

Para integrandos com funções peso os algoritmos usam regras de quadratura de Clenshaw-Curtis.

Uma regra de Clenshaw-Curtis inicia com uma aproximação de polinômio de Chebyshev de n-ésima ordem para o integrando. Esse polinômio pode ser integrado exatamente para fornecer uma aproximação para a integral da função original. A expansão de Chebyshev pode ser extendida para ordens mais altas para melhorar a aproximação e fornecer uma estimativa de erro.

17.1.3 Integrandos com funções ponderadas singulares

A presença de singularidades (ou outro comportamento) no integrando pode fazer com que a convergência seja lenta na aproximação de Chebyshev. As regras modificadas de Clenshaw-Curtis usadas em QUADPACK separam muitas funções de peso comuns que causam convergência lenta.

Essas funções de peso são integradas analiticamente contra os polinômios de Chebyshev para pré-calcular momentos de Chebyshev modificados. Combinando os momentos como a aproximação de Chebyshev para a função fornece a integral desejada. O uso de integração analítica para a parte singular da função permite cancelamentos exatos e melhora substancialmente o comportamento geral da convergência da integração.

17.2 Integração QNG Não Adaptativa de Gauss-Kronrod

O algoritmo QNG é um procedimento não adaptativo que usa abscissas fixas de Gauss-Kronrod-Patterson para fazer amostras do integrando em um máximo de 87 pontos. É fornecido para integração rápida de de funções simples.

Function: int gsl_integration_qng (const gsl_function * f, double a, double b, double epsabs, double epsrel, double * result, double * abserr, size_t * neval)

Essa função aplica regras de integração de Gauss-Kronrod em 10-pontos, 21-pontos, 43-pontos e 87-pontos em sucessão até que uma estimativa da integral de f sobre (a,b) seja encontrada dentro dos limites desejados de erro absoluto e erro relativo, epsabs e epsrel. A função retorna a aproximação final, result, uma estimativa de erro absoluto, abserr e o número de avalizações de função usadas, neval. As regras de Gauss-Kronrod são designadas em um tal caminho que cada regra usa todos os resultados de seus predecessores, com o objetivo de minimizar o número total de avalizações de função.

17.3 Integração QAG Adaptativa

O algoritmo QAG é um procedimento simples de quadratura adaptativa. A região de integração é dividida em subintervalos, e a cada iteração o subintervalo com maior estimativa de erro é bissectado. Isso reduz o erro médio rapidamente, uma vez que os subintervalos tornam-se concentrados em volta do local de dificuldades no integrando. Esses subintervalos são gerenciados por uma estrutura gsl_integration_workspace, a qual manipula a memória para a amplitudes de subintervalos, resultados e estimativas de erro.

Function: gsl_integration_workspace * gsl_integration_workspace_alloc (size_t n)

Essa função aloca um espaço de trabalho suficiente para manter n intervalos de precisão dupla, seu resultado de integração e estimativas de erro. Um espaço de trabalho pode ser usado multiplas vezes uma vez que toda a reinicialização necessária é executada automaticamente pelas rotinas de integração.

Function: void gsl_integration_workspace_free (gsl_integration_workspace * w)

Essa função libera a memória assocada ao espaço de trabalho w.

Function: int gsl_integration_qag (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, int key, gsl_integration_workspace * workspace, double * result, double * abserr)

Essa função aplica uma regra de integração adaptativa até que uma estimativa da integral de f no intervalo (a,b) seja encontrada dentro dos limites desejados de erro absoluto e erro relativo, epsabs e epsrel. A função retorna a aproximação final, result, e uma estimativa do erro absoluto, abserr. A regra de integração é determinada pelo valor de key, a qual deve ser escolhida a partir dos seguintes nomes simbólicos,

GSL_INTEG_GAUSS15  (key = 1)
GSL_INTEG_GAUSS21  (key = 2)
GSL_INTEG_GAUSS31  (key = 3)
GSL_INTEG_GAUSS41  (key = 4)
GSL_INTEG_GAUSS51  (key = 5)
GSL_INTEG_GAUSS61  (key = 6)

correspondendo a 15, 21, 31, 41, 51 e 61 pontos para regras de integração de Gauss-Kronrod. As regras de mais alta ordem fornecem melhor precisão para funções bem comportadas, enquanto regras de ordem mais baixa economizam tempo quando a função contém dificuldades localizadas, tais como descontinuidades.

A cada iteração a estratégia de integração adaptativa bisecciona o intervalo com o maior erro estimado. Os subintervalos e seus resultados são armazenados na memória fornecida por workspace. O número máximo def subintervalos é fornecido por limit, o qual pode não exceder o tamanho alocado do espaço de trabalho.

17.4 Integração adaptativa QAGS com singularidades

A presença de uma singularidade integrável na região de integração faz com que uma rotina adaptativa para concentrar novos subintervalos em torno da singularidade. Como o subintervalo decresce em tamanho as sucessivas aproximações da integral convergem para um limite. Essa aproximação de um limite pode ser acelerada usando um procedimento de extrapolação. O algoritmo QAGS combina a bissecção adaptativa ao algoritmo épsilon de Wynn para aumentar a velocidade da integração de muitos tipos de singularidades integráveis.

Function: int gsl_integration_qags (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)

Essa função aplica a regra de integração de 21 pontos de Gauss-Kronrod adaptativamente até que uma estimativa da integral de f no intervalo (a,b) seja encontrada dentro do erro absoluto desejado e dentro dos limites de erro relativo, epsabs e epsrel. Os resultados são extrapolados usando o algoritmo épsilon, o qual acelera a convergência da integral na presença de descontinuidades e de singularidades integráveis. A função retorna a aproximação final a partir da extrapolação, result, e uma estimativa de erro absoluto, abserr. Os subintervalos e seus resultados são armazenados na memória fornecida através de workspace. O número máximo de subintervalos é fornecido através de limit, o qual não pode exceder o tamanho alocado do espaço de trabalho.

17.5 Integração adaptativa QAGP com pontos singulares conhecidos

Function: int gsl_integration_qagp (const gsl_function * f, double * pts, size_t npts, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)

Essa função aplica o algoritmo adaptativo de integração QAGS levando em consideração localizações de pontos singulares fornecidas pelo usuário. O vetor estático pts de comprimento npts deve conter os pontos extremos dos intervalos de integração definidos pela região de integração e pelas localizações da singularidades. Por exemplo, para integrar sobre a região (a,b) com pontos de parada em x_1, x_2, x_3 (onde a < x_1 < x_2 < x_3 < b) o seguinte vetor estático pts deve ser usado

pts[0] = a
pts[1] = x_1
pts[2] = x_2
pts[3] = x_3
pts[4] = b

com npts = 5.

Se você conhece as localizações dos pontos singulares na região de integração então essa rotina irá ser mais rápida que QAGS.

17.6 Integração adaptativa QAGI sobre intervalos infinitos

Function: int gsl_integration_qagi (gsl_function * f, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)

Essa função calcula a integral da função f no intervalo infinito (-\infty,+\infty). A integral é mapeada sobre o intervalo semi-aberto (0,1] usando a transformação x = (1-t)/t,

⌠ ⌡ +∞ −∞  dx  f(x) = ⌠ ⌡ 1 0  dt  (f((1−t)/t) + f(−(1−t)/t))/t2.

É então integrada usando o algoritmo QAGS. A regra normal de 21 pontos de Gauss-Kronrod de QAGS é substituída por uma regra de 15 pontos, pelo fato de a transoformação poder gerar uma singularidade integrável na orígem. Nesse caso uma regra de ordem mais baixa é mais eficiente.

Function: int gsl_integration_qagiu (gsl_function * f, double a, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)

Essa função calcula a integral da função f sobre o intervalo semi-infinito (a,+\infty). A integral é mapeada no intervalo semi-aberto (0,1] usando a transformação x = a + (1-t)/t,

⌠ ⌡ +∞ a  dx  f(x) = ⌠ ⌡ 1 0  dt  f(a + (1−t)/t)/t2

e então integrada usando o algoritmo QAGS.

Function: int gsl_integration_qagil (gsl_function * f, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)

Essa função calcula a integral da função f sobre o intervalo semi-infinito (-\infty,b). A integral é mapeada no intervalo semi-aberto (0,1] usando a transformação x = b - (1-t)/t,

⌠ ⌡ b −∞  dx  f(x) = ⌠ ⌡ 1 0  dt  f(b − (1−t)/t)/t2

e então integrada usando o algoritmo QAGS.

17.7 Integração Adaptativa QAWC em valores principais de Cauchy

Function: int gsl_integration_qawc (gsl_function * f, double a, double b, double c, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)

Essa função calcula o valor principal de Cauchy da integral de f no intervalo (a,b), com uma singularidade em c,

I = ⌠ ⌡ b a  dx  f(x) x − c = lim ϵ→ 0  ⎧ ⎨ ⎩ ⌠ ⌡ c−ϵ a  dx  f(x) x − c + ⌠ ⌡ b c+ϵ  dx  f(x) x − c ⎫ ⎬ ⎭

O algoritmo de bissecção adaptativa QAG é usado, com modificações para garantir que subdivisões não ocorram no ponto singular x = c. Quando um subintervalo contém o ponto x = c ou está muito próximo a ele então uma regra modificada de 25 pontos de Clenshaw-Curtis é usada para controlar a singularidade. Adicionalmente afastando-se da singularidade o algoritmo usa uma regra de integração comum de 15 pontos de Gauss-Kronrod.

17.8 Integrador adaptativo QAWS em funções singulares

O algoritmo QAWS é pensado para integrandos com singularidades algébricas logaritmicas nos extremos de uma região de integração. Com o objetivo de trabalhar eficientemente o algoritmo requer uma tabela pré calculada de momentos de Chebyshev.

Function: gsl_integration_qaws_table * gsl_integration_qaws_table_alloc (double alpha, double beta, int mu, int nu)

Essa função aloca espaço para uma estrutura gsl_integration_qaws_table descrevendo uma função singular de peso W(x) com os parâmetros (\alpha, \beta, \mu, \nu),

W(x) = (x − a)α (b − x)β logμ (x − a) logν (b − x)

onde \alpha > -1, \beta > -1, e \mu = 0, 1, \nu = 0, 1. A função de peso pode receber quatro diferentes formas dependendo dos valores de \mu e \nu,

W(x) = (x − a)α (b − x)β   (μ = 0, ν = 0) W(x) = (x − a)α (b − x)β log(x − a)   (μ = 1, ν = 0) W(x) = (x − a)α (b − x)β log(b − x)   (μ = 0, ν = 1) W(x) = (x − a)α (b − x)β log(x − a) log(b − x)   (μ = 1, ν = 1)

Os pontos singulares (a,b) não são especificados até que a integral seja calculada, onde eles são extremidades do intervalo de integração.

A função retorna um apontador para a tabela mais recentemente alocada gsl_integration_qaws_table se erros não forem encontrados, e 0 em caso de erro.

Function: int gsl_integration_qaws_table_set (gsl_integration_qaws_table * t, double alpha, double beta, int mu, int nu)

Essa função modifica os parâmetros (\alpha, \beta, \mu, \nu) de uma já existente gsl_integration_qaws_table struct t.

Function: void gsl_integration_qaws_table_free (gsl_integration_qaws_table * t)

Essa função libera toda a memória associada a gsl_integration_qaws_table struct t.

Function: int gsl_integration_qaws (gsl_function * f, const double a, const double b, gsl_integration_qaws_table * t, const double epsabs, const double epsrel, const size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)

Essa função calcula a integral da função f(x) sobre o intervalo (a,b) com função de peso singular (x-a)^\alpha (b-x)^\beta \log^\mu (x-a) \log^\nu (b-x). Os parâmetros da função peso (\alpha, \beta, \mu, \nu) são recebidos a partir da tabela t. A integral é,

I = ⌠ ⌡ b a  dx f(x) (x − a)α (b − x)β logμ (x − a) logν (b − x).

O algoritmo de bisecção adaptativa QAG é usado. Quando um subintervalo contém um dos extremos então uma regra modificada especial de 25 pontos de Clenshaw-Curtis é usada para controlar as singularidades. Para subintervalos que não incluem as extremidades uma regra de integração comum de 15 pontos de Gauss-Kronrod é usada.

17.9 Integração adaptativa QAWO em funções oscilatórias

O algoritmo QAWO foi pensado para integrandos com um fator oscilatório, \sin(\omega x) ou \cos(\omega x). Com o objetivo de trabalhar eficientemente o algoritmo requer uma tabela de momentos de Chebyshev a qual deve ser pré-calculada com chamadas às funções abaixo.

Function: gsl_integration_qawo_table * gsl_integration_qawo_table_alloc (double omega, double L, enum gsl_integration_qawo_enum sine, size_t n)

Essa função aloca espaço para uma estrutura gsl_integration_qawo_table e seu espaço de trabalho associado descrevendo uma função de peso seno ou cosseno W(x) com os parâmetros (\omega, L),

W(x) = ⎧ ⎪ ⎨ ⎪ ⎩ sin(ωx) cos(ωx) ⎫ ⎪ ⎬ ⎪ ⎭

O parâmetro L deve ser de comprimento do intervalo sobre o qual a função irá ser integrada L = b - a. A escolha do seno ou do cosseno é feito com o parâmetro sine que deve ser escolhido a partir de um entre os dois seguintes valores simbólicos:

GSL_INTEG_COSINE
GSL_INTEG_SINE

A gsl_integration_qawo_table é uma tabela de coeficientes trigonométricos requerida no processo de integração. O parâmetro n determina o número de níveis de coeficientes que são calculados. Cada nível corresponde a uma bissecção do intervalo L, de forma que n níveis são suficientes para subintervalos abaixo do comprimento L/2^n. A rotina de integração gsl_integration_qawo retorna o erro GSL_ETABLE se o número de níveis for insuficiente para a precisão requisitada.

Function: int gsl_integration_qawo_table_set (gsl_integration_qawo_table * t, double omega, double L, enum gsl_integration_qawo_enum sine)

Essa função modifica os parâmetros omega, L e sine do espaço de trabalho existente t.

Function: int gsl_integration_qawo_table_set_length (gsl_integration_qawo_table * t, double L)

Essa função permite que o parâmetro de comprimento L do espaço de trabalho t seja modificado.

Function: void gsl_integration_qawo_table_free (gsl_integration_qawo_table * t)

Essa função libera toda a memória associada ao espaço de trabalho t.

Function: int gsl_integration_qawo (gsl_function * f, const double a, const double epsabs, const double epsrel, const size_t limit, gsl_integration_workspace * workspace, gsl_integration_qawo_table * wf, double * result, double * abserr)

Essa função usa um algorítimo adaptativo para calcular a integral de f no intervalo (a,b) com a função de peso \sin(\omega x) ou \cos(\omega x) definida pela tabela wf,

I = ⌠ ⌡ b a  dx f(x) ⎧ ⎪ ⎨ ⎪ ⎩ sin(ωx) cos(ωx) ⎫ ⎪ ⎬ ⎪ ⎭

Os resultados foram extrapolados usando o algoritmo épsilon para acelerar a convergência da integral. A função retorna a aproximação final a partir da extrapolação, result, e uma estimativa de erro absoluto, abserr. Os subintervalos e seus resultados são armazenados na memória fornecida por workspace. O número máximo de subintervalos é fornecido por limit, que não deve exceder o tamanho alocado do espaço de trabalho.

Aqueles subtintervalos com “grande” largura d onde d\omega > 4 são calculados usando uma regra de integração de 25 pontos de Clenshaw-Curtis, a qual controla o comportamento oscilatório. Subintervalos com uma “pequena” largura onde d\omega < 4 são calculados usando uma regra de integração de 15 pontos de Gauss-Kronrod.

17.10 Integração adaptativa QAWF de integrais de Fourier

Function: int gsl_integration_qawf (gsl_function * f, const double a, const double epsabs, const size_t limit, gsl_integration_workspace * workspace, gsl_integration_workspace * cycle_workspace, gsl_integration_qawo_table * wf, double * result, double * abserr)

Essa função tenta calcular uma integral de Fourier da função f sobre o intervalo semi-infinito [a,+\infty).

I = ⌠ ⌡ +∞ a  dx f(x) ⎧ ⎪ ⎨ ⎪ ⎩ sin(ωx) cos(ωx) ⎫ ⎪ ⎬ ⎪ ⎭

O parâmetro \omega e a escolha de \sin ou de \cos é feita a partir da tabela wf (o comprimento L pode tomar qualquer valor, uma vez que é sobrescrito por essa função assumindo um valor apropriado para a integração de Fourier). A integral é calculada usando o algoritmo de QAWO sobre cada um dos subintervalos,

C1 = [a, a + c] C2 = [a + c, a + 2c] ... = ... Ck = [a + (k−1) c, a + k c]

onde c = (2 floor(|\omega|) + 1) \pi/|\omega|. A largura c é escolhida para abranger um número ímpar de períodos de forma que as contribuições dos intervalos alternados em sinal e seja monotonamente decrescente quando f for positiva e monotonamente decrescente. O somatório dessa sequência de contribuições é acelerado usando o algoritmo épsilon.

Essa função funciona uma tolerância absoluta geral de abserr. A seguinte estratégia é usada: sobre cada intervalo C_k o algoritmo tenta encontrar a tolerância

TOLk = uk abserr

onde u_k = (1 - p)p^{k-1} e p = 9/10. O somatório das séries geométricas de contribuições a partir de cada intervalo fornece uma tolerância geral de abserr.

Se a integração de um subintervalo levar a dificuldades então a precisão requerida para o subintervalo subsequente é relaxada,

TOLk = uk max (abserr, max i < k  {Ei})

onde E_k é o erro estimado no intervalo C_k.

Os subintervalos e seus resultados são armazenados na memoria fornecida por workspace. O número máximo de subintervalos é fornecido por limit, o qual não pode exceder o tamanho alocado do espaço de trabalho. A integração sobre cada subintervalo usa a memória fornecida por cycle_workspace como espaço de trabalho para o algoritmo QAWO.

17.11 Integração duplamente adaptativa CQUAD

CQUAD é uma nova rotina de quadratura duplamente adaptativa de propósito geral a qual pode tratar a maioria dos tipos de singularidades, valores de função não numéricos tais como Inf ou NaN, bem como algumas integrais divergentes. Geralmente requerem mais avaliações de função que as rotinas de integração em QUADPACK, ainda falha menos e na maioria das vezes por integrandos difíceis.

O algoritmo básico usa um esquema duplamente adaptativo no qual regras de quadratura de Clenshaw-Curtis de incremento do grau são usadas para calcular a interal em cada intervalo. O módulo L_2 da diferença entre os polinômios interpolatórios básicos de duas regras sucessivas é usado como um erro estimado. O intervalo é subdividido se a diferença entre duas regras sucessivas for muito grande ou uma regra de grau máximo tiver sido alcançada.

Function: gsl_integration_cquad_workspace * gsl_integration_cquad_workspace_alloc (size_t n)

Essa função aloca um espaço de trabalho suficiente para manter os dados para n intervalos. O número n não é o número máximo de intervalos que irão ser avaliados. Se o espaço de trabalho estiver cheio, intervalos com memnor erro estimado irão ser descartados. Um mínimo de três intervalos é requerido e para a maioria das funções, um espaço de trabalho de tamanho 100 é suficiente.

Function: void gsl_integration_cquad_workspace_free (gsl_integration_cquad_workspace * w)

Essa função libera a memória assocada ao espaço de trabalho w.

Function: int gsl_integration_cquad (const gsl_function * f, double a, double b, double epsabs, double epsrel, gsl_integration_cquad_workspace * workspace, double * result, double * abserr, size_t * nevals)

Essa função calcula a integral de f no intervalo (a,b) dentro dos limites desejados de erro absoluto e de erro relativo, epsabs e epsrel usando o algoritmo CQUAD. A função retorna a aproximação final, result, uma estimativa do erro absoluto, abserr, e o número de avaliações de função requeridas, nevals.

O algorítimo CQUAD divide a região de integração em subintervalos, e em cada iteração, o subintervalo com o maior erro estimado é processado. O algoritmo usa regras de quadratura de Clenshaw-Curits de graus 4, 8, 16 e 32 sobre 5, 9, 17 e 33 nodos respectivamente. Cada intervalo é inicializado com a regra de menor grau. Quando um intervalo é processado, o maior grau seguinte é avaliado e um erro estimado é calculado baseado no módulo L_2 da diferença entre os polinômios básicos de interpolação de ambas as regras de menor grau e de maior grau seguinte. Se o maior grau disponível tiver sido alcançada, ou os polinômios de interpolação diferirem significativamente, o intervalo é bisectado.

Os subintervalos e seus resultados são armazenados na memória fornecida por workspace. Se o erro estimado ou o número de avaliações de função não forem necessários, os apontadores abserr e nevals podem ser ajustados para NULL.

17.12 Integração de Gauss-Legendre

As rotinas de integração de ordem fixa de Gauss-Legendre são fornecidas para integrações rápidas de funções bem comportadas com ordem polinomial conhecida. Os n pontos da regra de Gauss-Legendre são exatos para polinômios de ordem 2*n-1 ou menor. Por exemplo, essas regras são úteis na integração de funções base para formar matrizes densas para o método de Galerkin. Diferentemente de outras rotinas de integração numérica dentro da biblioteca, essas rotinas não aceitam erros absolutos ou relativos associados.

Function: gsl_integration_glfixed_table * gsl_integration_glfixed_table_alloc (size_t n)

Essa função determina a abscissa de Gauss-Legendre e pesos necessários para um esquema de integração de ordem fixa com n pontos. Se possível, coeficientes de alta precisão pré-calculados são usados. Se pesos pré-calculados não estiverem disponíveis, coeficientes de baixa precisão são calculados imediatamente.

Function: double gsl_integration_glfixed (const gsl_function * f, double a, double b, const gsl_integration_glfixed_table * t)

Essa função aplica a regra de integração de Gauss-Legendre contida na tabela t e retorna o resultado.

Function: int gsl_integration_glfixed_point (double a, double b, size_t i, double * xi, double * wi, const gsl_integration_glfixed_table * t)

Para i em [0, …, t->n - 1], essa função obtém o i-ésimo ponto de Gauss-Legendre xi e peso wi no intervalo [a,b]. Os pontos e pesos são ordenados por valor crescente de ponto. Uma função f pode ser integrada no intervalo [a,b] somando wi * f(xi) sobre i.

Function: void gsl_integration_glfixed_table_free (gsl_integration_glfixed_table * t)

Essa função libera a memória associada à tabela t.

17.13 Códigos de erro

Adicionalmente a códigos de erros padronizados para argumentos inválidos as funções retornam os seguintes valores,

GSL_EMAXITER

o número máximo de subdivisões foi extrapolado.

GSL_EROUND

não posso alcançar a tolerância devido a erro de arredondamento, ou erro de arredondamento foi detectado na tabela de extrapolação.

GSL_ESING

uma singularidade não integrável ou outro mau comportamento do integrado foi encontrado no intervalo de integração.

GSL_EDIVERGE

a intergral é divergente, ou de convergência extremamente lenta para ser integrada numericamente.

17.14 Exemplos

O integrador QAGS irá tratar uma larga classe de integrais definidas. Por exemplo, considere a seguinte integral, a qual tem uma singularidade logarítmica algébrica na orígem,

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f (double x, void * params) {
  double alpha = *(double *) params;
  double f = log(alpha*x) / sqrt(x);
  return f;
}

int
main (void)
{
  gsl_integration_workspace * w 
    = gsl_integration_workspace_alloc (1000);
  
  double result, error;
  double expected = -4.0;
  double alpha = 1.0;

  gsl_function F;
  F.function = &f;
  F.params = &alpha;

  gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
                        w, &result, &error); 

  printf ("result          = % .18f\n", result);
  printf ("exact result    = % .18f\n", expected);
  printf ("estimated error = % .18f\n", error);
  printf ("actual error    = % .18f\n", result - expected);
  printf ("intervals =  %d\n", w->size);

  gsl_integration_workspace_free (w);

  return 0;
}

Os resultados abaixo mostram que a precisão desejada é encontrada após 8 subdivisões.

$ ./a.out 
result          = -3.999999999999973799
exact result    = -4.000000000000000000
estimated error =  0.000000000000246025
actual error    =  0.000000000000026201
intervals =  8

De fato, o procedimento de extrapolação usado por QAGS produz uma precisão de quase o dobro de dígitos. O erro estimado retornado pelo procedimento de extrapolação é maior que o erro atual, fornecendo uma margem de segurança de uma ordem de magnitude.

17.15 Referências e Leituras Adicionais

Os seguinte livro é a referência definitiva para QUADPACK, e foi escrito pelos autores originais. Fornece descrições dos algoritmos, listagem de programas, programas de teste e exemplos. Também inclui dicas úteis sobre integração numérica e muitas referências sobre a literatura de integração numérica usada em desenvolvimento QUADPACK.

• R. Piessens, E. de Doncker-Kapenga, C.W. Ueberhuber, D.K. Kahaner. QUADPACK A subroutine package for automatic integration Springer Verlag, 1983.

O algoritmo de integração CQUAD é descrito no seguinte artigo:

• P. Gonnet, “Increasing the Reliability of Adaptive Quadrature Using Explicit Interpolants”, ACM Transactions on Mathematical Software, Volume 37 (2010), Impressão 3, Artigo 26.

Esse documento foi gerado em 23 de Julho de 2013 usando texi2html 5.0.