Numerical integration over non-uniform grid in matlab. Is there any function?

Question

I've got function values in a vector f and also the vector containing values of the argument x. I need to find the define integral value of f. But the argument vector x is not uniform. Is there any function in Matlab that deals with integration over non-uniform grids?

Answer 1

Taken from help :

Z = trapz(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, or X must be a column vector and Y an array whose first non-singleton dimension is length(X). trapz operates along this dimension.

As you can see x does not have to be uniform.

For instance:

x = sort(rand(100,1)); %# Create random values of x in [0,1]
y = x;
trapz( x, y) 

Returns:

ans =

    0.4990

Another example:

x = sort(rand(100,1)); %# Create random values of x in [0,1]
y = x.^2;
trapz( x, y) 

returns:

ans =

    0.3030

Answer 2

Depending on your function (and how x is distributed), you might get more accuracy by doing a spline interpolation through your data first:

pp  = spline(x,y);
quadgk(@(t) ppval(pp,t), [range]) 

That's the quick-n-dirty way. Ther is a faster and more direct approach, but that is fugly and much less transparent:

result = sum(sum(...
    bsxfun(@times, pp.coefs, 1./(4:-1:1)) .*...  % coefficients of primitive
    bsxfun(@power, diff(pp.breaks).', 4:-1:1)... % all 4 powers of shifted x-values
    ));

As an example why all this could be useful, I borrow the example from here. The exact answer should be

>> pi/2/sqrt(2)*(17-40^(3/4))
ans =
     1.215778726893561e+00

Defining

>> x = [0 sort(3*rand(1,5)) 3];
>> y = (x.^3.*(3-x)).^(1/4)./(5-x);

we find

>> trapz(x,y)
ans =
    1.142392438652055e+00

>> pp  = spline(x,y);
>> tic; quadgk(@(t) ppval(pp,t), 0, 3), toc
ans =
    1.213866446458034e+00
Elapsed time is 0.017472 seconds.

>> tic; result = sum(sum(...
    bsxfun(@times, pp.coefs, 1./(4:-1:1)) .*...  % coefficients of primitive
    bsxfun(@power, diff(pp.breaks).', 4:-1:1)... % all 4 powers of shifted x-values
    )), toc
result =
    1.213866467945575e+00
Elapsed time is 0.002887 seconds.

So trapz underestimates the value by more than 0.07. With the latter two methods, the error is an order of magnitude less. Also, the less-readable version of the spline approach is an order of magnitude faster.

So, armed with this knowledge: choose wisely :)