GSL - quadraturas

Computação Científica - Um Livro Colaborativo

13.1 GSL - quadraturas

A GSL (GNU Scientific Library) é uma biblioteca livre para C e C++ com milhares de rotinas matemáticas envolvendo funções especiais, mínimos quadrados, integração numéricas, autovalores e autovetores, etc. A página do projeto https://www.gnu.org/software/gsl/ descreve a documentação sobre instalação e uso do pacote e o manual https://www.gnu.org/software/gsl/doc/html/index.html detalha o uso das rotinas.

Vamos começar trabalhando com as rotinas de integração numérica: https://www.gnu.org/software/gsl/doc/html/integration.html. O uso das rotinas de integração exigem a inclusão do pacote gsl_integration.h.

Uma das rotinas mais simplesde integração é a gsl_integration_qng. O nome da rotina tem um significado:

q $\to$ quadrature routine
n $\to$ non-adaptive integrator
g $\to$ general integrand

A sintaxe dessa rotina é:

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

onde

*f é um ponteiro para função do tipo gsl (gsl_function). Essa função aceita passagem de parâmetros.
a é o limite inferior de integração;
b é o limite superior de integração;
epsabs é a tolerância absoluta;
epsrel é o tolerância relativa;
result é o resultado da integral
abserr estimativa do erro absoluto;
neval é o número de vezes que a função foi calculada.

O esquema de integração é Gauss-Konrad com número fixo de pontos: 10, 21, 43 ou 87. Se a integral apresentar estimativa de erro grande, escolha outra rotina de integração.

Exemplo 13.1.1. Calcule a integral $\int_{0}^{1} e^{- x^{2}} d x$

usando a rotina qng da biblioteca gsl.

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

double f(double x, void *params)
{
return exp(-x*x);
}

void main()
{
double params,result,abserr;
size_t neval;
gsl_function F;
F.function=&f;
F.params=&params;
gsl_integration_qng(&F, 0, 1, 1e-10, 1e-10,&result, &abserr, &neval);
printf("resultado=%f, erro=%e, n_calculos=%zu\n",result,abserr,neval);
}

Compile o código com a linha

gcc teste.c -lm -lgsl -lgslcblas

Uma função para ser integrada no gsl deve ter a formato double f(double x, void *params), onde x é a variável de integração e *params é um ponteiro para uma lista de parâmetros da função. A definição da funçao do tipo gsl se dá pela linha gsl_function F;. Em seguida, usa-se a linha F.function=&f; para configurar a função gsl F, associando ao ponteiro para a função f. A linha F.params=&params; atribui o vetor params a lista de parâmetro da função F. No caso acima, nós definimos double params;, mas não atribuímos valor algum, pois a função não tem parâmetros.

Exemplo 13.1.2. Calcule a integral $\int_{0}^{1} e^{- λ x^{2}} d x$

usando $λ = . 5$ e $λ = 2$ .

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

double f(double x, void *params)
{
double lambda=*((double *) params);
return exp(-lambda*x*x);
}

void main()
{
double params,result,abserr;
size_t neval;
gsl_function F;
F.function=&f;
F.params=&params;

params=.5;
gsl_integration_qng(&F, 0, 1, 1e-10, 1e-10,&result, &abserr, &neval);
printf("resultado=%f, erro=%f, n_calculos=%zu\n",result,abserr,neval);

params=2;
gsl_integration_qng(&F, 0, 1, 1e-10, 1e-10,&result, &abserr, &neval);
printf("resultado=%f, erro=%f, n_calculos=%zu\n",result,abserr,neval);
}

A rotina autoadaptativa mais simples é a gsl_integration_qag:

q $\to$ quadrature routine
a $\to$ adaptive integrator
g $\to$ general integrand

Ela consiste em dividir o intervalo de integração ao meio e, depois, cada um dos subintervalos ao meio outra vez, até a convergência. A sintaxe é

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)

onde

*f é um ponteiro para função do tipo gsl (gsl_function). Essa função aceita passagem de parâmetros.
a é o limite inferior de integração;
b é o limite superior de integração;
epsabs é a tolerância absoluta;
epsrel é o tolerância relativa;
limit é o máximo de subintervalos;
key é a regra de integração
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);
* workspace é o espaço na memória para dividir o intervalo e recalcular a integral;
result é o resultado da integral
abserr estimativa do erro absoluto;

Exemplo 13.1.3. Calcule a integral $\int_{0}^{1} \frac{1}{μ} e^{- \frac{r}{μ}} d μ$

para $r = 0.5$ e $r = 1$ .

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

