integral vs trapz vs sum

Question

Luqman Saleem 2024-1-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2075746-integral-vs-trapz-vs-sum

评论： John D'Errico 2024-1-29

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] 

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2024-1-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2075746-integral-vs-trapz-vs-sum#answer_1399086

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Luqman Saleem 2024-1-29

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

thank you very much for the reply.

John D'Errico 2024-1-29

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.

请先登录，再进行评论。