Function over a vector which involves a sum over 3D coordinates

2022 9 1
2 个回答
回答已采纳
更新时间：2022 9 1
4 次查看（30 天）

0 个投票

Hi all, I am working on a numerical approximation for heat capacity as a function of temperature. The temperate is a 1x500 vector, and so the heat capacity vector needs to be 1x500 as well. My problem is that the function at each temperature involves a sum over a 100x100x100 grid. EDIT: The function in question is:
Where the vector runs over the 100x100x100 values in the grid. The question is how can I deal with creating the discrete values as a function of T with the sum over the grid in each entry? The code I tried to run is:
%%Heat Capacity
%Definining relevant quanitites
k_B = 1.38*10^(-23);
N = 100;
a = 4*10^(-10); %m
L_i = a*N;
w_0 = 10*10^12; %Hz
hbar = 1.0545*10^(-34);
c_v_a_ht = 3*k_B*N^3/a^3;
del_k = (2*pi)^3/L_i^3;
del_k_i = 2*pi/L_i;
bz_1 = pi/a;
%First Brillouin Zone
k_x = -bz_1:del_k_i:bz_1-del_k_i;
k_y = -bz_1:del_k_i:bz_1-del_k_i;
k_z = -bz_1:del_k_i:bz_1-del_k_i;
[X,Y,Z] = meshgrid(k_x,k_y,k_z);
%Temp vector
debye_temp = (hbar/k_B)*w_0*(3*pi^2/4)^(1/3); %K
T = 0.001*debye_temp:debye_temp/500:debye_temp;
%Heat capacity approx
abs_k = sqrt(X.^2+Y.^2+Z.^2);
beta = 1./(k_B*T);
c_v_a = ones(1,500);
for i=1:500
ex_k = exp(beta(i)*hbar*w_0*a.*abs_k/2)
summand = (hbar^2* w_0^2 .*(abs_k.^2)*a^2.*ex_k)./(4*(ex_k-1).^2*k_B*(T(i)^2))
c_v_a(i) = sum(summand,'all')
end

11 个评论
显示 9更早的评论 隐藏 9更早的评论

Star Strider
Star Strider 2022-9-1
It might be best to interpolate the temperature vector to (1x100) to match the matrix. Doing the opposite (interpolating the matrix to be 500 in the chosen dimension) risks creating data.
Jared Goldberg
Jared Goldberg 2022-9-1
I edited the original post to include the function in question and the code I tried to run. I don't think the size of the matrix is the issue here, I'm trying to sum over all of its values in a function for each value of T...
Torsten
Torsten 2022-9-1
Did you make sure that ex_k is not 1 and T(i) is not 0 ?
Jared Goldberg
Jared Goldberg 2022-9-1
T(i) isn't zero. ex_k shouldn't obtain values of one as this is a low temp approximation and so the exponent should be large. without a snag of it obtaining unity values, do you think this code should work? is there a more efficient way to build the vector for the heat capacity?
Torsten
Torsten 2022-9-1
If I exclude the denominator in your expression for "summand", I get no NaN values ...
Jared Goldberg
Jared Goldberg 2022-9-1
I also don't, but I can't just exclude the denominator, it is part of the function I am trying to approximate.
Bruno Luong
Bruno Luong 2022-9-1
编辑：Torsten 2022-9-1
Ran in:
Ran in:
>See what MATLAB server returns:
%%Heat Capacity
%Definining relevant quanitites
k_B = 1.38*10^(-23);
N = 100;
a = 4*10^(-10); %m
L_i = a*N;
w_0 = 10*10^12; %Hz
hbar = 1.0545*10^(-34);
c_v_a_ht = 3*k_B*N^3/a^3;
del_k = (2*pi)^3/L_i^3;
del_k_i = 2*pi/L_i;
bz_1 = pi/a;
%First Brillouin Zone
k_x = -bz_1:del_k_i:bz_1-del_k_i;
k_y = -bz_1:del_k_i:bz_1-del_k_i;
k_z = -bz_1:del_k_i:bz_1-del_k_i;
[X,Y,Z] = meshgrid(k_x,k_y,k_z);
%Temp vector
debye_temp = (hbar/k_B)*w_0*(3*pi^2/4)^(1/3); %K
T = 0.001*debye_temp:debye_temp/500:debye_temp;
%Heat capacity approx
abs_k = sqrt(X.^2+Y.^2+Z.^2);
beta = 1./(k_B*T);
c_v_a = ones(1,500);
for i=1:500
ex_k = exp(beta(i)*hbar*w_0*a.*abs_k/2);
if any(ex_k(:) == 1)
fprintf('0 denominator for sure\n')
return
end
summand = (hbar^2* w_0^2 .*(abs_k.^2)*a^2.*ex_k)./(4*(ex_k-1).^2*k_B*(T(i)^2));
c_v_a(i) = sum(summand,'all');
end
0 denominator for sure
Torsten
Torsten 2022-9-1
I also don't, but I can't just exclude the denominator, it is part of the function I am trying to approximate.
But that's no argument to divide by 0 ...
Jared Goldberg
Jared Goldberg 2022-9-1
it also wasn't division by zero, as the exponent is rather large in the low temp limit...

请先登录，再进行评论。

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

 采纳的回答

