Numerical Integration

integration.jl

using Base.Test
using Iterators: filter
using Distributions


function h_xs_ys(f, a, b, n_parts)
    h = (b - a) / n_parts
    xs = linspace(a, b, n_parts + 1)
    ys = [f(x) for x in xs]
    h, xs, ys
end


function torapezoidal_integration(f, a, b, n_parts=100)
    assert(b > a)
    h, xs, ys = h_xs_ys(f, a, b, n_parts)
    h * ((ys[1] + ys[end]) / 2 + sum(ys[2:end-1]))
end


function simpson_integration(f, a, b, n_parts=100)
    assert(b > a)
    h, xs, ys = h_xs_ys(f, a, b, 2*n_parts)
    s₁ = (ys[1]+ys[end])  #  y₀ + y₂ₙ
    s₂ = 4 * sum(ys[2:2:end-1])  # 4 * \sum y₁, y₃ .. y₂ₙ₋₁
    s₃ = 2 * sum(ys[3:2:end-2])  # 2 * \sum y₂, y₄ .. y₂ₙ₋₂
    h * (s₁ + s₂ + s₃) / 3
end


function simpson_integration2(f, a, b, n_parts=100)
    assert(b > a)
    h, xs, ys = h_xs_ys(f, a, b, 3*n_parts)
    s₁ = ys[1] + ys[end]  # y₀ + y₃ₙ
    s₂ = 3 * sum(ys[2:3:end-2])
    s₃ = 3 * sum(ys[3:3:end-1])
    s₄ = 2 * sum(ys[4:3:end-3])
    3h * (s₁ + s₂ + s₃ + s₄) / 8
end


is_in_positive_area(f, x, y) = f(x) > 0 && 0 <= y <= f(x)
is_in_negative_area(f, x, y) = f(x) < 0 && f(x) <= y <= 0


function monte_carlo_integration(f, xrange, yrange, n_samples=1e4)
    assert(xrange[2] > xrange[1])
    assert(yrange[2] > yrange[1])

    xs = rand(Uniform(xrange...), n_samples)
    ys = rand(Uniform(yrange...), n_samples)

    P = sum(is_in_positive_area(f, x, y) for (x, y) in zip(xs, ys))
    N = sum(is_in_negative_area(f, x, y) for (x, y) in zip(xs, ys))

    S = (xrange[2] - xrange[1]) * (yrange[2] - yrange[1])
    S * (P - N) / n_samples
end


function test_h_xs_ys()
    h, xs, ys = h_xs_ys(x -> x, 0, 10, 10)
    @test h == 1
    @test collect(xs) == collect(ys) == collect(0:1:10)
end


function test_integration()
    @test abs(torapezoidal_integration(sin, 0, 2π, 100) - 0) < 1e-8
    @test abs(torapezoidal_integration(cos, 0, π/2, 1e5) - 1) < 1e-8
    @test abs(simpson_integration(x->x, 0, 1, 100) - 1/2) < 1e-8
    @test abs(simpson_integration(sin, 0, 2π, 100) - 0) < 1e-8
    @test abs(simpson_integration(cos, 0, π/2, 1e5) - 1) < 1e-8
end


function test_is_in_area()
    f(x) = x

    @test is_in_positive_area(f, 6, 5)  # f(x) = 6 > 0 and 0 <= y = 5 < f(x)
    @test !is_in_positive_area(f, 3, 5)  # f(x) = 3 > 0 but 0 <= y = 5 ≰ f(x)
    @test !is_in_positive_area(f, -3, 5)  # f(x) = -3 ≯ 0


    @test is_in_negative_area(f, -5, -2)  # f(x) = -5 < 0 and f(x) <= y = -2 <= 0
    @test !is_in_negative_area(f, -3, -5)  # f(x) = -3 < 0 but f(x) ≰ y = -5 <= 0
    @test !is_in_negative_area(f, 3, -5) # f(x) = 3 ≮ 0
end


function run_homework()
    f₁(x) = 4 / (1+x^2)
    f₂(x) = sin(100x) + sin(57x) + 5
    range_f₁ = [0, 1]
    range_f₂ = [0, π/2]

    n_parts = 11  # divide the ranges into 10 chunks equally
    n_samples = 100000000

    println("Torapezoidal rule")
    println(
        "Integration of f₁ in range $range_f₁   : " *
        "$(torapezoidal_integration(f₁, range_f₁..., n_parts))"
    )
    println(
        "Integration of f₂ in range $range_f₂ : " *
        "$(torapezoidal_integration(f₂, range_f₂..., n_parts))"
    )

    println("")

    println("Simpson's 1/3 rule")
    println(
        "Integration of f₁ in range $range_f₁   : " *
        "$(simpson_integration(f₁, range_f₁..., n_parts))"
    )
    println(
        "Integration of f₂ in range $range_f₂ : " *
        "$(simpson_integration(f₂, range_f₂..., n_parts))"
    )

    println("")

    println("Simpson's 3/8 rule")
    println(
        "Integration of f₁ in range $range_f₁   : " *
        "$(simpson_integration2(f₁, range_f₁..., n_parts))"
    )
    println(
        "Integration of f₂ in range $range_f₂ : " *
        "$(simpson_integration2(f₂, range_f₂..., n_parts))"
    )


    println("")

    println("Monte Carlo")
    println(
        "Integration of f₁ in range $range_f₁   : " *
        "$(monte_carlo_integration(f₁, range_f₁, [0, 6], n_samples))"
    )
    println(
        "Integration of f₂ in range $range_f₂ : " *
        "$(monte_carlo_integration(f₂, range_f₂, [0, 8], n_samples))"
    )

end


test_h_xs_ys()
test_integration()
test_is_in_area()
run_homework()