Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
Derived types and interfaces | Functions/Subroutines
sll_m_gauss_lobatto_integration Module Reference

Description

Gauss-Lobatto integration.

Low-level mathematical utility that applies the Gauss-Lobatto method to compute numeric integrals.

Derived types and interfaces

interface  sll_o_gauss_lobatto_integrate_1d
 Integrate numerically with Gauss-Lobatto formula. More...
 

Functions/Subroutines

real(kind=f64) function gauss_lobatto_integral_1d (f, a, b, n)
 Gauss-Lobatto Quadrature. More...
 
real(kind=f64) function, dimension(2, n), public sll_f_gauss_lobatto_points_and_weights (n, a, b)
 Returns a 2d array of size (2,n) containing gauss-lobatto points and weights in the interval [a,b]. More...
 
real(kind=f64) function, dimension(n), public sll_f_gauss_lobatto_points (n, a, b)
 Returns a 1d array of size (n) containing gauss-lobatto points in the interval [a,b]. More...
 
real(kind=f64) function, dimension(n), public sll_f_gauss_lobatto_weights (n, a, b)
 Returns a 1d array of size (n) containing gauss-lobatto weights in the interval [a,b]. More...
 
function, public sll_f_gauss_lobatto_derivative_matrix (n, y)
 Construction of the derivative matrix for Gauss-Lobatto, The matrix must be already allocated of size \( n^2 \). More...
 
subroutine dlob (n, dalpha, dbeta, dleft, dright, dzero, dweigh, ierr, de, da, db)
 This comes from http://dl.acm.org, Algorithme 726 : ORTHPOL, appendices and supplements To use those functions, READ the documentation beside and find more information about coefficients in paper Algorithm 726 - ORTHPOL: A package of routines for generating orthogonal polynomials and Gauss-type quadrature rules by Walter Gautschi (here xxx is 726 in other references) formulas (1.1) to (1.3) page 2, and book Numerical Mathematics by Alfio Quarteroni, Riccardo Sacco and Fausto Saleri section 10. More...
 
subroutine dgauss (n, dalpha, dbeta, deps, dzero, dweigh, ierr, de)
 Given n and a measure dlambda, this routine generates the n-point Gaussian quadrature formula. More...
 

Function/Subroutine Documentation

◆ dgauss()

subroutine sll_m_gauss_lobatto_integration::dgauss ( integer(kind=i32)  n,
real(kind=f64), dimension(n)  dalpha,
real(kind=f64), dimension(n)  dbeta,
real(kind=f64)  deps,
real(kind=f64), dimension(n)  dzero,
real(kind=f64), dimension(n)  dweigh,
integer(kind=i32)  ierr,
real(kind=f64), dimension(n)  de 
)
private

Given n and a measure dlambda, this routine generates the n-point Gaussian quadrature formula.

integral over supp(dlambda) of f(x)dlambda(x)

   = sum from k=1 to k=n of w(k)f(x(k)) + R(n;f).

The nodes are returned as zero(k)=x(k) and the weights as weight(k)=w(k), k=1,2,...,n. The user has to supply the recursion coefficients alpha(k), beta(k), k=0,1,2,...,n-1, for the measure dlambda. The routine computes the nodes as eigenvalues, and the weights in term of the first component of the respective normalized eigenvectors of the n-th order Jacobi matrix associated with dlambda. It uses a translation and adaptation of the algol procedure imtql2, Numer. Math. 12, 1968, 377-383, by Martin and Wilkinson, as modified by Dubrulle, Numer. Math. 15, 1970, 450. See also Handbook for Autom. Comput., vol. 2 - Linear Algebra, pp.241-248, and the eispack routine imtql2.

   Input:
         - n - - the number of points in the Gaussian quadrature
                 formula; type integer
         - alpha,beta - - arrays of dimension  n  to be filled
                 with the values of  alpha(k-1), beta(k-1), k=1,2,
                 ...,n
         - eps - the relative accuracy desired in the nodes
                 and weights

   Output:
         - zero- array of dimension  n  containing the Gaussian
                 nodes (in increasing order)  zero(k)=x(k), k=1,2,
                 ...,n
         - weight - array of dimension  n  containing the
                 Gaussian weights  weight(k)=w(k), k=1,2,...,n
         - ierr- an error flag equal to  0  on normal return,
                 equal to  i  if the QR algorithm does not
                 converge within 30 iterations on evaluating the
                 i-th eigenvalue, equal to  -1  if  n  is not in
                 range, and equal to  -2  if one of the beta's is
                 negative.

The array e is needed for working space.

Definition at line 388 of file sll_m_gauss_lobatto_integration.F90.

Here is the caller graph for this function:

◆ dlob()

subroutine sll_m_gauss_lobatto_integration::dlob (   n,
  dalpha,
  dbeta,
  dleft,
  dright,
  dzero,
  dweigh,
  ierr,
  de,
  da,
  db 
)
private

This comes from http://dl.acm.org, Algorithme 726 : ORTHPOL, appendices and supplements To use those functions, READ the documentation beside and find more information about coefficients in paper Algorithm 726 - ORTHPOL: A package of routines for generating orthogonal polynomials and Gauss-type quadrature rules by Walter Gautschi (here xxx is 726 in other references) formulas (1.1) to (1.3) page 2, and book Numerical Mathematics by Alfio Quarteroni, Riccardo Sacco and Fausto Saleri section 10.

