File Exchange

image thumbnail

Fast Trapezoidal Integration

version 1.0 (3.09 KB) by

Similar to TRAPZ but (up to) 16x-times faster. Equally spaced data NOT required.

7 Downloads

Updated

No License

F=trapzf(X,Y) computes the integral of Y with respect to X using the trapezoidal method. X and Y must be vectors of the same length.
For trapzf to work properly:
X must be a STRICTLY monotonically increasing vector.

Example:
>> x=[1:1:1000];
>> y=log(sqrt(x+1.001)-1.001);
>> area=trapzf(x,y);

Warning: not much in the way of error checking, since this slows things down, so pay attention to the argument passed to the function!!!

Installation (Windows users):
------------------------------

Simply copy the trapzf.dll file in a directory recognizable by MATLAB. You won't need trapzf.c to perform computations, but it is useful if you want to read the user instructions or to customize your code.
If you modify trapzf.c then you must run the following command

mex trapzf.c

in order to obtain the new dll file.

Installation (other OS users):
-------------------------------

A precompiled c-mex file is not provided.
You have to compile trapzf.c with your favourite C compiler (write "mex -setup" at the MATLAB prompt), then write

mex trapzf.c

in order to obtain the corresponding mex-file (see the table below)

        OS ---> Extension
      ------------------------------
        sol2, SunOS 5.x ---> .mexsol
        hpux ---> .mexhpux
        hp700 ---> .mexhp7
        ibm_rs ---> .mexrs6
        sgi ---> .mexsg
        alpha ---> .mexaxp
        glnx86 ---> .mexglx

PLEASE RATE MY SUBMISSION!

Comments and Ratings (5)

Lance Optican

trapzf does not allow for unevenly sampled data, even though the comment claims it does. Matlab's trapz does allow for uneven data. To correct this, the following change must be made to the C code of trapzf, here called trapzfu. Note that the speedup on a 64-bit WinXP computer of trapfzu over trapz was about 10-fold. Unrolling the loop in trapzfu2.c provides another 20% improvement in speed.

--- trapzfu.c ------
/* F=trapzfu(X,Y) computes the integral of Y with respect to X using
    the trapezoidal method. X and Y must be vectors of the same length.
    For trapzf to work properly:
    X must be a STRICTLY monotonically increasing vector.
    X does not need to be evenly spaced.
 
    C. Quaia, National Eye Institute, NIH, DHHS
 */

#include "mex.h"

void trapzfu( double *area, double *xv, double *yv, int N) {
        register int i;
        
        *area = 0;
        for (i = 0; i < (N-1); i++) {
                *area = *area + (xv[i+1] - xv[i]) * (yv[i] + yv[i+1]);
        }
        *area = *area / 2;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        int N, N1, N2;
        double *area;
        double *xv, *yv ;
        
        N1 = mxGetNumberOfElements(prhs[0]);
        N2 = mxGetNumberOfElements(prhs[1]);
        if (N1 != N2)
                mexErrMsgTxt("Input x and y must have the same length");
        else
                N = N1;
        if (N < 2)
                mexErrMsgTxt("The minimum length is 2");
        
        xv = mxGetPr(prhs[0]);
        yv = mxGetPr(prhs[1]);
        plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
        area = mxGetPr(plhs[0]);
        
        trapzfu(area, xv, yv, N);
}

------------ trapzfu2.c --------
/* F=trapzfu2(X,Y) computes the integral of Y with respect to X using
    the trapezoidal method. X and Y must be vectors of the same length.
    For trapzf to work properly:
    X must be a STRICTLY monotonically increasing vector.
    X does not need to be evenly spaced.

    L. Optican, NEI/NIH/DHHS

 */

#include "mex.h"

void trapzfu2( double *area, double *xv, double *yv, int N) {
        register int i, j,r;
        
        *area = 0;
        i = 0;
        r = (N - 1) % 4;
        for(; i < r; i++) {
                // *area += (xv[i+1] - xv[i]) * (yv[i] + yv[i+1]);
                *area += (xv[i+1] - xv[i]) * (yv[i+1] + yv[i]);
        }
        
        for(; i <= (N-4); i += 4) {
                *area += (xv[i+1] - xv[i]) * (yv[i+1] + yv[i]);
                *area += (xv[i+2] - xv[i+1]) * (yv[i+2] + yv[i+1]);
                *area += (xv[i+3] - xv[i+2]) * (yv[i+3] + yv[i+2]);
                *area += (xv[i+4] - xv[i+3]) * (yv[i+4] + yv[i+3]);
        }
        *area = *area / 2;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        int N, N1, N2;
        double *area;
        double *xv, *yv ;
        
        N1 = mxGetNumberOfElements(prhs[0]);
        N2 = mxGetNumberOfElements(prhs[1]);
        if (N1 != N2)
                mexErrMsgTxt("Input x and y must have the same length");
        else
                N = N1;
        if (N < 2)
                mexErrMsgTxt("The minimum length is 2");
        
        xv = mxGetPr(prhs[0]);
        yv = mxGetPr(prhs[1]);
        plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
        area = mxGetPr(plhs[0]);
        
        trapzfu2(area, xv, yv, N);
}

lll hhh

not good

mm prabhu2003

thanks

sind baad

Jay Leno

Updates

1) Bug corrected: now trapzf(x,y) returns zero if length(x)=length(y)=1.
2) A bit of effort on dimension check of x and y has been added (just a little...to not decrease the algorithm speed)

no significative changes

No changes

no significative changes

Requirements added

this version has been created with the LCC compiler (instead of Visual Studio 6); further, the installation steps description have been improved.

no significant changes

Now the inputs can be either row or column arrays.

no modifications: I just added my request to give a rate to this file.

MATLAB Release
MATLAB 6.5 (R13)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video