Approximating Definite Integrals as Sums

22 views (last 30 days)
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

Accepted Answer

Matt J
Matt J on 2 Nov 2020
Edited: Matt J on 2 Nov 2020
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
  1 Comment
Luqman Saleem
Luqman Saleem on 3 Nov 2020
Thank you very much, it worked very well.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by