+rombergs_method, ⬆️ ver->0.5.0

Author: fusion809

.gitignore

@@ -5,3 +5,4 @@
 *.pdf
 *.tex
 *.adoc
+*.ipynb

Project.toml

@@ -1,7 +1,7 @@
 name = "FunctionIntegrator"
 uuid = "7685536e-2581-4f83-bef1-2ba363c9cb91"
 authors = ["Brenton Horne <[email protected]>"]
-version = "0.4.1"
+version = "0.5.0"
 
 [deps]
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"

src/FunctionIntegrator.jl

@@ -2,7 +2,7 @@ module FunctionIntegrator
 
 using FastGaussQuadrature
 
-export chebyshev_quadrature, hermite_quadrature, jacobi_quadrature, laguerre_quadrature, legendre_quadrature, lobatto_quadrature, radau_quadrature, rectangle_rule_left, rectangle_rule_midpoint, rectangle_rule_right, simpsons_rule, simpsons38_rule, adaptive_simpsons_rule, trapezoidal_rule
+export chebyshev_quadrature, hermite_quadrature, jacobi_quadrature, laguerre_quadrature, legendre_quadrature, lobatto_quadrature, radau_quadrature, rectangle_rule_left, rectangle_rule_midpoint, rectangle_rule_right, rombergs_method, simpsons_rule, simpsons38_rule, adaptive_simpsons_rule, trapezoidal_rule
 
 include("Chebyshev.jl")
 include("Hermite.jl")
@@ -12,6 +12,7 @@ include("Legendre.jl")
 include("Lobatto.jl")
 include("Radau.jl")
 include("Rectangle.jl")
+include("Romberg.jl")
 include("Simpson.jl")
 include("Trapezoidal.jl")
 end # module

src/Romberg.jl

@@ -0,0 +1,28 @@
+"""
+    rombergs_method(f::Function, N::Number, a::Number, b::Number)
+
+computes the integral:
+
+``\\displaystyle \\int_a^b f(x) dx``
+
+using [Romberg's method](https://en.wikipedia.org/wiki/Romberg's_method#Method) with the ``n`` mentioned in the linked article being equal to ``N`` in the function arguments. 
+"""
+function rombergs_method(f::Function, N::Number, a::Number, b::Number)
+    N = convert(Int64, N);
+    n = 1:N;
+    h = (2.0 .^ (-n)).*(b-a);
+    n = nothing;
+    R = zeros(N+1,N+1);
+    # This is really R[0,0]
+    R[1,1] = h[1]*(f.(a)+f.(b));
+    for n=2:N+1
+        upper_bounds = convert(Int64, 2^(n-2));
+        k = 1:upper_bounds;
+        R[n,1] = 1/2*R[n-1,1] + h[n-1]*sum(f.((2*k.-1)*h[n-1].+a));
+        k = nothing;
+        for m=2:n
+            R[n,m] = (4.0^(m-1) - 1).^(-1)*((4^(m-1)*R[n,m-1]-R[n-1,m-1]));
+        end
+    end
+    return R[N+1,N+1]
+end

test/airyai.jl

