Added Legendre polynomials

README.md

@@ -82,5 +82,5 @@ combination with DifferentialEquations.jl) see
 
 ## License
 
-Written by David A.W. Barton ([email protected]) 2016-2018 and released
+Written by David A.W. Barton ([email protected]) 2016-2021 and released
 under the MIT license <https://opensource.org/licenses/MIT>.

src/BarycentricInterpolation.jl

@@ -7,31 +7,36 @@ equispaced points and Chebyshev points of the first and second kind. The
 formulae used are taken from the paper of Berrut and Trefethen, SIAM Review,
 2004.
 
+Additional Barycentric weights for Legendre points are taken from the paper of Wang,
+Huybrechs, and Vandewalle, Mathematics of Computation, 83 (290) 2893-2914, 2014.
+
 Other packages that may be of interest are FastGaussQuadrature and ApproxFun.
 
-Written by David A.W. Barton ([email protected]) 2016-2018 and licensed
+Written by David A.W. Barton ([email protected]) 2016-2021 and licensed
 under the MIT license <https://opensource.org/licenses/MIT>
 """
 module BarycentricInterpolation
 
-export Chebyshev1, Chebyshev2, Equispaced, ArbitraryPolynomial,
-    weights, nodes, interpolate, interpolation_matrix, differentiation_matrix, degree
+using FastGaussQuadrature: gausslegendre
+
+export Chebyshev1, Chebyshev2, Legendre, Equispaced, ArbitraryPolynomial
+
+export weights, nodes, interpolate, interpolation_matrix, differentiation_matrix, degree
 
 #-- Node distributions
 
 abstract type AbstractPolynomial{N, T <: Number} end
 
-for name in [:Chebyshev1, :Chebyshev2, :Equispaced]
+for name in [:Chebyshev1, :Chebyshev2, :Legendre, :Equispaced]
     @eval struct $name{N, T <: Number} <: AbstractPolynomial{N, T}
         shift::T
         scale::T
         nodes::Vector{T}
         weights::Vector{T}
         function $name{N, T}(start::T, stop::T) where {N, T <: Number}
-            _shift = (stop + start)/2
-            _scale = (stop - start)/2
-            _nodes = nodes($name{N, T}, _shift, _scale)
-            _weights = weights($name{N, T})
+            _shift = T(stop + start)/2
+            _scale = T(stop - start)/2
+            (_nodes, _weights) = nodes_weights($name{N, T}, _shift, _scale)
             new{N, T}(_shift, _scale, _nodes, _weights)
         end
     end
@@ -52,14 +57,31 @@ end
 (poly::ArbitraryPolynomial)(y) = interpolate(poly, y)
 (poly::ArbitraryPolynomial)(y, x) = interpolate(poly, y, x)
 
-
 """
     degree(poly)
 
 Return the degree of the polynomial specified.
 """
 degree(poly::AbstractPolynomial{N}) where N = N
 
+"""
+    nodes_weights(poly)
+
+Return the nodes and weights of the polynomial specified.
+"""
+function nodes_weights(poly::Type{<:AbstractPolynomial}, shift=0, scale=1)
+    return (nodes(poly, shift, scale), weights(poly))
+end
+
+nodes_weights(poly::AbstractPolynomial) = (poly.nodes, poly.weights)
+
+function nodes_weights(::Type{<:Legendre{N, Float64}}, shift=0, scale=1) where N
+    (x, w) = gausslegendre(N+1)
+    _nodes = [xᵢ * scale + shift for xᵢ in x]
+    _weights = [(2*isodd(i)-1)*sqrt((1-x[i]^2)*w[i]) for i in eachindex(x)]
+    return (_nodes, _weights)
+end
+
 """
     weights(poly)
 
@@ -94,6 +116,8 @@ function weights(poly::Type{<: ArbitraryPolynomial{N, T}}, x::AbstractVector{T})
     return w
 end
 
+function weights(poly::Type{<:Legendre}, x::AbstractVector{T}) where T end
+
 weights(poly::AbstractPolynomial) = poly.weights
 
 """
@@ -223,4 +247,4 @@ function differentiation_matrix(poly::AbstractPolynomial{N, T}) where {N, T}
     return D
 end
 
-end
+end  # module

test/runtests.jl

@@ -11,6 +11,8 @@ cheb1₅ = Chebyshev1{5}()
 cheb1₆ = Chebyshev1{6}()
 cheb2₅ = Chebyshev2{5}()
 cheb2₆ = Chebyshev2{6}()
+leg₅ = Legendre{5}()
+leg₆ = Legendre{6}()
 
 x5 = [0.814723686393179,1.7205156234688,1.8475024397623,2.76087829590132,3.39323754212673,3.49077794712614]
 x6 = [0.278498218867048,0.825379738072032,1.78288657350633,2.74777510870561,2.90538819038315,3.87598097214377,4.83314792038672]
@@ -32,6 +34,9 @@ end
     @test nodes(cheb2₅) ≈ -[1,0.809016994374947,0.309016994374947,-0.309016994374947,-0.809016994374947,-1]
     @test nodes(cheb2₆) ≈ -[1,0.866025403784439,0.5,6.12323399573677e-17,-0.5,-0.866025403784439,-1]
 
+    @test nodes(leg₅) ≈ [-0.932469514203152,-0.661209386466264,-0.238619186083197,0.238619186083197,0.661209386466264,0.932469514203152]
+    @test nodes(leg₆) ≈ [-0.949107912342758,-0.741531185599394,-0.405845151377397,0,0.405845151377397,0.741531185599394,0.949107912342758]
+
     @test nodes(arb₅) ≈ x5
     @test nodes(arb₆) ≈ x6
 
@@ -49,6 +54,15 @@ end
     @test weights(cheb2₅) ≈ -[0.5,-1,1,-1,1,-0.5]
     @test weights(cheb2₆) ≈ -[0.5,-1,1,-1,1,-1,0.5]
 
+    leg₅w_arb = weights(ArbitraryPolynomial(nodes(leg₅)))
+    leg₅w = weights(leg₅)
+    leg₅w .*= leg₅w_arb[3] / leg₅w[3]
+    @test leg₅w_arb ≈ leg₅w
+    leg₆w_arb = weights(ArbitraryPolynomial(nodes(leg₆)))
+    leg₆w = weights(leg₆)
+    leg₆w .*= leg₆w_arb[3] / leg₆w[3]
+    @test leg₆w_arb ≈ leg₆w
+
     @test weights(arb₅) ≈ [-0.079601734468871,2.82206348049722,-3.28654127480991,1.17155138701322,-2.4317744194857,1.80430256125404]
     @test weights(arb₆) ≈ [0.0114361913589559,-0.0390633084690722,0.100394094241955,-0.588770483980979,0.552871013562898,-0.0415349717143145,0.00466746500055649]