matlabsplinenumerical-integration

Double integration using experimental data


I need to perform a double integration using experimental data, but my integration limits are the same for each integral, in this case, the time. Mathematically I need to calculate:

E [ a0+ ∫0T a(t)dt ] = a + limTx → ∞ (1/T) ∫0T0t a dt dT

After some search, I reached to:

T = 0:0.1:600;
x = T;
A = rand(1,length(T)); % my data
pp_int = spline(T,A );
DoubleIntegration = integral(@(x)arrayfun(@(T )(integral(@(T ) ppval(pp_int,T ),0,  T  )),x),0,T(end)  );

The code took so long to run and give huge values. I think my problem is that Matlab might have trouble dealing with the splines, but I have no idea how to solve that.


Solution

  • First of all, you should not use the same letter for a bunch of things; your code is pretty hard to read because one has to figure out what T means at every instance.

    Second, pure math helps: after a change of variables and simple computation, the double integral becomes a single integral:

    0T0x a(t) dt dx = ∫0TtT a(t) dx dt = ∫0T (T-t)*a(t) dt

    I used non-random data on a smaller range for demonstration:

    T = 0:0.1:60;
    x = T;
    A = sin(T.^2); % my data
    pp_int = spline(T,A );
    tic
    DoubleIntegration = integral(@(x) arrayfun(@(T )(integral(@(T ) ppval(pp_int,T ),0,  T  )),x),0,T(end)  );
    toc
    tic
    SingleIntegration = integral(@(t) (T(end)-t).*ppval(pp_int,t), 0, T(end));
    toc
    disp(DoubleIntegration)
    disp(SingleIntegration)
    

    Output:

    Elapsed time is 27.751744 seconds.
    Elapsed time is 0.007223 seconds.
       51.3593
    
       51.3593
    

    Same result, 3800 times faster.