legendre polynomials and gaussian quadrature

doc/specs/stdlib_quadrature.md

@@ -186,3 +186,95 @@ program demo_simps_weights
 ! 64.0
 end program demo_simps_weights
 ```
+
+## `gauss_legendre` - Gauss-Legendre quadrature (a.k.a. Gaussian quadrature) nodes and weights
+
+### Status
+
+Experimental
+
+### Description
+
+Computes Gauss-Legendre quadrature (also known as simply Gaussian quadrature) nodes and weights,
+ for any `N` (number of nodes).
+Using the nodes `x` and weights `w`, you can compute the integral of some function `f` as follows:
+`integral = sum(f(x) * w)`.
+
+Only double precision is supported - if lower precision is required, you must do the appropriate conversion yourself.
+Accuracy has been validated up to N=64 by comparing computed results to tablulated values known to be accurate to machine precision
+(maximum difference from those values is 2 epsilon).
+
+### Syntax
+
+`subroutine [[stdlib_quadrature(module):gauss_legendre(interface)]] (x, w[, interval])`
+
+### Arguments
+
+`x`: Shall be a rank-one array of type `real(real64)`. It is an *output* argument, representing the quadrature nodes.
+
+`w`: Shall be a rank-one array of type `real(real64)`, with the same dimension as `x`. 
+It is an *output* argument, representing the quadrature weights.
+
+`interval`: (Optional) Shall be a two-element array of type `real(real64)`. 
+If present, the nodes and weigts are calculated for integration from `interval(1)` to `interval(2)`.
+If not specified, the default integral is -1 to 1.
+
+### Example
+
+```fortran
+program integrate
+	use iso_fortran_env, dp => real64
+	implicit none
+
+	integer, parameter :: N = 6
+	real(dp), dimension(N) :: x,w
+	call gauss_legendre(x,w)
+	print *, "integral of x**2 from -1 to 1 is",  sum(x**2 * w)
+end program
+```
+
+## `gauss_legendre_lobatto` - Gauss-Legendre-Lobatto quadrature nodes and weights
+
+### Status
+
+Experimental
+
+### Description
+
+Computes Gauss-Legendre-Lobatto quadrature nodes and weights,
+ for any `N` (number of nodes).
+Using the nodes `x` and weights `w`, you can compute the integral of some function `f` as follows:
+`integral = sum(f(x) * w)`.
+
+Only double precision is supported - if lower precision is required, you must do the appropriate conversion yourself.
+Accuracy has been validated up to N=64 by comparing computed results to tablulated values known to be accurate to machine precision
+(maximum difference from those values is 2 epsilon).
+
+### Syntax
+
+`subroutine [[stdlib_quadrature(module):gauss_legendre_lobatto(interface)]] (x, w[, interval])`
+
+### Arguments
+
+`x`: Shall be a rank-one array of type `real(real64)`. It is an *output* argument, representing the quadrature nodes.
+
+`w`: Shall be a rank-one array of type `real(real64)`, with the same dimension as `x`. 
+It is an *output* argument, representing the quadrature weights.
+
+`interval`: (Optional) Shall be a two-element array of type `real(real64)`. 
+If present, the nodes and weigts are calculated for integration from `interval(1)` to `interval(2)`.
+If not specified, the default integral is -1 to 1.
+
+### Example
+
+```fortran
+program integrate
+	use iso_fortran_env, dp => real64
+	implicit none
+
+	integer, parameter :: N = 6
+	real(dp), dimension(N) :: x,w
+	call gauss_legendre_lobatto(x,w)
+	print *, "integral of x**2 from -1 to 1 is",  sum(x**2 * w)
+end program
+```

doc/specs/stdlib_specialfunctions.md

@@ -0,0 +1,63 @@
+---
+title: specialfunctions
+---
+
+# Special functions
+
+[TOC]
+
+## `legendre` - Calculate Legendre polynomials
+
+### Status
+
+Experimental
+
+### Description
+
+Computes the value of the n-th Legendre polynomial at a specified point.
+Currently only 64 bit floating point is supported.
+
+This is an `elemental` function.
+
+### Syntax
+
+`result = [[stdlib_specialfunctions(module):legendre(interface)]] (n, x)`
+
+### Arguments
+
+`n`: Shall be a scalar of type `real(real64)`. 
+
+`x`: Shall be a scalar or array (this function is elemental) of type `real(real64)`. 
+
+### Return value
+
+The function result will be the value of the `n`-th Legendre polynomial, evaluated at `x`.
+
+
+
+## `dlegendre` - Calculate first derivatives of Legendre polynomials
+
+### Status
+
+Experimental
+
+### Description
+
+Computes the value of the first derivative of the n-th Legendre polynomial at a specified point.
+Currently only 64 bit floating point is supported.
+
+This is an `elemental` function.
+
+### Syntax
+
+`result = [[stdlib_specialfunctions(module):dlegendre(interface)]] (n, x)`
+
+### Arguments
+
+`n`: Shall be a scalar of type `real(real64)`. 
+
+`x`: Shall be a scalar or array (this function is elemental) of type `real(real64)`. 
+
+### Return value
+
+The function result will be the value of the first derivative of the `n`-th Legendre polynomial, evaluated at `x`.

src/CMakeLists.txt

@@ -45,6 +45,9 @@ set(SRC
     stdlib_logger.f90
     stdlib_strings.f90
     stdlib_system.F90
+    stdlib_specialfunctions.f90
+    stdlib_specialfunctions_legendre.f90
+    stdlib_quadrature_gauss.f90
     ${outFiles}
 )
 

src/Makefile.manual

@@ -25,8 +25,12 @@ SRCFYPP =\
 
 SRC = f18estop.f90 \
       stdlib_error.f90 \
+      stdlib_specialfunctions.f90 \
+      stdlib_specialfunctions_legendre.f90 \
+      stdlib_io.f90 \
       stdlib_kinds.f90 \
       stdlib_logger.f90 \
+      stdlib_quadrature_gauss.f90 \
       stdlib_strings.f90 \
       $(SRCGEN)
 
@@ -61,6 +65,8 @@ stdlib_bitsets.o: stdlib_kinds.o
 stdlib_bitsets_64.o: stdlib_bitsets.o
 stdlib_bitsets_large.o: stdlib_bitsets.o
 stdlib_error.o: stdlib_optval.o
+stdlib_specialfunctions.o: stdlib_kinds.o
+stdlib_specialfunctions_legendre.o: stdlib_kinds.o stdlib_specialfunctions.o
 stdlib_io.o: \
 	stdlib_error.o \
 	stdlib_optval.o \
@@ -73,6 +79,9 @@ stdlib_linalg_diag.o: \
 stdlib_logger.o: stdlib_ascii.o stdlib_optval.o
 stdlib_optval.o: stdlib_kinds.o
 stdlib_quadrature.o: stdlib_kinds.o
+
+stdlib_quadrature_gauss.o: stdlib_kinds.o stdlib_quadrature.o
+
 stdlib_quadrature_simps.o: \
 	stdlib_quadrature.o \
 	stdlib_error.o \

src/stdlib_quadrature.fypp

@@ -12,6 +12,8 @@ module stdlib_quadrature
     public :: trapz_weights
     public :: simps
     public :: simps_weights
+    public :: gauss_legendre
+    public :: gauss_legendre_lobatto
 
 
     interface trapz
@@ -90,6 +92,28 @@ module stdlib_quadrature
     end interface simps_weights
 
 
+    interface gauss_legendre
+        !! version: experimental
+        !!
+        !! Computes Gauss-Legendre quadrature nodes and weights.
+        pure module subroutine gauss_legendre_fp64 (x, w, interval)
+            real(dp), intent(out) :: x(:), w(:)
+            real(dp), intent(in), optional :: interval(2)
+        end subroutine
+    end interface gauss_legendre
+
+
+    interface gauss_legendre_lobatto
+        !! version: experimental
+        !!
+        !! Computes Gauss-Legendre-Lobatto quadrature nodes and weights.
+        pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
+            real(dp), intent(out) :: x(:), w(:)
+            real(dp), intent(in), optional :: interval(2)
+        end subroutine
+    end interface gauss_legendre_lobatto
+
+
     ! Interface for a simple f(x)-style integrand function.
     ! Could become fancier as we learn about the performance
     ! ramifications of different ways to do callbacks.
@@ -103,4 +127,6 @@ module stdlib_quadrature
         #:endfor
     end interface
 
+
+
 end module stdlib_quadrature

src/stdlib_quadrature_gauss.f90

@@ -0,0 +1,122 @@
+submodule (stdlib_quadrature) stdlib_quadrature_gauss
+    use stdlib_specialfunctions, only: legendre, dlegendre
+    implicit none
+
+    real(dp), parameter :: pi = acos(-1._dp)
+    real(dp), parameter :: tolerance = 4._dp * epsilon(1._dp)
+    integer, parameter :: newton_iters = 100
+
+contains
+
+    pure module subroutine gauss_legendre_fp64 (x, w, interval)
+        real(dp), intent(out) :: x(:), w(:)
+        real(dp), intent(in), optional :: interval(2)
+
+        associate (n => size(x)-1 )
+        select case (n)
+            case (0)
+                x = 0
+                w = 2
+            case (1)
+                x(1) = -sqrt(1._dp/3._dp)
+                x(2) = -x(1)
+                w = 1
+            case default
+                block
+                integer :: i,j
+                real(dp) :: leg, dleg, delta
+
+                do i = 0, (n+1)/2 - 1
+                    ! use Gauss-Chebyshev points as an initial guess
+                    x(i+1) = -cos((2*i+1)/(2._dp*n+2._dp) * pi)
+                    do j = 1, newton_iters
+                        leg  = legendre(n+1,x(i+1))
+                        dleg = dlegendre(n+1,x(i+1))
+                        delta = -leg/dleg
+                        x(i+1) = x(i+1) + delta
+                        if ( abs(delta) <= tolerance * abs(x(i+1)) )  exit
+                    end do
+                    x(n-i+1) = -x(i+1)
+
+                    dleg = dlegendre(n+1,x(i+1))
+                    w(i+1)   = 2._dp/((1-x(i+1)**2)*dleg**2) 
+                    w(n-i+1) = w(i+1)
+                end do
+
+                if (mod(n,2) == 0) then
+                    x(n/2+1) = 0
+
+                    dleg = dlegendre(n+1, 0.0_dp)
+                    w(n/2+1) = 2._dp/(dleg**2) 
+                end if
+                end block
+        end select
+        end associate
+
+        if (present(interval)) then
+            associate ( a => interval(1) , b => interval(2) )
+            x = 0.5_dp*(b-a)*x+0.5_dp*(b+a)
+            x(1)       = interval(1)
+            x(size(x)) = interval(2)
+            w = 0.5_dp*(b-a)*w
+            end associate
+        end if
+    end subroutine
+
+    pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
+        real(dp), intent(out) :: x(:), w(:)
+        real(dp), intent(in), optional :: interval(2)
+
+        associate (n => size(x)-1)
+        select case (n)
+            case (1)
+                x(1) = -1
+                x(2) =  1
+                w = 1
+            case default
+                block
+                integer :: i,j
+                real(dp) :: leg, dleg, delta
+
+                x(1)   = -1._dp
+                x(n+1) =  1._dp
+                w(1)   =  2._dp/(n*(n+1._dp))
+                w(n+1) =  2._dp/(n*(n+1._dp))
+
+                do i = 1, (n+1)/2 - 1
+                    ! initial guess from an approximate form given by SV Parter (1999)
+                    x(i+1) = -cos( (i+0.25_dp)*pi/n  - 3/(8*n*pi*(i+0.25_dp)))
+                    do j = 1, newton_iters
+                        leg  = legendre(n+1,x(i+1)) - legendre(n-1,x(i+1))
+                        dleg = dlegendre(n+1,x(i+1)) - dlegendre(n-1,x(i+1))
+                        delta = -leg/dleg
+                        x(i+1) = x(i+1) + delta
+                        if ( abs(delta) <= tolerance * abs(x(i+1)) )  exit
+                    end do
+                    x(n-i+1) = -x(i+1)
+
+                    leg = legendre(n, x(i+1))
+                    w(i+1)   = 2._dp/(n*(n+1._dp)*leg**2) 
+                    w(n-i+1) = w(i+1)
+                end do
+
+                if (mod(n,2) == 0) then
+                    x(n/2+1) = 0
+
+                    leg = legendre(n, 0.0_dp)
+                    w(n/2+1)   = 2._dp/(n*(n+1._dp)*leg**2) 
+                end if
+                end block
+        end select
+        end associate
+
+        if (present(interval)) then
+            associate ( a => interval(1) , b => interval(2) )
+            x = 0.5_dp*(b-a)*x+0.5_dp*(b+a)
+            x(1)       = interval(1)
+            x(size(x)) = interval(2)
+            w = 0.5_dp*(b-a)*w
+            end associate
+        end if
+    end subroutine
+end submodule