Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simpson's rule #222

Merged
merged 4 commits into from
Jan 27, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions docs/src/basics/SampledIntegralProblem.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
```
4 changes: 2 additions & 2 deletions src/Integrals.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -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
Expand Down
24 changes: 24 additions & 0 deletions src/algorithms.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
62 changes: 62 additions & 0 deletions src/simpsons.jl
Original file line number Diff line number Diff line change
@@ -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)

Check warning on line 21 in src/simpsons.jl

View check run for this annotation

Codecov / codecov/patch

src/simpsons.jl#L21

Added line #L21 was not covered by tests
# 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 *

Check warning on line 28 in src/simpsons.jl

View check run for this annotation

Codecov / codecov/patch

src/simpsons.jl#L24-L28

Added lines #L24 - L28 were not covered by tests
(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 *

Check warning on line 32 in src/simpsons.jl

View check run for this annotation

Codecov / codecov/patch

src/simpsons.jl#L31-L32

Added lines #L31 - L32 were not covered by tests
(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]) *

Check warning on line 35 in src/simpsons.jl

View check run for this annotation

Codecov / codecov/patch

src/simpsons.jl#L35

Added line #L35 was not covered by tests
(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 *

Check warning on line 38 in src/simpsons.jl

View check run for this annotation

Codecov / codecov/patch

src/simpsons.jl#L38

Added line #L38 was not covered by tests
(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 &&

Check warning on line 43 in src/simpsons.jl

View check run for this annotation

Codecov / codecov/patch

src/simpsons.jl#L43

Added line #L43 was not covered by tests
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) &&

Check warning on line 49 in src/simpsons.jl

View check run for this annotation

Codecov / codecov/patch

src/simpsons.jl#L49

Added line #L49 was not covered by tests
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) &&

Check warning on line 52 in src/simpsons.jl

View check run for this annotation

Codecov / codecov/patch

src/simpsons.jl#L52

Added line #L52 was not covered by tests
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)

Check warning on line 61 in src/simpsons.jl

View check run for this annotation

Codecov / codecov/patch

src/simpsons.jl#L61

Added line #L61 was not covered by tests
end
7 changes: 5 additions & 2 deletions test/sampled_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading