CODEWITHSUNDEEP
matlabnumericalnumerical-integration
rebatoma
rebatoma

Reputation: 742

integral or trapz, which one is the more appropriate in MATLAB?

I'm computing multiple integrals using MATLAB.

I'm using the integral function to compute the integral but I was wondering is it faster to use trapz instead of using integral?

I know that trapz introduces a bit of error in the computation, but despite that, with is the best function to compute integrals in MATLAB?

Upvotes: 1

Views: 7412

Answers (1)

TroyHaskin
TroyHaskin

Reputation: 8401

Short and sweet:

  • Use trapz for discrete data or for selected functional data if you don't care about (potentially extremely) low accuracy of the integral value
  • Use integral for integrands that have a functional form, adjusting tolerances as needed for speed.

As mentioned by the MATLAB documentation, trapz is intended "to perform numerical integrations on discrete data sets" and leverages the trapezoidal rule for the integrations. The error between the true integral and the trapz approximation is almost entirely dependent on the input x vector (sometimes called the abscissa in integration parlance) with no automatic adaptability. The good part is that if the underlying function is "nice" (i.e., continuous, smooth, no sharp peaks or excessive oscillations, etc.), trapz will likely be the fastest function to approximate the integral since it

  1. Doesn't have to call a function for values (they're input)

  2. Doesn't automatically adapt (which takes time and can be complex to implement).

However, for general integrals, trapz may also be the most inaccurate and may require a denser x vector to calculate a low-error value.

For discrete data, this is a short-coming that must be lived with, but if the integrand has a functional form, integral and its family is highly recommended.

The black-box numeric integrators in MATLAB have evolved over the years, and MathWorks co-founder Clever Moler has a nice blog post going over some of the evolutions. The post discusses the quad, quadl, and quadgk functions and how quadgk became the core for integral and its ilk. The basic breakdown of the three functions is

  • quad uses a three-point and five-point Simpson's Rule
  • quadl uses a four-seven-thirteen point1 Lobatto-Kronrod2 rule
  • quadgk a uses seven-fifteen point Gauss-Kronrod2 rule

to acquire both an approximation of the integral and an error approximation for adaptive quadrature. The summary of the history lesson and test problems is that quadgk was written with vectorization incorporated3, uses a higher-order rule which excludes end-points, and gives extremely accurate answers faster than its competitors. As a result, quadgk is the core of the new and highly-recommended integral family.

1 Adaptive quadrature usually lists the number of points used to form its approximation of the value and the error. Typically, there are two numbers that indicate the number of points to form the low-order and high-order approximations. quadl is interesting in that it uses a four-point Gauss-Lobatto rule and seven-point and thirteen-point Kronrod extensions for its error handling.

2 Gaussian Quadrature, which is an integration technique that chooses it abscissa to exactly integrate a family of polynomials over a given interval instead of prescribing them as in Newton-Cotes, has a lot of names associated with it to indicate a lot of "stuff" that's going on without being explicit about it (which can be very annoying to newcomers). "Gauss" refers to the aforementioned method of choosing abscissa and associated weights for the integration. "Lobatto" indicates an extension to Gauss-Legendre integration methods that incorporates end-points (others may not like my link between these two, but I find the parallels pleasing). "Kronrod" indicates an extension to any particular Gauss rule that creates a high-order rule using a given set of abscissa and adding to it; this creates a "nesting" (the low-order points are part of the high-order point set) that results in fewer function evaluations overall.

3 Since vectorization is written into integral, integrands or limits that are vector-valed must use the 'ArrayValued' flag to tell the program to make functional evaluations differently so as not to create a size-mismatch error. It might be possible to program around this to a certain extent, but the MathWorks decided not to.

Upvotes: 6

Related Questions