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

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

[KSSV]

Show us your code.

[Jared Goldberg]

Will do once I'm home, thanks.

[Star Strider]

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]

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]

Did you make sure that ex_k is not 1 and T(i) is not 0 ?

[Jared Goldberg]

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]

If I exclude the denominator in your expression for "summand", I get no NaN values ...

[Jared Goldberg]

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]

>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]

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]

it also wasn't division by zero, as the exponent is rather large in the low temp limit...

Answer 1

[Jared Goldberg]

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!

[Torsten]

This cannot solve the NaN values in "summand", but if you are happy with the result, we are too.

[Bruno Luong]

This explanation doesn't make sense.

Mistake on mistake...

[Jared Goldberg]

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]

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).

Answer 2

[Bruno Luong]

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) ...

[Jared Goldberg]

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]

"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

[Jared Goldberg]

K>> abs_k(1,1,1)

ans =

1.3603e+10

[Bruno Luong]

%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]

check the code again... there was an earlier edit which had an accidental a typed twice.

[Bruno Luong]

>> abs_k(505051)

ans =

0

[Jared Goldberg]

we shall see, right now it's running and all is well