From 72e3e9e3c2c8b6ee8df0aad093bc4c7c5753d9eb Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Wed, 24 Jan 2024 23:26:57 +0100 Subject: [PATCH 1/3] Add simpson's rule --- docs/src/basics/SampledIntegralProblem.md | 5 ++- src/Integrals.jl | 3 +- src/algorithms.jl | 26 +++++++++++- src/simpsons.jl | 51 +++++++++++++++++++++++ test/sampled_tests.jl | 7 +++- 5 files changed, 86 insertions(+), 6 deletions(-) create mode 100644 src/simpsons.jl diff --git a/docs/src/basics/SampledIntegralProblem.md b/docs/src/basics/SampledIntegralProblem.md index 02317cd4..ebc5ec4f 100644 --- a/docs/src/basics/SampledIntegralProblem.md +++ b/docs/src/basics/SampledIntegralProblem.md @@ -60,14 +60,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 ``` \ No newline at end of file diff --git a/src/Integrals.jl b/src/Integrals.jl index 02b5bc5d..0865f840 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -15,6 +15,7 @@ include("infinity_handling.jl") include("quadrules.jl") include("sampled.jl") include("trapezoidal.jl") +include("simpsons.jl") abstract type QuadSensitivityAlg end struct ReCallVJP{V} @@ -150,5 +151,5 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, lb, ub, p; SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success) end -export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre, QuadratureRule, TrapezoidalRule +export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre, QuadratureRule, TrapezoidalRule, SimpsonsRule end # module diff --git a/src/algorithms.jl b/src/algorithms.jl index 0adefdbd..153565ba 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -137,13 +137,37 @@ f = x -> x^2 x = range(0, 1, length=20) y = f.(x) problem = SampledIntegralProblem(y, x) -method = TrapezoidalRul() +method = TrapezoidalRule() 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) diff --git a/src/simpsons.jl b/src/simpsons.jl new file mode 100644 index 00000000..5da98b70 --- /dev/null +++ b/src/simpsons.jl @@ -0,0 +1,51 @@ +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 diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl index ef60ded8..35511c45 100644 --- a/test/sampled_tests.jl +++ b/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(grid2); 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 From 5beccfb9aafe977b1bd662ec304cb26b9f686161 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Fri, 26 Jan 2024 18:48:20 +0100 Subject: [PATCH 2/3] typo --- test/sampled_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl index 35511c45..b5e2d676 100644 --- a/test/sampled_tests.jl +++ b/test/sampled_tests.jl @@ -9,7 +9,7 @@ using Integrals, Test grid2 = [lb; sort(grid2); ub] grid3 = rand(npoints+1).*(ub-lb) .+ lb # also test odd number of points - grid3 = [lb; sort(grid2); ub] + grid3 = [lb; sort(grid3); ub] exact_sols = [1 / 6 * (ub^6 - lb^6), sin(ub) - sin(lb)] for method in [TrapezoidalRule, SimpsonsRule] From fd6ece1e66b3cd7183c8c1a5a500d97984628321 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Fri, 26 Jan 2024 19:08:27 +0100 Subject: [PATCH 3/3] format --- src/simpsons.jl | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/src/simpsons.jl b/src/simpsons.jl index 5da98b70..7c45b0f0 100644 --- a/src/simpsons.jl +++ b/src/simpsons.jl @@ -24,25 +24,36 @@ end 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." + @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])) + (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 * + 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 + 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])) + 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)