Given n and a measure dlambda, this routine generates the (n+2)-point Gauss-Lobatto quadrature formula

integral over supp(dlambda) of f(x)dlambda(x)

= w(0)f(x(0)) + sum from k=1 to k=n of w(k)f(x(k))

    + w(n+1)f(x(n+1)) + R(n;f).

The nodes are returned as zero(k)=x(k), the weights as weight(k) =w(k), k=0,1,...,n,n+1. The user has to supply the recursion coefficients alpha(k), beta(k), k=0,1,...,n,n+1, for the measure dlambda. The nodes and weights are computed in terms of the eigenvalues and first component of the normalized eigenvectors of a slightly modified Jacobi matrix of order n+2. The routine calls upon the subroutine gauss and the function subroutine r1mach.

Input:

  • n - - the number of interior points in the Gauss-Lobatto formula; type integer
  • alpha,beta - arrays of dimension n+2 to be supplied with the recursion coefficients alpha(k-1), beta(k-1), k=1,2,...,n+2, of the underlying measure; the routine does not use alpha(n+2), beta(n+2)
  • aleft,right - the prescribed left and right endpoints x(0) and x(n+1) of the Gauss-Lobatto formula

Output:

  • zero - an array of dimension n+2 containing the nodes (in increasing order) zero(k)=x(k), k=0,1,...,n,n+1
  • weight-an array of dimension n+2 containing the weights weight(k)=w(k), k=0,1,...,n,n+1
  • ierr - an error flag inherited from the routine gauss

The arrays e,a,b are needed for working space.

Definition at line 312 of file sll_m_gauss_lobatto_integration.F90.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ gauss_lobatto_integral_1d()

real(kind=f64) function sll_m_gauss_lobatto_integration::gauss_lobatto_integral_1d ( procedure(function_1d)  f,
real(kind=f64), intent(in)  a,
real(kind=f64), intent(in)  b,
integer(kind=i32), intent(in)  n 
)
private

Gauss-Lobatto Quadrature.

To integrate the function \( f(x) \) (real-valued and of a single, real-valued argument x) over the interval \( [a,b] \), we use the Gauss-Lobatto formula

\[ \int_{-1}^1 f(x)dx \approx \sum_{k=1}^{n} w_k f(x_k) \]

where n represents the desired number of Gauss points.

the function maps the interval \( [-1,1] \) into the arbitrary interval \( [a,b] \).

To be considered is to split this function into degree-specific functions to avoid the select statement.

Parameters
fthe function to be integrated
[in]aleft-bound of the definition interval of f
[in]bright-bound of the definition interval of f
[in]nthe desired number of Gauss points
Returns
The value of the integral

Definition at line 61 of file sll_m_gauss_lobatto_integration.F90.

Here is the caller graph for this function:

◆ sll_f_gauss_lobatto_derivative_matrix()

function, public sll_m_gauss_lobatto_integration::sll_f_gauss_lobatto_derivative_matrix ( intent(in)  n,
dimension(n), optional  y 
)

Construction of the derivative matrix for Gauss-Lobatto, The matrix must be already allocated of size \( n^2 \).

\[ der(i,j)=(Phi_i.Phi'_j) \]

Definition at line 217 of file sll_m_gauss_lobatto_integration.F90.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ sll_f_gauss_lobatto_points()

real(kind=f64) function, dimension(n), public sll_m_gauss_lobatto_integration::sll_f_gauss_lobatto_points ( integer(kind=i32), intent(in)  n,
real(kind=f64), intent(in), optional  a,
real(kind=f64), intent(in), optional  b 
)

Returns a 1d array of size (n) containing gauss-lobatto points in the interval [a,b].

Parameters
[in]nNumber of gauss points.
[in]aOPTIONAL Minimum value of the interval.
[in]bOPTIONAL Maximun value of the interval.
Returns
xk array containing points

Definition at line 132 of file sll_m_gauss_lobatto_integration.F90.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ sll_f_gauss_lobatto_points_and_weights()

real(kind=f64) function, dimension(2, n), public sll_m_gauss_lobatto_integration::sll_f_gauss_lobatto_points_and_weights ( integer(kind=i32), intent(in)  n,
real(kind=f64), intent(in), optional  a,
real(kind=f64), intent(in), optional  b 
)

Returns a 2d array of size (2,n) containing gauss-lobatto points and weights in the interval [a,b].

Parameters
[in]nNumber of gauss points.
[in]aOPTIONAL Minimum value of the interval.
[in]bOPTIONAL Maximun value of the interval.
Returns
array containing points (1,:) and weights (2,:)
wx points and weights

Definition at line 115 of file sll_m_gauss_lobatto_integration.F90.

Here is the call graph for this function:

◆ sll_f_gauss_lobatto_weights()

real(kind=f64) function, dimension(n), public sll_m_gauss_lobatto_integration::sll_f_gauss_lobatto_weights ( integer(kind=i32), intent(in)  n,
real(kind=f64), intent(in), optional  a,
real(kind=f64), intent(in), optional  b 
)

Returns a 1d array of size (n) containing gauss-lobatto weights in the interval [a,b].

Parameters
[in]nNumber of gauss points.
[in]aOPTIONAL Minimum value of the interval.
[in]bOPTIONAL Maximun value of the interval.
Returns
wk array containing points

Definition at line 179 of file sll_m_gauss_lobatto_integration.F90.

Here is the call graph for this function:
Here is the caller graph for this function:
    Report Typos and Errors