Improving efficiency of algorithms for rect, simps and trapez

Author: fusion809

.travis.yml

@@ -3,6 +3,5 @@ julia:
   - nightly
   - 1.4
   - 1.5
-  - 1.6
 os:
   - linux

Project.toml

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

src/Rectangle.jl

@@ -13,7 +13,7 @@ function rectangle_rule_left(f::Function, N::Number, a::Number, b::Number)
     y = 0;
     x = a;
     for i=1:N
-        y += h*f(x)
+        y += h*f(x);
         x += h;
     end
     return y
@@ -34,7 +34,7 @@ function rectangle_rule_midpoint(f::Function, N::Number, a::Number, b::Number)
     y = 0;
     x = a+h/2;
     for i=1:N
-        y += h*f(x)
+        y += h*f(x);
         x += h;
     end
     return y
@@ -55,7 +55,7 @@ function rectangle_rule_right(f::Function, N::Number, a::Number, b::Number)
     y = 0;
     x = a+h;
     for i=1:N
-        y += h*f(x)
+        y += h*f(x);
         x += h;
     end
     return y

src/Simpson.jl

@@ -1,11 +1,5 @@
-function stepwise_simpsons(f::Function, h::Number, x::Number, i::Integer, N::Number)
-    if i == 1 || i == N + 1
-        return h / 3 * f(x)
-    elseif (i % 2) == 1
-        return 2 * h / 3 * f(x)
-    else
-        return 4 * h / 3 * f(x)
-    end
+function stepwise_simpsons(f::Function, h::Number, x::Number)
+    return h/6 * (f(x) + 4*f(x+h/2) + f(x+h));
 end
 
 function stepwise_simpsons38(f::Function, h::Number, x::Number, i::Integer, N::Number)