Jared Goldberg
Jared Goldberg 2022-9-1

0 个投票

Found my problem, had to use three sums over the summand variable instead of using sum(summand, all). This was the proper way to preform the sum. Thanks all!

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Torsten
Torsten 2022-9-1
This cannot solve the NaN values in "summand", but if you are happy with the result, we are too.
Bruno Luong
Bruno Luong 2022-9-1
This explanation doesn't make sense.
Mistake on mistake...
Jared Goldberg
Jared Goldberg 2022-9-1
编辑：Jared Goldberg 2022-9-1
To clarify, forgot to write it here I wrote it below: the NaN values in the summand were due to an erroneous *a I corrected in the edit. It still didn't work after that, though, but this solved the problem. Thanks so much, really appreciate the help in finding the errors!
Currently it's running, we'll see how the result comes out!
Bruno Luong
Bruno Luong 2022-9-1
编辑：Bruno Luong 2022-9-1
So the "solution" is not sum(sum(sum(...))).
This answer is non-sense.
By changing a you just make you grid differently and by chance it does not contain (0,0,0).

请先登录，再进行评论。

更多回答（1 个）

Bruno Luong
Bruno Luong 2022-9-1
编辑：Bruno Luong 2022-9-1

0 个投票

You get 0/0 expression when abs_k is small (or 0), returning NaN.
summand = (hbar^2* w_0^2 .*(abs_k.^2)*a^2.*ex_k)./(4*(ex_k-1).^2*k_B*(T(i)^2))
You should make use of Taylor expansion of exp function to simplify the ratio
... (abs_k).^2/(ex_k-1) ...

7 个评论
显示 5更早的评论 隐藏 5更早的评论

Jared Goldberg
Jared Goldberg 2022-9-1
abs_k should never be small or zero. they should all be relatively large values, I just checked the code to make sure of this and they are...
also - the whole point of the exercise is to evaluate the function numerically to check the validity of the Debye approximation, so approximating the function is off the table.
Bruno Luong
Bruno Luong 2022-9-1
编辑：Bruno Luong 2022-9-1
Ran in:
Ran in:
"abs_k should never be small or zero"
Not in your script.
How do you check? See what the MATLAB server tell me, there is clearly index where abs_k == 0
%%Heat Capacity
%Definining relevant quanitites
k_B = 1.38*10^(-23);
N = 100;
a = 4*10^(-10); %m
L_i = a*N;
w_0 = 10*10^12; %Hz
hbar = 1.0545*10^(-34);
c_v_a_ht = 3*k_B*N^3/a^3;
del_k = (2*pi)^3/L_i^3;
del_k_i = 2*pi/L_i;
bz_1 = pi/a;
%First Brillouin Zone
k_x = -bz_1:del_k_i:bz_1-del_k_i;
k_y = -bz_1:del_k_i:bz_1-del_k_i;
k_z = -bz_1:del_k_i:bz_1-del_k_i;
[X,Y,Z] = meshgrid(k_x,k_y,k_z);
%Temp vector
debye_temp = (hbar/k_B)*w_0*(3*pi^2/4)^(1/3); %K
T = 0.001*debye_temp:debye_temp/500:debye_temp;
%Heat capacity approx
abs_k = sqrt(X.^2+Y.^2+Z.^2);
beta = 1./(k_B*T);
c_v_a = ones(1,500);
for i=1:500
ex_k = exp(beta(i)*hbar*w_0*a.*abs_k/2);
% Detect pb in denominator
b = ex_k == 1;
if any(b(:))
ibad = find(b(:));
abs_k = abs_k(ibad)
beta(i)*hbar*w_0*a.*abs_k/2
exp(beta(i)*hbar*w_0*a.*abs_k/2)
return
end
summand = (hbar^2* w_0^2 .*(abs_k.^2)*a^2.*ex_k)./(4*(ex_k-1).^2*k_B*(T(i)^2));
c_v_a(i) = sum(summand,'all');
end
abs_k = 0
ans = 0
ans = 1
Bruno Luong
Bruno Luong 2022-9-1
Ran in:
Ran in:
%Definining relevant quanitites
k_B = 1.38*10^(-23);
N = 100;
a = 4*10^(-10); %m
L_i = a*N;
w_0 = 10*10^12; %Hz
hbar = 1.0545*10^(-34);
c_v_a_ht = 3*k_B*N^3/a^3;
del_k = (2*pi)^3/L_i^3;
del_k_i = 2*pi/L_i;
bz_1 = pi/a;
%First Brillouin Zone
k_x = -bz_1:del_k_i:bz_1-del_k_i;
k_y = -bz_1:del_k_i:bz_1-del_k_i;
k_z = -bz_1:del_k_i:bz_1-del_k_i;
[X,Y,Z] = meshgrid(k_x,k_y,k_z);
abs_k = sqrt(X.^2+Y.^2+Z.^2);
find(abs_k==0)
ans = 505051
Jared Goldberg
Jared Goldberg 2022-9-1
check the code again... there was an earlier edit which had an accidental a typed twice.
Jared Goldberg
Jared Goldberg 2022-9-1
we shall see, right now it's running and all is well

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 MATLAB 的更多信息

产品

版本

R2020a

标签

Community Treasure Hunt

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

Start Hunting!

Translated by