Merge pull request #264 from kishore-nori/array-gl-fix

Small fix for composite Gauss-Legendre quadrature with array-valued integrands

Author: ChrisRackauckas

ext/IntegralsFastGaussQuadratureExt.jl

@@ -26,9 +26,9 @@ function gauss_legendre(f::F, p, lb, ub, nodes, weights) where {F}
 end
 function composite_gauss_legendre(f::F, p, lb, ub, nodes, weights, subintervals) where {F}
     h = (ub - lb) / subintervals
-    I = zero(h)
-    for i in 1:subintervals
-        _lb = lb + (i - 1) * h
+    I = gauss_legendre(f, p, lb, lb + h, nodes, weights)
+    for i in 1:subintervals-1
+        _lb = lb + i * h
         _ub = _lb + h
         I += gauss_legendre(f, p, _lb, _ub, nodes, weights)
     end

test/derivative_tests.jl

@@ -15,6 +15,8 @@ alg_req = Dict(
         allows_iip = false),
     GaussLegendre(n = 50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
         allows_iip = false),
+    GaussLegendre(n = 50, subintervals=3) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
+        allows_iip = false),
     QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1,
         allows_iip = true),
     HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1,

test/gaussian_quadrature_tests.jl

@@ -97,3 +97,21 @@ sol = solve(prob, alg)
 alg = GaussLegendre(n = 500, subintervals = 17)
 sol = solve(prob, alg)
 @test sol.u ≈ -240.25235266303063249920743158729
+
+# Gauss-Legendre with array-valued integrands
+prob = IntegralProblem((x,p)->[1,2], (-1,1))
+sol = solve(prob, GaussLegendre(;n=5))
+@test sol.u ≈ 2*[1,2]
+
+prob = IntegralProblem((x,p)->[1 2; 3 4], (-1,1))
+sol = solve(prob, GaussLegendre(;n=5))
+@test sol.u ≈ 2*[1 2; 3 4]
+
+# Composite Gauss-Legendre with array-valued integrands
+prob = IntegralProblem((x,p)->[1,2], (-1,1))
+sol = solve(prob, GaussLegendre(;n=5, subintervals=5))
+@test sol.u ≈ 2*[1,2]
+
+prob = IntegralProblem((x,p)->[1 2; 3 4], (-1,1))
+sol = solve(prob, GaussLegendre(;n=5, subintervals=5))
+@test sol.u ≈ 2*[1 2; 3 4]