double f(double mu, void *params)
{
double r=*((double *) params);
return exp(-r/mu)/mu;
}

void main()
{
double params,result,abserr;

gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
gsl_function F;
F.function=&f;
F.params=&params;

params=.5;
gsl_integration_qag(&F, 0, 1, 1e-10, 1e-10, 1000,6,w ,&result,&abserr);
printf("resultado=%f, erro=%e \n",result,abserr);

params=1;
gsl_integration_qag(&F, 0, 1, 1e-10, 1e-10, 1000,6,w ,&result,&abserr);
printf("resultado=%f, erro=%e \n",result,abserr);

gsl_integration_workspace_free(w);
}

A função gsl_integration_qags é adequada para integrandos com singularidades e gsl_integration_qagp é adequada para integrandos com singularidades conhecidas. Na sintaxe da primeira função não passa a regra de integração e, na segunda, passa os pontos onde a função é singular junto com os limites de integração. Observe a resolução do exemplo 13.1.3 com essas funções.

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

double f(double mu, void *params)
{
double r=*((double *) params);
return exp(-r/mu)/mu;
}

void main()
{
double params,result,abserr,pontos[2]={0,1};

gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
gsl_function F;
F.function=&f;
F.params=&params;

params=.5;
gsl_integration_qags(&F, 0, 1, 0, 1e-10, 1000,w ,&result,&abserr);
printf("resultado=%f, erro=%e \n",result,abserr);

params=1;
gsl_integration_qagp(&F,pontos,2, 0, 1e-10, 1000,w ,&result,&abserr);
printf("resultado=%f, erro=%e \n",result,abserr);

gsl_integration_workspace_free(w);
}

Agora, estudar as rotinas do gsl que aplicam a quadratura de Gauss-Legendre. As quatro funções têm as seguintes sintaxes:

gsl_integration_glfixed_table * gsl_integration_glfixed_table_alloc(size_t n)

void gsl_integration_glfixed_table_free(gsl_integration_glfixed_table * t)

A primeira função calcula as abscissas e pesos da quadratura de Gauss-Legendre com $n$ pontos. A segunda libera o espaço alocado na memória. A terceira função tem a sintaxe

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

onde o retorno da função é o valor da integral no intervalo $[a, b]$ usando os pesos e abscissas salvos na tabela t. A última função tem sintaxe

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

onde

i, i=0,1,2,...n-1 é o índice que indica o i-ésimo ponto;
xi é a i-ésima abscissa;
wi é o i-ésima peso;
gsl_integration_glfixed_table * t é a tabela com os pesos e abscissas.

Exemplo 13.1.4. Calcule a integral $\int_{0}^{1} e^{- λ x^{2}} d x$

usando a quadratura de Gauss-Legendre.

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

double f(double mu, void *params)
{
double r=*((double *) params);
return exp(-r/mu)/mu;
}

void main()
{
  int N=200;
  gsl_integration_glfixed_table * table;
  table = gsl_integration_glfixed_table_alloc (N);
  double result,params;

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

  params=.5;
  result=gsl_integration_glfixed(&F, 0, 1, table);
  printf("O valor da integral é %f\n",result);

  gsl_integration_glfixed_table_free(table);
}

Alternativamente, podemos salvar os pesos e abscissas em vetores integrar usando o somatório $\sum_{i = 0}^{n - 1} w_{i} f (x_{i}) .$

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

double f(double mu, void *params)
{
double r=*((double *) params);
return exp(-r/mu)/mu;
}

void main()
{
  int N=200,flag;
  gsl_integration_glfixed_table * table;
  table = gsl_integration_glfixed_table_alloc (N);
  double result=0,a[N],p[N],params=.5;
  size_t i;
  for (i=0; i<N;i++)
  {
    flag = gsl_integration_glfixed_point (0, 1, i, &a[i], &p[i], table);
    printf("%zu %f %f\n", i, a[i], p[i]);
  }
  for (i=0; i<N;i++) result+=p[i]*f(a[i],&params);
  printf("O valor da integral é %f\n",result);

  gsl_integration_glfixed_table_free(table);
}

Este texto é disponibilizado nos termos da licença Creative Commons Atribuição-CompartilhaIgual 3.0 Não Adaptada (CC-BY-SA 3.0). Página gerada em 19/8/2020 às 17:33:17.

Recursos

Projeto

IME - UFRGS
- Página do IME
- Página da UFRGS

UFRGS - IME - Recursos Educacionais Abertos de Matemática. Contato: reamat@ufrgs.br.

Informe erros ou edite você mesmo!