Integration of these kind of functions
4 次查看(过去 30 天)
显示 更早的评论
I want to integrate
is a 2 by 2 matrix (given below), is identity matrix of order 2, and Eand η and positive real numbers. I want to calculate this integration (for each component of the matrix) with .
Note that the integrand does not diverges as long as η is non-zero, however when η is significantly small, the integrand can be very large for some points in plane. As an example, I have plot the imaginary part of (1,2) componet of the integrand matrix for and below. (the other components also follow the same behavior):
It can be seen that when η is small the integrand is fairly smooth and can be integrated easily even with integrate2 function. But when η is small, the integrand has very sharp peaks due to which integral2 fails. To calculate the correct integration with small η by trapz, we need a lot of grid points which makes everything very slow.
Is there a way to calculate this integration with ?
I have analytically calculated the condition for which the integrand shows sharp peaks, that condition is:
where are components of matrix H. We could also find the exact points on which the condition is satisfied.
where x are the roots of equation and y is all values between -1 and +1. And
So, for each y there will be 4 roots, x, which will give eight pionts (four and four . From these eight pionts, we will take only the real ; imaginary must be neglected. Also for each y there will be two points (one and one .
Is there any way to tell MATLAB to take a lot of gird pionts near points and then calculate the integration? Or is there any way to calculate this type of integration?
Code:
clear; clc;
% parameters
E = 4.4;
eta = 0.1;5e-4;
% kx and ky limits and points
dkx = 0.006;
dky = dkx;
xmin = -2*pi/(3*sqrt(3));
xmax = 4*pi/(3*sqrt(3));
ymin = -2*pi/3;
ymax = 2*pi/3;
kxs = xmin:dkx:xmax;
kys = ymin:dky:ymax;
NBZx = length(kxs);
NBZy = length(kys);
% H matrix:
J = 1;
Dz = 0.5;
S = 1;
s0 = eye(2);
sx = [0, 1; 1, 0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
h0 = 3 * J * S;
hx = @(kx, ky) -J * S * (cos(ky / 2 - (3^(1/2) * kx) / 2) + cos(ky / 2 + (3^(1/2) * kx) / 2) + cos(ky));
hy = @(kx, ky) -J * S * (sin(ky / 2 - (3^(1/2) * kx) / 2) + sin(ky / 2 + (3^(1/2) * kx) / 2) - sin(ky));
hz = @(kx, ky) -2 * Dz * 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 = @(kx, ky) s0 * h0 + sx * hx(kx, ky) + sy * hy(kx, ky) + sz * hz(kx, ky);
%integrand:
G00 = @(kx, ky) inv(E*eye(2) - H(kx,ky) + 1i*eye(2)*eta); %integrand
2 个评论
David Goodmanson
2024-1-17
Hi Luqman,
since G00 appears to be a 2x2 matrix, what quantity are you plotting?
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!