integral vs trapz vs sum
21 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Luqman Saleem
el 29 de En. de 2024
Comentada: John D'Errico
el 29 de En. de 2024
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 comentarios
Respuesta aceptada
John D'Errico
el 29 de En. de 2024
SIMPLE. READ THE HELP.
help 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]
2 comentarios
John D'Errico
el 29 de En. de 2024
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.
Más respuestas (0)
Ver también
Categorías
Más información sobre Numerical Integration and Differentiation en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!