Approximating Definite Integrals as Sums

34 次查看（过去 30 天）
Luqman Saleem
Luqman Saleem 2020-11-2
评论： Luqman Saleem 2020-11-3
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

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

采纳的回答

Matt J
Matt J 2020-11-2
编辑：Matt J 2020-11-2
Ran in:
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by