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; endendanswer = 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
Best Answer