 # [ANN] Trapz 2.0 - Fast integration routines for discrete domains

Hi everyone, this package is not exactly new, but i never announced it, and recently i had the time&need to make it even faster and flexible

# Trapz.jl

A simple Julia package to perform trapezoidal integration over common Julia arrays.

the package is now registered on Julia Registry, so it can be added as follows

``````
import Pkg
``````

## Example Usage:

``````using Trapz

vx=range(0,1,length=100)
vy=range(0,2,length=200)
vz=range(0,3,length=300)
M=[x^2+y^2+z^2 for x=vx,y=vy,z=vz]
I=trapz((vx,vy,vz),M)
print("result: ",I)
``````
``````result: 28.000303707970264
``````

# Benchmarks

``````using BenchmarkTools
@btime trapz(\$(vx,vy,vz),\$M);
``````
``````2.793 ms (4 allocations: 157.30 KiB)
``````
``````@btime trapz(\$(:,vy, vz),\$M);
``````
``````3.046 ms (3 allocations: 157.20 KiB)
``````
``````@btime trapz(\$(:,vy,:),\$M);
``````
``````4.022 ms (2 allocations: 234.45 KiB)
``````

# Comparison to Numpy

``````using PyCall
np=pyimport("numpy")

timejulia = @belapsed trapz(\$(vx,vy,vz),\$M)

print("Trapz.jl is ~ ",how_faster," times faster than numpy's trapz")
``````
``````Trapz.jl is ~ 7.2802292929009 times faster than numpy's trapz
``````
13 Likes

How fast is this relative to quadgk?

2 Likes

They have entirely different usecase, QuadGK needs a function, this operates on arrays. comparing their speed does not make much sense. in any case, quadgk can be pretty heavy on resources if you are not careful.

5 Likes

Nice! I like the new multidimensional integration syntax.

1 Like

If it is used appropriately, Gaussian quadrature can require exponentially fewer integrand evaluations than the trapezoidal rule as the required accuracy increases, assuming smooth functions. `quadgk` (and, in multiple dimensions, the HCubature package) are also adaptive (they put more function evaluations where the function is more oscillatory etc.).

As @francesco.alemanno says, the main utility of the trapezoidal rule is for cases where you cannot evaluate your function at arbitrary points in the domain, and only have a predetermined grid of points to work with.

6 Likes

Your PR to Trapz is to be thanked, i made some enhancements (like avoiding PermutedDimsArray) but it is almost like you designed it, thank you again @lstagner My experience is that for non-smooth/spiky functions the advantage is also enormous. Some functions seem more or less impossible to integrate accurately on a regular grid.