Approximating Definite Integrals as Sums

Question

Luqman Saleem 2020-11-2

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/633819-approximating-definite-integrals-as-sums

评论： Luqman Saleem 2020-11-3

采纳的回答： Matt J

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2020-11-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/633819-approximating-definite-integrals-as-sums#answer_531714

编辑：Matt J 2020-11-2

在 MATLAB Online 中打开

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