integral vs trapz vs sum

The result obtained from trapz() differs in sign compared to the approximation from integral() and the summation approach for the following integration (note that integration is taken over E). From the following code, I get:

1.7136 from integral(),

-1.7099 from trapz(),

1.7136 from sum().

Why trapz() sign is different?

clear; clc;

%----------------------------- some constants -----------------------------
kx = -2*pi/(3*sqrt(3));
ky = 1*pi/3;
T = 10;
J = 1;
D = 0.8;
S = 1;
eta = 0.01;

beta = 1/(0.0863*T);

s0 = eye(2);
sx = [0, 1; 1, 0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];

h0 = 3 * J * S;
hx = -J * S * (cos(ky / 2 - (3^(1/2) * kx) / 2) + cos(ky / 2 + (3^(1/2) * kx) / 2) + cos(ky));
hy = -J * S * (sin(ky / 2 - (3^(1/2) * kx) / 2) + sin(ky / 2 + (3^(1/2) * kx) / 2) - sin(ky));
hz = -2 * D * S * (sin(3^(1/2) * kx) + sin((3 * ky) / 2 - (3^(1/2) * kx) / 2) - sin((3 * ky) / 2 + (3^(1/2) * kx) / 2));

H = s0*h0 + sx*hx + sy*hy + sz*hz;

