7

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?

4

3 回答 3

6

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
于 2012-11-15T11:27:24.623 回答
3

根据您的功能（以及x分布方式），您可以spline通过首先对数据进行插值来获得更高的准确性：

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

那是快速n肮脏的方式。还有一种更快、更直接的方法，但这种方法很丑陋，而且不太透明：

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
    ));

作为一个为什么所有这些都可能有用的例子，我从这里借用了这个例子。确切的答案应该是

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

定义

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

我们发现

>> 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.

所以trapz低估了价值超过0.07. 使用后两种方法，误差要小一个数量级。此外，该spline方法的可读性较差的版本要快一个数量级。

所以，有了这些知识：明智地选择:)

于 2012-11-15T11:59:32.647 回答
0

您可以对每对分段进行高斯求积x并将它们相加以获得完整的积分。

于 2012-11-15T13:28:32.203 回答