Numerical integration only with evaluated point

New to Julia
1

Hi,
I’d like to integrate a function numerically.
I found some packages, e.g., QuadGK.jl, it seems only supports numerical integration with a given function.

In my case, suppose we cannot access the function at arbitrary points.
That is, we can access only some given data points.

2

Check this post out about Romberg method or NumericalIntegration.jl.

1 Like
3

You just need to define quadrature weights.
For example, in FastGaussQuadrature.jl

using FastGaussQuadrature
p, w = gausslegendre(10)

Let’s say you have 10 values at these points, then do an integration

vals = randn(10)
integral = sum(vals .* w)
4

Suppose a change of interval is required when using the Gauss-Legendre quadrature and that it will not handle arbitrary input points, as it might be the case in OP’s question?

5

You can always regularize your Gaussian weights according to the geometric transformation on points, right?
I have not used the library you mentioned. Maybe it provides a straightforward solution though.

6

Romerg, mentioned by @rafael.guerra above, is surprisingly accurate for integrating smooth functions but it requires the values being integrated be evenly sampled over the independent variable.

You can always fall back to trapezoidal integration if the points are not evenly sampled, but it’s low accuracy.

I"m not sure how the quadrature weighting method would compare in performance and accuracy. Need to try it out.

4 Likes
7

Thanks for the tip, @vavrines! One resource I found was slide 32 of these slides from Janet Peterson on numerical quadrature – they explain how to compute the transformed points and weights.

I have implemented these transformations to the quadrature points and weights from FastGaussQuadrature.jl.

using FastGaussQuadrature, LinearAlgebra

"""Transform the Gauss quadrature points x and weights w from [-1,1] to [a,b]. 

   Ref: https://people.sc.fsu.edu/~jpeterson/numerical_quadrature
"""
function transform_gauss_xw(x, w, a, b)
    x1 = a .+ (b-a)/2 * (x .+ 1)
    w1 = (b-a)/2 * w
    return (x=x1, w=w1)
end

"""Calculate the Gauss quadrature integral of f(x) from `a` to `b`, using `order` quadrature points.

    Example:
    ```
    f(x) = sin(x) * x^2
    LB = 0; UB = 10; n = 10;
    I = GL_integral(f, LB, UB, n)
    ```

    Ref: https://juliaapproximation.github.io/FastGaussQuadrature.jl/stable/gaussquadrature/
"""
function GL_integral(f, a, b, order)
    x, w = gausslegendre(order)
    T = transform_gauss_xw(x, w, a, b)
    x1 = T.x; w1 = T.w;
    I = dot(w1, f.(x1))
    return I
end
1 Like