Numerical integration of 2d array

[Solomon]

Hi everyone

I have a numerical integration problem. Here is the problem. I have a 6X6 matrix which is a function of 2 dimensional momenta (kx,ky). I want to numerically integrate the sum of the 3 largest eigenvalues, E4+E5+E6.

To calculate these eigenvalues, I loop over kx and ky in the range kx=(-4pi/3, 4pi/3) and ky=(-2pi/sqrt(3), 2pi/sqrt(3)) and calculate the sum of the eigenvalues, which is saved as a 2d array for each point in kx and ky, i.e fn(l,j)= E(4)+E(5)+E(6), where l runs over kx and j over ky.

Then, to integrate this function I did trapz(ky,trapz(kx,fn,2)). Unfortunately I am not getting the correct answer I should get. Are my missing something here. Can anyone help me resolve this problem.

Thanks,
Solomon.

[Roger Stafford]

"Solomon" <alao...@gmail.com> wrote in message <mpeets$qe3$1...@newscl01ah.mathworks.com>...

@@@@@@@@@@@@@@
I can think of at least four things to check on in your computations.

1. Is your double integration to be taken with respect to kx and ky? That is what your formula calls for.

2. Are you sure that the three maximum eigenvalues are always being taken? Remember, you cannot assume that 'eig' will always arrange them in any particular order.

3. Is the second dimension - that is the columns of your matrix - the one that has kx varying? That is what your formula calls for when you write:

trapz(kx,fn,2)

4. And finally you should verify that your method of computing the 36 elements of your 6 x 6 matrix as functions of fx and fy is correct.

You haven't stated how large your discrepancy is. Remember, 'trapz' only gives you an approximate answer whose accuracy depends on how closely-spaced your matrix values of kx and ky are. You could possibly achieve greater accuracy with the use of the quadrature function 'integral2'.

It would be wise to give further details of your computation. We can't really spot the errors with only the information you have given.

Roger Stafford

[Solomon]

Hi Roger,

Here are the answers to your questions.

1) The integration is to be taken with respect to kx and ky.

2). The eigenvalues are being sorted within the loop, i.e sort(....).

3). The second dimension has kx varying.

4). I don't know how to verify this.

The value of the integration gives 0.01348. I should get 0.2184. I think the error is too big with regard to what I am doing. I have tried it with 1000 points each for kx and ky and yet the answer did not change. Let me know if you have further questions.

"Roger Stafford" wrote in message <mpehi3$2ml$1...@newscl01ah.mathworks.com>...

[Solomon]

Hi Roger,

Here is the code

steps=1000;
kx = linspace(-2*pi/3,2*pi/3,steps);
ky =linspace(-pi/sqrt(3),pi/sqrt(3),steps);
ev1=zeros(length(kx));
for l= 1:length(kx)
for j= 1:length(ky)
k3= kx(l); k1= -(kx(l) + sqrt(3)*ky(j))/2; % wave vectors
k2= (-kx(l) + sqrt(3)*ky(j))/2;
ga=(1/12)*(exp(1i*k1)+exp(1i*k2)+exp(1i*k3));
gas=(1/12)*(exp(-1i*k1)+exp(-1i*k2)+exp(-1i*k3));
m1=[0, ga, gas; gas, 0, ga; ga, gas, 0]; m2= eye(3) +m1;
xx=[0,1; 1,0]; id2= eye(2); id3= eye(3);
zz=[1,0; 0,-1]; id6= kron(zz,id3);
H1= kron(m2, id2)-3*kron(m1, xx); H= id6*H1;
[V,eg]=eig(H);
egd=sort(real(diag(eg)));
ev1(l,j)= egd(4)+egd(5)+egd(6);
end
end
ev = 1-(1/(4*pi^2))*trapz(ky,trapz(kx,ev1,2))

"Roger Stafford" wrote in message <mpehi3$2ml$1...@newscl01ah.mathworks.com>...

[Roger Stafford]

"Solomon" <alao...@gmail.com> wrote in message <mpelhk$b8v$1...@newscl01ah.mathworks.com>...

>
> Here is the code
>
> steps=1000;
> kx = linspace(-2*pi/3,2*pi/3,steps);
> ky =linspace(-pi/sqrt(3),pi/sqrt(3),steps);
> ev1=zeros(length(kx));
> for l= 1:length(kx)
> for j= 1:length(ky)
> k3= kx(l); k1= -(kx(l) + sqrt(3)*ky(j))/2; % wave vectors
> k2= (-kx(l) + sqrt(3)*ky(j))/2;
> ga=(1/12)*(exp(1i*k1)+exp(1i*k2)+exp(1i*k3));
> gas=(1/12)*(exp(-1i*k1)+exp(-1i*k2)+exp(-1i*k3));
> m1=[0, ga, gas; gas, 0, ga; ga, gas, 0]; m2= eye(3) +m1;
> xx=[0,1; 1,0]; id2= eye(2); id3= eye(3);
> zz=[1,0; 0,-1]; id6= kron(zz,id3);
> H1= kron(m2, id2)-3*kron(m1, xx); H= id6*H1;
> [V,eg]=eig(H);
> egd=sort(real(diag(eg)));
> ev1(l,j)= egd(4)+egd(5)+egd(6);
> end
> end
> ev = 1-(1/(4*pi^2))*trapz(ky,trapz(kx,ev1,2))

The code you give with ev1(l,j) has kx varying along the first dimension of ev1, not its second, so you are incorrect when you say "The second dimension has kx varying." This means that trapz(kx,ev1,2) is incorrect by a scale factor of 2*sqrt(3). Moreover if the number of steps in kx and ky were different you would get an error message for the trapz call because the number of elements in kx would not be equal to the number of columns in ev1. However, with the number of steps the same, the outer trapz would have the opposite scale factor error which compensates for the inner scale factor error, giving you essentially the same answer as if you had written trapz(kx,trapz(ky,ev1,2)), so this error does not account for the numerical discrepancy between what you expect and what you get.

This apparently narrows the problem down to my fourth point, namely the way the elements of ev1 are computed, and I can't help you with this at present since I don't know the theory behind these computations. If you would give a reference to some accessible internet article that carries out this computation in non-matlab form, perhaps it would be possible for me or someone else to detect where your matlab code may go wrong.

[Solomon]

Hi Roger,

I computed the same exact integral of the eigenvalues in mathematica. It happens that the result of the trapz function is actually correct. It's a scientific article I wanted to confirm their results. I have sent the authors an e-mail with my mathematica file asking them to figure out how they obtained the result in their paper.

Thanks for the help anyway.




"Roger Stafford" wrote in message <mpg4l9$mko$1...@newscl01ah.mathworks.com>...