How to integrate matrices in MATLAB ?

I want to do the following integration in MATLAB, given by this equation

The data stored in ue is from a numerical simulation (i.e Discrete Data). T and Ly are given and are 5 and 80, respectively. The total simulation time is 500 s ( dy = 0.1 and dt = 0.01).

ue - is a matrix with 3x3 dimensions where u (t,y,x), t being the time of the points stored in the matrix, y is the points in the y-drection, x is the points in the x-direction.

My code is :

% Inner Integration

for t = 1:totaltime
    for i=1:x
        u_integrated(t,:,i) = (trapz(0:dy:Ly,ue(t,:,i))/Ly);
    end
end

% Outer Integration

for t = 1:totaltime
   for j = 1:y
      for i = 1:x
        u_outint(t,j,i) = (trapz(0:dt:T,u_integrated(t,j,i)))/(T);
      end
   end
end

I have tried to use the trapz function in MATLAB, but I am not sure whether my implementation is correct or not.

[Jan]

What is "u2"? Do you mean "u_integrated"?

[Abdulhalim Saeed]

Hi Jan, I have made a mistake u2 should be replaced by u_integrated, u_integrated is the result from the first integration ( the inner integration in the equation that I have written).

[Torsten]

You integrate a three-dimensional array with respect to two of its dimensions. So the result should be a one-dimensional array only depending on the x-coordinate. Thus your solution cannot be correct.

Answer 1

[Torsten]

% Generate a data matrix
fun = @(t,y,x) x.^2+y+sin(t);
X = 0:0.1:1;
Y = 1:0.05:2;
T = 1:1:10;
for it = 1:numel(T)
    for iy = 1:numel(Y)
        for ix = 1:numel(X)
            Z(it,iy,ix) = fun(T(it),Y(iy),X(ix));
        end
    end
end
Lend = Y(end);
Tend = T(end);

% Integrate with respect to y
for it = 1:numel(T)
    for ix = 1:numel(X)
        integral_y(it,ix) = trapz(Y,Z(it,:,ix));
    end
end
%Integrate the result with respect to t
for ix = 1:numel(X)
    integral_ty(ix) = trapz(T,integral_y(:,ix));
end
integral_ty = integral_ty/(Lend*Tend);

% Integrate analytically
for ix = 1:numel(X)
  integral_ty_ana(ix) = integral2(@(t,y)fun(t,y,X(ix)),T(1),Tend,Y(1),Lend);
end
integral_ty_ana = integral_ty_ana/(Lend*Tend);

% Compare results
plot(X,[integral_ty;integral_ty_ana])