Approximating Definite Integrals as Sums

I have a function that I want to integrate over and . I know the answer of this integration is . When I use function, I get the correct result but when I try to evaluate this integral using loops, my result is always . I increased points to but still get the same result. By using loops for integration. The sum approximation is . I am following a research paper which claims that they used loops with gaussian meshes and get as answer.

My question is why exactly I am not getting the correct result and what exactly are gaussian meshes? The function and the code runner.m that I am using are given below:

runner.m

%runner.m
clear; clc;

NBZ = 500; %number of x and y points. total points = NBZ^2

JNN = 0; %a parameter
a = 1;

% x and y limits... and x and y points. 
xmin = -2*pi/(3*a);
xmax = 4*pi/(3*a);
ymin = -2*pi/(sqrt(3)*a);
ymax = 2*pi/(sqrt(3)*a);
dx = (xmax-xmin)/(NBZ-1); %Delta x
dy = dx; %Delta y
xs = xmin:dx:xmax; %array of x points
ys = ymin:dy:ymax; %array of y points


dsum= 0;

for ny = 1:NBZ
    y = ys(ny);
    for nx = 1:NBZ
        x = xs(nx);
        out = F(x,y,JNN);
        dsum = dsum + out*dx*dy;
    end
end
answer = dsum
%it gives: answer = -0.6778

F(x,y)

function out = F(x,y,JNN)

JN = 4;
D = 1;
s = 1;

h11 = 0;
h12 = -(JN+1i*D)*s*(1+exp(1i*(-x-sqrt(3)*y)/2))...
    -JNN*s*(exp(1i*x)+exp(1i*(+x-sqrt(3)*y)/2));
h13 = -(JN-1i*D)*s*(1+exp(1i*(+x-sqrt(3)*y)/2))...
    -JNN*s*(exp(1i*x)+exp(1i*(-x-sqrt(3)*y)/2));
h23 = -(JN+1i*D)*s*(1+exp(1i*x))-2*JNN*s*exp(1i*x/2)*cos(sqrt(3)/2*y);
h = [h11 h12 h13;
    conj(h12) h11 h23;
    conj(h13) conj(h23) h11];
[evecs, evals] = eig(h);
v1 = evecs(:,1);
v2 = evecs(:,2);
v3 = evecs(:,3);
e1 = evals(1,1);
e2 = evals(2,2);
e3 = evals(3,3);

X = [                                                                                                                                                  0, - exp(- (x*1i)/2 - (3^(1/2)*y*1i)/2)*(1/2 - 2i) - JNN*(exp(x*1i)*1i + (exp((x*1i)/2 - (3^(1/2)*y*1i)/2)*1i)/2), - exp((x*1i)/2 - (3^(1/2)*y*1i)/2)*(1/2 + 2i) - JNN*(exp(x*1i)*1i - (exp(- (x*1i)/2 - (3^(1/2)*y*1i)/2)*1i)/2)
 - exp((conj(x)*1i)/2 + (3^(1/2)*conj(y)*1i)/2)*(1/2 + 2i) + JNN*((exp(- (conj(x)*1i)/2 + (3^(1/2)*conj(y)*1i)/2)*1i)/2 + exp(-conj(x)*1i)*1i),                                                                                                                   0,                                                     exp(x*1i)*(1 - 4i) - JNN*exp((x*1i)/2)*cos((3^(1/2)*y)/2)*1i
 - exp(- (conj(x)*1i)/2 + (3^(1/2)*conj(y)*1i)/2)*(1/2 - 2i) - JNN*((exp((conj(x)*1i)/2 + (3^(1/2)*conj(y)*1i)/2)*1i)/2 - exp(-conj(x)*1i)*1i),                                 exp(-conj(x)*1i)*(1 + 4i) + JNN*exp(-(conj(x)*1i)/2)*cos((3^(1/2)*conj(y))/2)*1i,                                                                                                                   0];
Y = [                                                                                                                                         0, - 3^(1/2)*exp(- (x*1i)/2 - (3^(1/2)*y*1i)/2)*(1/2 - 2i) + (3^(1/2)*JNN*exp((x*1i)/2 - (3^(1/2)*y*1i)/2)*1i)/2, 3^(1/2)*exp((x*1i)/2 - (3^(1/2)*y*1i)/2)*(1/2 + 2i) + (3^(1/2)*JNN*exp(- (x*1i)/2 - (3^(1/2)*y*1i)/2)*1i)/2
 - 3^(1/2)*exp((conj(x)*1i)/2 + (3^(1/2)*conj(y)*1i)/2)*(1/2 + 2i) - (3^(1/2)*JNN*exp(- (conj(x)*1i)/2 + (3^(1/2)*conj(y)*1i)/2)*1i)/2,                                                                                                                 0,                                                                  3^(1/2)*JNN*exp((x*1i)/2)*sin((3^(1/2)*y)/2)
   3^(1/2)*exp(- (conj(x)*1i)/2 + (3^(1/2)*conj(y)*1i)/2)*(1/2 - 2i) - (3^(1/2)*JNN*exp((conj(x)*1i)/2 + (3^(1/2)*conj(y)*1i)/2)*1i)/2,                                                       3^(1/2)*JNN*sin((3^(1/2)*conj(y))/2)*exp(-(conj(x)*1i)/2),                                                                                                               0];

o1a = ( ((v2'*X*v1)*(v1'*Y*v2)) - ((v2'*Y*v1)*(v1'*X*v2)))/((e2-e1)^2);
o1b = ( ((v3'*X*v1)*(v1'*Y*v3)) - ((v3'*Y*v1)*(v1'*X*v3)))/((e3-e1)^2);
o1 = 1i*(o1a+o1b);

out = real(o1/(2*pi));
end

Answer 1

[Matt J]

I don't know what gaussian meshes are, but the Riemann sum works with a just few fixes to the way the sampling is set up:

NBZ = 600; %number of x and y points. total points = NBZ^2

JNN = 0; %a parameter
a = 1;

% x and y limits... and x and y points. 
xmin = -2*pi/(3*a);
xmax = 4*pi/(3*a);
ymin = -2*pi/(sqrt(3)*a);
ymax = 2*pi/(sqrt(3)*a);

xs = linspace(xmin,xmax,NBZ+1); %array of x points
ys = linspace(ymin,ymax,NBZ+1); %array of y points
 xs(end)=[]; ys(end)=[];
dx=xs(2)-xs(1);
dy=ys(2)-ys(1);


dsum= 0;

for ny = 1:NBZ
    y = ys(ny);
    for nx = 1:NBZ
        x = xs(nx);
        out = F(x,y,JNN);
        dsum = dsum + out;
    end
end
answer = dsum*dx*dy

answer = -1.0000

[Luqman Saleem]

Thank you very much, it worked very well.