From 982c5964b4411067cf76ee3a454b63170ffaedfb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bernat=20Guill=C3=A9n?= Date: Tue, 18 Sep 2018 11:52:44 -0400 Subject: [PATCH] Added Romberg's method (#14) * Added Romberg's method, no support for SVector yet * Changed RombergEven so that it is compatible with static vectors * Changed README to add Romberg * Added warning in case final step reached without reaching accuracy, added test for warning --- README.md | 5 +++-- src/NumericalIntegration.jl | 42 +++++++++++++++++++++++++++++++++++++ test/runtests.jl | 18 +++++++++++----- 3 files changed, 58 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index 87bcc24..b98c7f8 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ [![codecov.io](http://codecov.io/github/dextorious/NumericalIntegration.jl/coverage.svg?branch=master)](http://codecov.io/github/dextorious/NumericalIntegration.jl?branch=master) -This is a simple package to provide functionality for numerically integrating presampled data (meaning you can't choose arbitrary nodes). If you have the ability to evaluate your integrand at arbitrary points, please consider using better tools for the job (such as the excellent [FastGaussQuadrature.jl](https://github.com/ajt60gaibb/FastGaussQuadrature.jl)). +This is a simple package to provide functionality for numerically integrating presampled data (meaning you can't choose arbitrary nodes). If you have the ability to evaluate your integrand at arbitrary points, please consider using better tools for the job (such as the excellent [FastGaussQuadrature.jl](https://github.com/ajt60gaibb/FastGaussQuadrature.jl)). Do note that while the code is trivial, it has not been extensively tested and does not focus on numerical precision. Issues, suggestions and pull requests are welcome. @@ -32,5 +32,6 @@ The currently available methods are: - TrapezoidalEvenFast - SimpsonEven - SimpsonEvenFast +- RombergEven -All methods containing "Even" in the name assume evenly spaced data. All methods containing "Fast" omit basic correctness checks and focus on performance. Consequently, the fast methods will segfault or produce incorrect results if you supply incorrect data (vectors of different lengths, etc.). +All methods containing "Even" in the name assume evenly spaced data. All methods containing "Fast" omit basic correctness checks and focus on performance. Consequently, the fast methods will segfault or produce incorrect results if you supply incorrect data (vectors of different lengths, etc.). RombergEven needs a power of 2 + 1 points (so 9, 17, 33, 65, 129, 257, 513, 1025...) evenly spaced for it to work. Useful when control over accuracy is needed. diff --git a/src/NumericalIntegration.jl b/src/NumericalIntegration.jl index 5f0b820..b368ecf 100644 --- a/src/NumericalIntegration.jl +++ b/src/NumericalIntegration.jl @@ -1,8 +1,12 @@ module NumericalIntegration +using LinearAlgebra +using Logging + export integrate export Trapezoidal, TrapezoidalEven, TrapezoidalFast, TrapezoidalEvenFast export SimpsonEven, SimpsonEvenFast +export RombergEven export IntegrationMethod abstract type IntegrationMethod end @@ -13,6 +17,10 @@ struct TrapezoidalFast <: IntegrationMethod end struct TrapezoidalEvenFast <: IntegrationMethod end struct SimpsonEven <: IntegrationMethod end # https://en.wikipedia.org/wiki/Simpson%27s_rule#Alternative_extended_Simpson.27s_rule struct SimpsonEvenFast <: IntegrationMethod end +struct RombergEven{T<:AbstractFloat} <: IntegrationMethod + acc::T +end # https://en.wikipedia.org/wiki/Romberg%27s_method +RombergEven() = RombergEven(1e-12) const HALF = 1//2 @@ -71,6 +79,40 @@ function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEvenFast) @inbounds return (x[2] - x[1]) * retval end +function integrate(x::AbstractVector, y::AbstractVector, m::RombergEven) + @assert length(x) == length(y) "x and y vectors must be of the same length!" + @assert ((length(x) - 1) & (length(x) - 2)) == 0 "Need length of vector to be 2^n + 1" + maxsteps::Integer = Int(log2(length(x)-1)) + rombaux = zeros(eltype(y), maxsteps, 2) + prevrow = 1 + currrow = 2 + @inbounds h = x[end] - x[1] + @inbounds rombaux[prevrow, 1] = (y[1] + y[end])*h*HALF + @inbounds for i in 1 : (maxsteps-1) + h *= HALF + npoints = 1 << (i-1) + jumpsize = div(length(x)-1, 2*npoints) + c = 0.0 + for j in 1 : npoints + c += y[1 + (2*j-1)*jumpsize] + end + rombaux[1, currrow] = h*c + HALF*rombaux[1, prevrow] + for j in 2 : (i+1) + n_k = 4^(j-1) + rombaux[j, currrow] = (n_k*rombaux[j-1, currrow] - rombaux[j-1, prevrow])/(n_k - 1) + end + + if i > maxsteps//3 && norm(rombaux[i, prevrow] - rombaux[i+1, currrow], Inf) < m.acc + return rombaux[i+1, currrow] + end + + prevrow, currrow = currrow, prevrow + end + finalerr = norm(rombaux[maxsteps-1, prevrow] - rombaux[maxsteps, currrow], Inf) + @warn "RombergEven :: final step reached, but accuracy not: $finalerr > $(m.acc)" + @inbounds return rombaux[maxsteps, prevrow] +end + integrate(x::AbstractVector, y::AbstractVector) = integrate(x, y, TrapezoidalFast()) end diff --git a/test/runtests.jl b/test/runtests.jl index 7ab3fa6..d8aa826 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,9 +4,9 @@ using InteractiveUtils # for subtypes using StaticArrays @testset "compare with analytic result" begin - x = collect(-π : π/1000 : π) + x = collect(-π : 2*π/2048 : π) y = sin.(x) - p = collect(0 : π/1000 : π) + p = collect(0 : π/1024 : π) q = sin.(p) for M in subtypes(IntegrationMethod) for T in [Float32, Float64, BigFloat] @@ -23,9 +23,9 @@ using StaticArrays end @testset "SVector" begin - xs = range(0, stop=1, length=10) - ys1 = randn(10) - ys2 = randn(10) + xs = range(0, stop=1, length=9) + ys1 = randn(9) + ys2 = randn(9) ys = map(SVector, ys1, ys2) for M in subtypes(IntegrationMethod) m = M() @@ -35,3 +35,11 @@ end @test res ≈ @SVector [res1, res2] end end + +@testset "Raising Warnings" begin + xs = collect(-1.0 : 0.5 : 1.0) + ys = xs.^2 + m = RombergEven() + expwarn = "RombergEven :: final step reached, but accuracy not: 1.0 > 1.0e-12" + @test_logs (:warn, expwarn) integrate(xs, ys, m) +end