Numerical integral of a vector of discrete points

688 Views Asked by At

I was trying to obtain the integral of a vector of discrete points (with MATLAB), that do not follow a known function. For instance, I have a vector with discrete points defining a sine function and I want as output a vector of discrete points defining a (-) cosine function.

My idea was to implement the inverse thing of finite differences, solving a linear system in which I know the derivatives writing it in matrix form (A*x=b).

Example with 4th order central scheme, handling the edges with second order schemes, one-sided and central.

A = zeros(length(b), length(b));

A(1, 1:3) = [-3 4 -1] ./ (2*dx); % Fwd 2nd order
A(2, 1:3) = [-1 0 1] ./ (2*dx); % Central 2nd order

for ii = 3 : length(b) - 2
    A(ii, (1:5)+ii-3) = [1 -8 0 8 -1] ./ (12*dx); % Central 4th order
end
A(end-1, end-2:end) = [-1 0 1] ./ (2*dx); % Central 2nd order
A(end, end-2:end) = [1 -4 3] ./ (2*dx); % Bwd 2nd order

x = A \ b; % Solves the linear system

I'm getting singular or near-singular matrix. I'm not mathematician and I don't really know if this makes sense, but I think it should.

(Maybe BCs needed for this to work?)

Any ideas for fixing this or another approach are highly appreciated! Regards and thanks in advance.

2

There are 2 best solutions below

2
On

Your matrix is singular for the same reason that the derivative operation is non-invertible: it has a nontrivial nullspace, the space of constant functions. Primitive functions are only defined up to a constant term.

You can likely fix this by using the pseudoinverse (matlab pinv), which will give you the least-norm solution of your underdetermined equations.

Or use a quadrature rule, as pointed out in the comments.

0
On

In case someone is interested, I got the best results with MATLAB by first detrending the vector and then using cumtrapz which I guess it is equivalent to the usage of cumsum described by Ian above.

cumtrapz(detrend(b,'constant'))