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

Author: bernatguillen

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. 

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

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