dkx_hx = -J*S*((3^(1/2)*sin(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*sin(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hy = J*S*((3^(1/2)*cos(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*cos(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hz = 2*D*S*((3^(1/2)*cos((3*ky)/2 - (3^(1/2)*kx)/2))/2 - 3^(1/2)*cos(3^(1/2)*kx) + (3^(1/2)*cos((3*ky)/2 + (3^(1/2)*kx)/2))/2);

dky_hx = J*S*(sin(ky/2 - (3^(1/2)*kx)/2)/2 + sin(ky/2 + (3^(1/2)*kx)/2)/2 + sin(ky));
dky_hy = -J*S*(cos(ky/2 - (3^(1/2)*kx)/2)/2 + cos(ky/2 + (3^(1/2)*kx)/2)/2 - cos(ky));
dky_hz = -2*D*S*((3*cos((3*ky)/2 - (3^(1/2)*kx)/2))/2 - (3*cos((3*ky)/2 + (3^(1/2)*kx)/2))/2);

X = sx*dkx_hx + sy*dkx_hy + sz*dkx_hz;
Y = sx*dky_hx + sy*dky_hy + sz*dky_hz;

G0R = @(E) inv(E*s0 - H + 1i*eta*s0);
fE = @(E) 1/( exp(beta*E) - 1 );


%--------------------------------------------------------------------------
%function:
fun = @(E) real(trace(1/2*fE(E)*E^2*(G0R(E)*X*G0R(E)*Y*G0R(E) - G0R(E)*Y*G0R(E)*X*G0R(E)) + (1/2*fE(E)*E^2*(G0R(E)*X*G0R(E)*Y*G0R(E) - G0R(E)*Y*G0R(E)*X*G0R(E)))'));
funn = @(E) arrayfun( @(E)fun(E), E ); %arrayfun


%limits:
Emin = -10;
Emax = 30;
dE = 1e-4;
Es = Emin:dE:Emax;


%integration using integral()
I_int = integral(funn,Emin,Emax); 


%integration using trapz()
funn_values = funn(Es);
[~,indx] = find(isnan(funn_values)); % replacing NaN with mean value
for idx = indx
    funn_values(idx) = ( funn_values(idx-1)+funn_values(idx+1) )/2;
end
I_trapz = trapz(funn_values,Es);


%integration using sum()
I_sum  = sum(funn_values(:))*dE;


I = [I_int, I_trapz, I_sum] 


% output: [1.7136   -1.7099    1.7136] 

Answer 1

[John D'Errico]

SIMPLE. READ THE HELP.

help trapz

 TRAPZ  Trapezoidal numerical integration.
    Z = TRAPZ(Y) computes an approximation of the integral of Y via
    the trapezoidal method (with unit spacing).  To compute the integral
    for spacing different from one, multiply Z by the spacing increment.
 
    For vectors, TRAPZ(Y) is the integral of Y. For matrices, TRAPZ(Y)
    is a row vector with the integral over each column. For N-D
    arrays, TRAPZ(Y) works across the first non-singleton dimension.
 
    Z = TRAPZ(X,Y) computes the integral of Y with respect to X using the
    trapezoidal method. X can be a scalar or a vector with the same length
    as the first non-singleton dimension in Y. TRAPZ operates along this
    dimension. If X is scalar, then TRAPZ(X,Y) is equivalent to X*TRAPZ(Y).
 
    Z = TRAPZ(X,Y,DIM) or TRAPZ(Y,DIM) integrates across dimension DIM
    of Y. The length of X must be the same as size(Y,DIM)).
 
    Example:
        Y = [0 1 2; 3 4 5]
        trapz(Y,1)
        trapz(Y,2)
 
    Class support for inputs X, Y:
       float: double, single
 
    See also SUM, CUMSUM, CUMTRAPZ, INTEGRAL.

    Documentation for trapz
       doc trapz

    Other uses of trapz

       codistributed/trapz    gpuArray/trapz

So the call to trapz should have been

I_trapz = trapz(Es,funn_values);

Changing that,

kx = -2*pi/(3*sqrt(3));
ky = 1*pi/3;
T = 10;
J = 1;
D = 0.8;
S = 1;
eta = 0.01;
Emin = -10;
Emax = 30;
dE = 1e-4;
Es = Emin:dE:Emax;
NE = length(Es);
beta = 1/(0.0863*T);
s0 = eye(2);
sx = [0, 1; 1, 0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
h0 = 3 * J * S;
hx = -J * S * (cos(ky / 2 - (3^(1/2) * kx) / 2) + cos(ky / 2 + (3^(1/2) * kx) / 2) + cos(ky));
hy = -J * S * (sin(ky / 2 - (3^(1/2) * kx) / 2) + sin(ky / 2 + (3^(1/2) * kx) / 2) - sin(ky));
hz = -2 * D * S * (sin(3^(1/2) * kx) + sin((3 * ky) / 2 - (3^(1/2) * kx) / 2) - sin((3 * ky) / 2 + (3^(1/2) * kx) / 2));
H = s0*h0 + sx*hx + sy*hy + sz*hz;
dkx_hx = -J*S*((3^(1/2)*sin(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*sin(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hy = J*S*((3^(1/2)*cos(ky/2 - (3^(1/2)*kx)/2))/2 - (3^(1/2)*cos(ky/2 + (3^(1/2)*kx)/2))/2);
dkx_hz = 2*D*S*((3^(1/2)*cos((3*ky)/2 - (3^(1/2)*kx)/2))/2 - 3^(1/2)*cos(3^(1/2)*kx) + (3^(1/2)*cos((3*ky)/2 + (3^(1/2)*kx)/2))/2);
dky_hx = J*S*(sin(ky/2 - (3^(1/2)*kx)/2)/2 + sin(ky/2 + (3^(1/2)*kx)/2)/2 + sin(ky));
dky_hy = -J*S*(cos(ky/2 - (3^(1/2)*kx)/2)/2 + cos(ky/2 + (3^(1/2)*kx)/2)/2 - cos(ky));
dky_hz = -2*D*S*((3*cos((3*ky)/2 - (3^(1/2)*kx)/2))/2 - (3*cos((3*ky)/2 + (3^(1/2)*kx)/2))/2);
X = sx*dkx_hx + sy*dkx_hy + sz*dkx_hz;
Y = sx*dky_hx + sy*dky_hy + sz*dky_hz;
G0R = @(E) inv(E*s0 - H + 1i*eta*s0);
fE = @(E) 1/( exp(beta*E) - 1 );
%--------------------------------------------------------------------------
%function:
fun = @(E) real(trace(1/2*fE(E)*E^2*(G0R(E)*X*G0R(E)*Y*G0R(E) - G0R(E)*Y*G0R(E)*X*G0R(E)) + (1/2*fE(E)*E^2*(G0R(E)*X*G0R(E)*Y*G0R(E) - G0R(E)*Y*G0R(E)*X*G0R(E)))'));
funn = @(E) arrayfun( @(E)fun(E), E ); %arrayfun
%integration using integral()
I_int = integral(funn,Emin,Emax); 
%integration using trapz()
funn_values = funn(Es);
[~,indx] = find(isnan(funn_values)); % replacing NaN with mean value
for idx = indx
    funn_values(idx) = ( funn_values(idx-1)+funn_values(idx+1) )/2;
end
I_trapz = trapz(Es,funn_values);
%integration using sum()
I_sum  = sum(funn_values(:))*dE;
I = [I_int, I_trapz, I_sum] 

I = 1×3

1.7136 1.7136 1.7136

[Luqman Saleem]

Lol. How did I miss that... I feel like an idiot now.

thank you very much for the reply.

[John D'Errico]

Lol. Honestly, I've tripped on it a few times in the past too. At least until you trip enough times that you no longer stub your big toe there.

Logically, if I were to write trapz, then if y is the first argument if called as trapz(y), then if you call it with TWO arguments, it SHOULD be trapz(y,x). SERIOUSLY, trapz is just begging for new users to trip over that. Sadly, I did not write trapz.