# Approximating Definite Integrals as Sums

22 views (last 30 days)
Luqman Saleem on 2 Nov 2020
Commented: Luqman Saleem on 3 Nov 2020
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
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

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
Luqman Saleem on 3 Nov 2020
Thank you very much, it worked very well.

### Categories

Find more on Robotics in Help Center and File Exchange

R2020b

### Community Treasure Hunt

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

Start Hunting!

Translated by