Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
|
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... | |
|
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.
|
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:
Output:
The arrays e,a,b are needed for working space.
Definition at line 312 of file sll_m_gauss_lobatto_integration.F90.
|
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.
f | the function to be integrated | |
[in] | a | left-bound of the definition interval of f |
[in] | b | right-bound of the definition interval of f |
[in] | n | the desired number of Gauss points |
Definition at line 61 of file sll_m_gauss_lobatto_integration.F90.
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.
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].
[in] | n | Number of gauss points. |
[in] | a | OPTIONAL Minimum value of the interval. |
[in] | b | OPTIONAL Maximun value of the interval. |
Definition at line 132 of file sll_m_gauss_lobatto_integration.F90.
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].
[in] | n | Number of gauss points. |
[in] | a | OPTIONAL Minimum value of the interval. |
[in] | b | OPTIONAL Maximun value of the interval. |
Definition at line 115 of file sll_m_gauss_lobatto_integration.F90.
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].
[in] | n | Number of gauss points. |
[in] | a | OPTIONAL Minimum value of the interval. |
[in] | b | OPTIONAL Maximun value of the interval. |
Definition at line 179 of file sll_m_gauss_lobatto_integration.F90.