@@ -28,15 +22,12 @@ uses [Simpson's rule](https://en.wikipedia.org/wiki/Simpson%27s_rule) to approxi
 """
 function simpsons_rule(f::Function, N::Number, a::Number, b::Number)
     N = convert(Int64, N);
-    iseven(N) || error("N must be even in order for Simspon's rule to work properly.")
     h = (b-a)/N;
     y = 0;
     x = a;
-    for i=1:N+1
-        y = y + stepwise_simpsons(f, h, x, i, N);
-        if i < N+1
-            x += h;
-        end
+    for i=1:N
+        y += stepwise_simpsons(f, h, x);
+        x += h;
     end
     return y
 end

src/Trapezoidal.jl

@@ -1,9 +1,5 @@
 function stepwise_trapezoidal(f, h, x, i, N)
-    if i == 1 || i == N + 1
-        return h / 2 * f(x)
-    else
-        return h * f(x)
-    end
+    return h/2 * (f(x) + f(x+h));
 end
 
 """
@@ -20,11 +16,9 @@ function trapezoidal_rule(f::Function, N::Number, a::Number, b::Number)
     h = (b-a)/N;
     y = 0;
     x = a;
-    for i=1:N+1
-        y = y + stepwise_trapezoidal(f, h, x, i, N);
-        if i < N+1
-            x += h;
-        end
+    for i=1:N
+        y += stepwise_trapezoidal(f, h, x, i, N);
+        x += h;
     end
     return y
 end

test/airyai.jl

@@ -30,7 +30,7 @@ printstyled("Performing the Airy Ai(x) test, where Ai(x) is integrated on the se
     printstyled("Running: rombergs_method\n"; color = :magenta)
     @time @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
+    @time @test simpsons_rule(x -> airyai(x), 1256, 0, 100) ≈ 1.0/3.0
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(x -> airyai(x), 3075, 0, 100) ≈ 1.0/3.0
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/besselj.jl

@@ -32,7 +32,7 @@ printstyled("Performing the besselj test, where BesselJ_1(2) is approximated and
     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)
+    @time @test simpsons_rule(besselj_integrand, 3, 0, pi/2) ≈ besselj(1,2)
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(besselj_integrand, 9, 0, pi/2) ≈ besselj(1,2)
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/cos_cot_integral.jl

@@ -34,7 +34,7 @@ printstyled("Integrating cos^2(x)/(1+cot(x)) from 0 to pi/2 and comparing the re
     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
+    @time @test simpsons_rule(cos_cot_fn, 39, 0, pi/2) ≈ 0.25
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(cos_cot_fn, 96, 0, pi/2) ≈ 0.25
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/cosine.jl

@@ -28,7 +28,7 @@ printstyled("Integrating cosine from 0 to pi/2 and comparing the result to the a
     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
+    @time @test simpsons_rule(x -> cos(x), 20, 0, pi/2) ≈ 1
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(x -> cos(x), 48, 0, pi/2) ≈ 1
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/expnx2datan.jl

@@ -34,7 +34,7 @@ printstyled("Integrating e^(-x^2)/(1+x^2) on the infinite domain [-inf, inf], or
     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
+    @time @test simpsons_rule(x -> exp(-x^2)/(x^2+1), 620, -100, 100) ≈ sol_9
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(x -> exp(-x^2)/(x^2+1), 1767, -100, 100) ≈ sol_9
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/gaussian.jl

@@ -34,7 +34,7 @@ printstyled("Integrate exp(-x^2) from minus infinity to positive infinity and co
     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)
+    @time @test simpsons_rule(x -> exp(-x^2), 268, -100, 100) ≈ sqrt(pi)
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(x -> exp(-x^2), 780, -100, 100) ≈ sqrt(pi)
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/logarithm.jl

@@ -28,7 +28,7 @@ printstyled("Integrating 1/x from 1 to e and comparing the result to the analyti
     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
+    @time @test simpsons_rule(x -> x^(-1), 34, 1, exp(1)) ≈ 1
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(x -> x^(-1), 84, 1, exp(1)) ≈ 1
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/logoverx.jl

@@ -27,7 +27,7 @@ printstyled("Integrating log(x)/x from 1 to e and comparing the result to the an
     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
+    @time @test simpsons_rule(x -> log(x)/x, 50, 1, exp(1)) ≈ 0.5
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(x -> log(x)/x, 114, 1, exp(1)) ≈ 0.5
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/modbessel0.jl

@@ -36,7 +36,7 @@ printstyled("Approximating the modified bessel function I_1(1) and comparing it
     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)
+    @time @test simpsons_rule(dmodified_bessel_1, 3, 0, pi/2) ≈ besseli(x,1)
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(dmodified_bessel_1, 9, 0, pi/2) ≈ besseli(x,1)
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/simppen.jl

@@ -30,7 +30,7 @@ printstyled("Integrating 1/sqrt(-19.6 sin(x)) from -pi to 0 and comparing the re
     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
+    @time @test abs(simpsons_rule(x -> (-19.6*sin(x))^(-0.5), 5e7, -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)

test/sinexpox.jl

@@ -40,7 +40,7 @@ printstyled("Integrating sin(x^2)e^(-x)/x from 0 to infinity, with the approxima
     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
+    @time @test simpsons_rule(sinexpox, 2101, 0, 100) ≈ sol_11
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(sinexpox, 5148, 0, 100) ≈ sol_11
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/sinxx.jl

@@ -36,7 +36,7 @@ printstyled("Integrating sin(x)/x from 0 to 100 and comparing it to the exact re
     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
+    @time @test simpsons_rule(sinxx, 339, 0, 100) ≈ sol_8
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(sinxx, 831, 0, 100) ≈ sol_8
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)

test/test_7.jl

@@ -35,7 +35,7 @@ printstyled("Integrating (x^3+1)/(x^4 (x+1)(x^2+1)) from 1 to e and comparing th
     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
+    @time @test simpsons_rule(partfrac, 100, 1, exp(1)) ≈ sol_7
     printstyled("Running: simpsons38_rule\n"; color = :magenta)
     @time @test simpsons38_rule(partfrac, 234, 1, exp(1)) ≈ sol_7
     printstyled("Running: trapezoidal_rule\n"; color = :magenta)