Merge pull request #222 from IlianPihlajamaa/Simpsons-rule

Simpson's rule

docs/src/basics/SampledIntegralProblem.md

@@ -61,14 +61,15 @@ f3 = x -> x^4
 x = range(0, 1, length=20)
 y = [f1.(x) f2.(x) f3.(x)]
 problem = SampledIntegralProblem(y, x; dim=1)
-method = TrapezoidalRule()
+method = SimpsonsRule()
 solve(problem, method)
 ```
 
 ### Supported methods
 
-Right now, only the `TrapezoidalRule` is supported, [see wikipedia](https://en.wikipedia.org/wiki/Trapezoidal_rule).
+Right now, only the [`TrapezoidalRule`](https://en.wikipedia.org/wiki/Trapezoidal_rule) and [`SimpsonsRule`](https://en.wikipedia.org/wiki/Simpson%27s_rule) are supported.
 
 ```@docs
 TrapezoidalRule
+SimpsonsRule
 ```

src/Integrals.jl

@@ -14,6 +14,7 @@ include("infinity_handling.jl")
 include("quadrules.jl")
 include("sampled.jl")
 include("trapezoidal.jl")
+include("simpsons.jl")
 
 abstract type QuadSensitivityAlg end
 struct ReCallVJP{V}
@@ -217,8 +218,7 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p;
     SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success)
 end
 
-
-export QuadGKJL, HCubatureJL, VEGAS, VEGASMC, GaussLegendre, QuadratureRule, TrapezoidalRule
+export QuadGKJL, HCubatureJL, VEGAS, VEGASMC, GaussLegendre, QuadratureRule, TrapezoidalRule, SimpsonsRule
 export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre
 export CubatureJLh, CubatureJLp
 export ArblibJL

src/algorithms.jl

@@ -148,6 +148,30 @@ solve(problem, method)
 struct TrapezoidalRule <: SciMLBase.AbstractIntegralAlgorithm
 end
 
+"""
+    SimpsonsRule
+
+Struct for evaluating an integral via the Simpson's composite 1/3-3/8
+rule over `AbstractRange`s (evenly spaced points) and 
+Simpson's composite 1/3 rule for non-equidistant grids.
+
+
+Example with equidistant data:
+
+```
+using Integrals 
+f = x -> x^2
+x = range(0, 1, length=20)
+y = f.(x)
+problem = SampledIntegralProblem(y, x)
+method = SimpsonsRule()
+solve(problem, method)
+```
+"""
+struct SimpsonsRule <: SciMLBase.AbstractIntegralAlgorithm
+end
+
+
 """
   QuadratureRule(q; n=250)
 

src/simpsons.jl

@@ -0,0 +1,62 @@
+struct SimpsonUniformWeights{T} <: UniformWeights
+    n::Int
+    h::T
+end
+
+@inline function Base.getindex(w::SimpsonUniformWeights, i)
+    # evenly spaced simpson's 1/3, 3/8 rule 
+    h = w.h
+    n = w.n
+    (i == 1 || i == n) && return 17h / 48
+    (i == 2 || i == n - 1) && return 59h / 48
+    (i == 3 || i == n - 2) && return 43h / 48
+    (i == 4 || i == n - 3) && return 49h / 48
+    return h
+end
+
+struct SimpsonNonuniformWeights{X <: AbstractArray} <: NonuniformWeights
+    x::X
+end
+
+@inline function Base.getindex(w::SimpsonNonuniformWeights, i)
+    # composite 1/3 rule for irregular grids
+
+    checkbounds(w.x, i)
+    x = w.x
+    j = i - firstindex(x)
+    @assert length(w.x)>2 "The length of the grid must exceed 2 for simpsons rule."
+    i == firstindex(x) && return (x[begin + 2] - x[begin + 0]) / 6 *
+           (2 - (x[begin + 2] - x[begin + 1]) / (x[begin + 1] - x[begin + 0]))
+
+    if isodd(length(x)) # even number of subintervals
+        i == lastindex(x) && return (x[end] - x[end - 2]) / 6 *
+               (2 - (x[end - 1] - x[end - 2]) / (x[end] - x[end - 1]))
+    else # odd number of subintervals, we add additional terms
+        i == lastindex(x) && return (x[end] - x[end - 1]) *
+               (2 * (x[end] - x[end - 1]) + 3 * (x[end - 1] - x[end - 2])) /
+               (x[end] - x[end - 2]) / 6
+        i == lastindex(x) - 1 && return (x[end - 1] - x[end - 3]) / 6 *
+               (2 - (x[end - 2] - x[end - 3]) / (x[end - 1] - x[end - 2])) +
+               (x[end] - x[end - 1]) *
+               ((x[end] - x[end - 1]) + 3 * (x[end - 1] - x[end - 2])) /
+               (x[end - 1] - x[end - 2]) / 6
+        i == lastindex(x) - 2 &&
+            return (x[end - 1] - x[end - 3])^3 / (x[end - 2] - x[end - 3]) /
+                   (x[end - 1] - x[end - 2]) / 6 -
+                   (x[end] - x[end - 1])^3 / (x[end - 1] - x[end - 2]) /
+                   (x[end] - x[end - 2]) / 6
+    end
+    isodd(j) &&
+        return (x[begin + j + 1] - x[begin + j - 1])^3 / (x[begin + j] - x[begin + j - 1]) /
+               (x[begin + j + 1] - x[begin + j]) / 6
+    iseven(j) && 
+        return (x[begin + j] - x[begin + j - 2]) / 6 *
+           (2 - (x[begin + j - 1] - x[begin + j - 2]) / (x[begin + j] - x[begin + j - 1])) +
+           (x[begin + j + 2] - x[begin + j]) / 6 *
+           (2 - (x[begin + j + 2] - x[begin + j + 1]) / (x[begin + j + 1] - x[begin + j]))
+end
+
+function find_weights(x::AbstractVector, ::SimpsonsRule)
+    x isa AbstractRange && return SimpsonUniformWeights(length(x), step(x))
+    return SimpsonNonuniformWeights(x)
+end

test/sampled_tests.jl

@@ -8,9 +8,12 @@ using Integrals, Test
     grid2 = rand(npoints) .* (ub - lb) .+ lb
     grid2 = [lb; sort(grid2); ub]
 
+    grid3 = rand(npoints+1).*(ub-lb) .+ lb # also test odd number of points
+    grid3 = [lb; sort(grid3); ub]
+
     exact_sols = [1 / 6 * (ub^6 - lb^6), sin(ub) - sin(lb)]
-    for method in [TrapezoidalRule] # Simpson's later
-        for grid in [grid1, grid2]
+    for method in [TrapezoidalRule, SimpsonsRule]
+        for grid in [grid1, grid2, grid3]
             for (i, f) in enumerate([x -> x^5, x -> cos(x)])
                 exact = exact_sols[i]
                 # single dimensional y