# Matlab double sum over vectors

2 vues (au cours des 30 derniers jours)
CC SS le 17 Sep 2023
I want to calculate , in which k and k' are vectors with x and y components. Both k and k' lie in the 1st Brillouin zone: -pi<kx<pi, -pi<ky<pi, -pi<kx'<pi, -pi<ky'<pi. How to write Matlab code to perform the sum over k and k'? My code doesn't work:
sum = 0;
kvec = linspace(-pi,pi,N);
for j1=1:1:length(kvec)
kx = kvec(j1);
for j2=1:1:length(kvec)
ky = kvec(j2);
for j3=1:1:length(kvec)
kxprime = kvec(j3);
for j4=1:1:length(kvec)
kyprime = kvec(j4);
sum = sum + f(kx,ky,kxprime,kyprime)
end
end
end
end
##### 4 commentairesAfficher 2 commentaires plus anciensMasquer 2 commentaires plus anciens
Dyuman Joshi le 17 Sep 2023
Can you share the definition of f?
CC SS le 17 Sep 2023
f = 2* (cos(kx) + cos(ky)) - 0.5* (cos(kxprime) + cos(kyprime)).

Connectez-vous pour commenter.

### Réponse acceptée

Walter Roberson le 17 Sep 2023
Different versions, and their timings.
Except... in different runs, the timings especially for the first version varied by more than a factor of 10, so the values shown here for a single run might not be representative of real use.
format long g
tic
N = 10;
kvec = linspace(-pi,pi,N);
ck = cos(kvec);
ca = reshape(ck, [], 1);
cb = reshape(ck, 1, []);
cc = reshape(ck, 1, 1, []);
cd = reshape(ck, 1, 1, 1, []);
f = (ca + cb)*2 - (cc + cd)/2;
out1 = sum(f,'all')
out1 =
-2999.99999999997
toc
Elapsed time is 0.022545 seconds.
tic
f = @(a,b,c,d) 2*(cos(a) + cos(b)) - 0.5*(cos(c) + cos(d));
%Random value for N for example
N = 10;
kvec = linspace(-pi,pi,N);
[kx,ky,kxprime,kyprime] = ndgrid(kvec);
arr = f(kx,ky,kxprime,kyprime);
out2 = sum(arr,'all')
out2 =
-2999.99999999997
toc
Elapsed time is 0.007611 seconds.
tic
out3 = 0;
N = 10;
kvec = linspace(-pi,pi,N);
f = @(k,kprime) 2* (cos(k(1)) + cos(k(2))) - 0.5* (cos(kprime(1)) + cos(kprime(2)));
for j1=1:1:length(kvec)
kx = kvec(j1);
for j2=1:1:length(kvec)
ky = kvec(j2);
for j3=1:1:length(kvec)
kxprime = kvec(j3);
for j4=1:1:length(kvec)
kyprime = kvec(j4);
out3 = out3 + f([kx,ky],[kxprime,kyprime]);
end
end
end
end
out3
out3 =
-3000.00000000004
toc
Elapsed time is 0.011396 seconds.
tic
N = 10;
kvec = linspace(-pi,pi,N);
ck = cos(kvec);
ca = reshape(ck, [], 1);
cb = reshape(ck, 1, []);
cc = reshape(ck, 1, 1, []);
cd = reshape(ck, 1, 1, 1, []);
f = bsxfun(@plus, bsxfun(@plus, ca, cb), bsxfun(@plus, cc, cd));
out4 = sum(f(:))
out4 =
-3999.99999999997
toc
Elapsed time is 0.010282 seconds.
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Plus de réponses (2)

Torsten le 17 Sep 2023
Modifié(e) : Torsten le 17 Sep 2023
Note that -pi<=kx<=pi, -pi<=ky<=pi, -pi<=kx'<=pi, -pi<=ky'<=pi because of your linspace choice.
And most probably you need to normalize the sum somehow because at the moment, it depends strongly on N.
sum = 0;
N = 10;
kvec = linspace(-pi,pi,N);
f = @(k,kprime) 2* (cos(k(1)) + cos(k(2))) - 0.5* (cos(kprime(1)) + cos(kprime(2)));
for j1=1:1:length(kvec)
kx = kvec(j1);
for j2=1:1:length(kvec)
ky = kvec(j2);
for j3=1:1:length(kvec)
kxprime = kvec(j3);
for j4=1:1:length(kvec)
kyprime = kvec(j4);
sum = sum + f([kx,ky],[kxprime,kyprime]);
end
end
end
end
sum
sum = -3.0000e+03
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

Dyuman Joshi le 17 Sep 2023
Using 4 nested for loops will be take quite a good amount of time to run, specially if N is a comparetively big value.
You can vectorize your code -
f = @(a,b,c,d) 2*(cos(a) + cos(b)) - 0.5*(cos(c) + cos(d));
%Random value for N for example
N = 10;
kvec = linspace(-pi,pi,N);
[kx,ky,kxprime,kyprime] = ndgrid(kvec);
arr = f(kx,ky,kxprime,kyprime);
Also, it's not a good idea to use inbuilt function names as variables, in your case - sum
out = sum(arr,'all')
out = -3.0000e+03
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Data Type Identification dans Help Center et File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by