Double integration using experimental data

Question

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 [ a₀+ ∫₀^T a(t)dt ] = a + lim_{Tx → ∞} (1/T) ∫₀^T ∫₀^t 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.

Answer 1

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:

∫₀^T ∫₀^x a(t) dt dx = ∫₀^T ∫_t^T a(t) dx dt = ∫₀^T (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.