Numerically Integrating an Array

I have often integrated a function by using quad and dblquad. When using quad and dblquad, I can increase the tolerance parameter to the function to make sure that the answer is converging.

Instead of integrating a function, I need to now integrate an array. I know trapz can be used to do this. A simple example of this would be: t = linspace(0,10) y = cos(t) trapz(t,y)

However is there any way to adjust the tolerance or some other parameters, so that I can make sure that the answer is converging? If there is not any way to do this using the "trapz" function, is there another function you would recommend instead for numerically integrating the array?

Thank you, Kevin

P.S. The reason I need to integrate an array instead of a function is that I want to calculate Fourier coefficients for a curve I already have (which takes a lot of time for Matlab to generate). I don't want to have to re-generated the curve for each Fourier coefficient because that would be extremely computationally inefficient.

Answer 1

[John D'Errico]

Numerical analysis 101. The trapezoidal rule, when applied to a FIXED sequence of data points, will return a SINGLE estimate of the area contained. This estimate is based on the area underneath a piecewise linear function that in turn defines a set of trapezoids. Compute the area of the trapezoids, and form the sum. Thus is the origin of the name - trapezoidal rule. There is no control over a convergence tolerance here. None at all.

Yes, there are other methods. You might choose to read about other Newton-Cotes rules. Simpson's rule is one. If you look on the FEX, you will find Simpson's rule codes. I don't believe that you will find any higher order codes than this on the FEX. (Be careful, as Simpson's rule has in it an implicit requirement on the number of points used. This is true of any Newton-Cotes rule beyond trapezoidal rule. Some of those tools on the FEX forget to mention they have this flaw, or fail to deal with it creatively.)

Note that you have no control over the error tolerance for Simpson's rule either. It MAY often be more accurate than trapezoidal rule, however, I can show circumstances where too high of an order integration rule will degrade your results. Usually this problem results in the case of noise. If your data has noise in it, then the noise will become amplified. Of course, when working with floating point arithmetic, noise always enters your data, but it is really tiny noise, down in the least significant bits of your numbers.

Finally, nothing stops you from writing code for a higher order Newton-Cotes rule. Go too high in the order, and you may well be wasting your time. It depends on your data, the source of it, etc.

[Kevin Bachovchin]

Thanks for your response. Would increasing the number of data points in the curve be a way to test convergence using trapz, just as increasing the tolerance when using quad can check convergence?

I am just trying to avoid having to call my time-consuming function every time I calculate a different Fourier coefficient from the same curve. Any suggestions on the best way to do this?

Thank you,

Kevin

Answer 2

[Khaled Khairy]

My suggestion is to use Gaussian quadrature! (see gaussquad function in the matlab file exchange) The only restriction is that your function will have to be evaluated at the Gaussian base points. The integration is then just a weighted sum (no need for quad or dblquad). I hope this helps.

[John D'Errico]

You can't use Gaussian quadrature if the points are already given!

[Khaled Khairy]

True, however the way I understood the Post-Script of the original message, he can in fact generate the curve from scratch, it is just expensive to do so. Therefore I would generate it once using Gaussian quadrature basepoints and then do the weighted sum.