@@ -27,6 +27,8 @@ printstyled("Performing the Airy Ai(x) test, where Ai(x) is integrated on the se
     @time @test rectangle_rule_midpoint(x -> airyai(x), 147348, 0, 100) ≈ 1.0/3.0
     printstyled("Running: rectangle_rule_right. Only a rough approximation can be practically achieved using this function.\n"; color = :magenta)
     @time @test abs(rectangle_rule_right(x -> airyai(x), 1e8, 0, 100) - 1.0/3.0) < 1.775176824944687e-7
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @test rombergs_method(x -> airyai(x), 9, 0, 100) ≈ 1.0/3.0
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(x -> airyai(x), 2512, 0, 100) ≈ 1.0/3.0
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/besselj.jl

@@ -29,6 +29,8 @@ printstyled("Performing the besselj test, where BesselJ_1(2) is approximated and
     @time @test rectangle_rule_midpoint(besselj_integrand, 4, 0, pi/2) ≈ besselj(1,2)
     printstyled("Running: rectangle_rule_right\n"; color = :magenta)
     @time @test rectangle_rule_right(besselj_integrand, 52995612, 0, pi/2) ≈ besselj(1,2)
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(besselj_integrand, 5, 0, pi/2) ≈ besselj(1,2)
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(besselj_integrand, 6, 0, pi/2) ≈ besselj(1,2)
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/cos_cot_integral.jl

@@ -31,6 +31,8 @@ printstyled("Integrating cos^2(x)/(1+cot(x)) from 0 to pi/2 and comparing the re
     @time @test rectangle_rule_midpoint(cos_cot_fn, 5254, 0, pi/2) ≈ 0.25
     printstyled("Running: rectangle_rule_right\n"; color = :magenta)
     @time @test rectangle_rule_right(cos_cot_fn, 7430, 0, pi/2) ≈ 0.25
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(cos_cot_fn, 4, 0, pi/2) ≈ 0.25
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(cos_cot_fn, 78, 0, pi/2) ≈ 0.25
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/cosine.jl

@@ -25,6 +25,8 @@ printstyled("Integrating cosine from 0 to pi/2 and comparing the result to the a
     @time @test rectangle_rule_midpoint(x -> cos(x), 2627, 0, pi/2) ≈ 1
     printstyled("Running: rectangle_rule_right\n"; color = :magenta)
     @time @test rectangle_rule_right(x -> cos(x), 5.2441522e7, 0, pi/2) ≈ 1
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(x -> cos(x), 3, 0, pi/2) ≈ 1
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(x -> cos(x), 40, 0, pi/2) ≈ 1
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/expnx2datan.jl

@@ -31,6 +31,8 @@ printstyled("Integrating e^(-x^2)/(1+x^2) on the infinite domain [-inf, inf], or
     @time @test rectangle_rule_midpoint(x -> exp(-x^2)/(x^2+1), 655, -100, 100) ≈ sol_9
     printstyled("Running: rectangle_rule_right\n"; color = :magenta)
     @time @test rectangle_rule_right(x -> exp(-x^2)/(x^2+1), 655, -100, 100) ≈ sol_9
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(x -> exp(-x^2)/(x^2+1), 12, -100, 100) ≈ sol_9
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(x -> exp(-x^2)/(x^2+1), 1240, -100, 100) ≈ sol_9
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/gaussian.jl

@@ -31,6 +31,8 @@ printstyled("Integrate exp(-x^2) from minus infinity to positive infinity and co
     @time @test rectangle_rule_midpoint(x -> exp(-x^2), 276, -100, 100) ≈ sqrt(pi)
     printstyled("Running: rectangle_rule_right\n"; color = :magenta)
     @time @test rectangle_rule_right(x -> exp(-x^2), 276, -100, 100) ≈ sqrt(pi)
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(x -> exp(-x^2), 12, -100, 100) ≈ sqrt(pi)
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(x -> exp(-x^2), 536, -100, 100) ≈ sqrt(pi)
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/logarithm.jl

@@ -25,6 +25,8 @@ printstyled("Integrating 1/x from 1 to e and comparing the result to the analyti
     @time @test rectangle_rule_midpoint(x -> x^(-1), 2672, 1, exp(1)) ≈ 1
     printstyled("Running: rectangle_rule_right\n"; color = :magenta)
     @time @test rectangle_rule_right(x -> x^(-1), 3.6607201e7, 1, exp(1)) ≈ 1
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(x -> 1/x, 5, 1, exp(1)) ≈ 1
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(x -> x^(-1), 70, 1, exp(1)) ≈ 1
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/logoverx.jl

@@ -24,6 +24,8 @@ printstyled("Integrating log(x)/x from 1 to e and comparing the result to the an
     @time @test rectangle_rule_midpoint(x -> log(x)/x, 4064, 1, exp(1)) ≈ 0.5
     printstyled("Running: rectangle_rule_right\n"; color = :magenta)
     @time @test rectangle_rule_right(x -> log(x)/x, 4.372011e7, 1, exp(1)) ≈ 0.5
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(x -> log(x)/x, 5, 1, exp(1)) ≈ 0.5
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(x -> log(x)/x, 100, 1, exp(1)) ≈ 0.5
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/modbessel0.jl

@@ -33,6 +33,8 @@ printstyled("Approximating the modified bessel function I_1(1) and comparing it
     @time @test rectangle_rule_midpoint(dmodified_bessel_1, 3, 0, pi/2) ≈ besseli(x,1)
     printstyled("Running: rectangle_rule_right\n"; color = :magenta)
     @time @test rectangle_rule_right(dmodified_bessel_1, 59275151, 0, pi/2) ≈ besseli(x,1)
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(dmodified_bessel_1, 5, 0, pi/2) ≈ besseli(x,1)
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(dmodified_bessel_1, 6, 0, pi/2) ≈ besseli(x,1)
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/simppen.jl

@@ -27,10 +27,12 @@ printstyled("Integrating 1/sqrt(-19.6 sin(x)) from -pi to 0 and comparing the re
     @time @test abs(rectangle_rule_midpoint(x -> (-19.6*sin(x))^(-0.5), 1e8, -pi+1e-8, -1e-8) - ellipk(1/2)/sqrt(2.45)) < 1.01e-4
     printstyled("Running: rectangle_rule_right on [-pi+1e-8, -1e-8]. Only a rough approximation can be realistically achieved with this function, due to the singularities.\n"; color = :magenta)
     @time @test abs(rectangle_rule_right(x -> (-19.6*sin(x))^(-0.5), 1e8, -pi+1e-8, -1e-8) - ellipk(1/2)/sqrt(2.45)) < 8.649e-5
+    printstyled("Running: rombergs_method on [-pi+1e-8, -1e-8]. Only a rough approximation can be realistically achieved with this function, due to the singularities.\n"; color = :magenta)
+    @time @test abs(rombergs_method(x -> (-19.6*sin(x))^(-0.5), 24, -pi+1e-8, -1e-8) - ellipk(1/2)/sqrt(2.45)) < 1e-4
     printstyled("Running: simpsons_rule on [-pi+1e-8, -1e-8]. Only a rough approximation can be realistically achieved with this function, due to the singularities.\n"; color = :magenta)
     @time @test abs(simpsons_rule(x -> (-19.6*sin(x))^(-0.5), 1e8, -pi+1e-8, -1e-8) - ellipk(1/2)/sqrt(2.45)) < 1e-4
     printstyled("Running: simpsons38_rule on [-pi+1e-8, -1e-8]. Only a rough approximation can be realistically achieved with this function, due to the singularities.\n"; color = :magenta)
     @time @test abs(simpsons38_rule(x -> (-19.6*sin(x))^(-0.5), 1e8+2, -pi+1e-8, -1e-8) - ellipk(1/2)/sqrt(2.45)) < 1e-4
     printstyled("Running: trapezoidal_rule on [-pi+1e-8, -1e-8]. Only a rough approximation can be realistically achieved with this function, due to the singularities.\n"; color = :magenta)
     @time @test abs(trapezoidal_rule(x -> (-19.6*sin(x))^(-0.5), 1e8, -pi+1e-8, -1e-8) - ellipk(1/2)/sqrt(2.45)) < 1e-4
-end
+end

test/sinexpox.jl

@@ -37,6 +37,8 @@ printstyled("Integrating sin(x^2)e^(-x)/x from 0 to infinity, with the approxima
     @time @test rectangle_rule_midpoint(sinexpox, 278981, 0, 100) ≈ sol_11
     printstyled("Running: rectangle_rule_right\n"; color = :magenta)
     @time @test rectangle_rule_right(sinexpox, 394538, 0, 100) ≈ sol_11
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(sinexpox, 12, 0, 100) ≈ sol_11
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(sinexpox, 4202, 0, 100) ≈ sol_11
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/sinxx.jl

@@ -31,6 +31,8 @@ printstyled("Integrating sin(x)/x from 0 to 100 and comparing it to the exact re
     @time @test rectangle_rule_midpoint(sinxx, 12461, 0, 100) ≈ sol_8
     printstyled("Running: rectangle_rule_right; only a rough approximation can be practically achieved using this function\n"; color = :magenta)
     @time @test abs(rectangle_rule_right(sinxx, 1e8, 0, 100) - sol_8) < 5.015e-7
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(sinxx, 9, 0, 100) ≈ sol_8
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(sinxx, 678, 0, 100) ≈ sol_8
     printstyled("Running: simpsons38_rule\n"; color = :magenta)

test/test_7.jl

@@ -30,6 +30,8 @@ printstyled("Integrating (x^3+1)/(x^4 (x+1)(x^2+1)) from 1 to e and comparing th
     @time @test rectangle_rule_midpoint(partfrac, 9887, 1, exp(1)) ≈ sol_7
     printstyled("Running: rectangle_rule_right; only a rough approximation can be practically achieved using this function.\n"; color = :magenta)
     @time @test abs(rectangle_rule_right(partfrac, 1e8, 1, exp(1)) - sol_7) < 1e-8
+    printstyled("Running: rombergs_method\n"; color = :magenta)
+    @time @test rombergs_method(partfrac, 6, 1, exp(1)) ≈ sol_7
     printstyled("Running: simpsons_rule\n"; color = :magenta)
     @time @test simpsons_rule(partfrac, 200, 1, exp(1)) ≈ sol_7
     printstyled("Running: simpsons38_rule\n"; color = :magenta)