原程序代码:
function f=tripleintL1(u,n)
a=3.25*10^(-10);t=1;U=8*t;
x0=0.3551*2*pi/a;
z0=0.18821*2*pi/a;
v0=2*(a^3)/(sqrt(3)*((0.3551)^2)*0.1882);k=1.725*10^(-4);T=10;
f=(3*v0/((2*pi)^3))*(quad(@(x)arrayfun(@(xx)quad2d(@(y,z)(1./(exp((-2*t*(cos(y*a)...
+2*cos(y*a/2).*cos(sqrt(3)*xx*a/2)+cos((sqrt(3)/3)*xx*a+sqrt(6)*z*a/3)...
+2*cos(y*a/2).*cos(sqrt(3)*xx*a/6-sqrt(6)*z*a/3))+U*n-u)./(k*T))+1)...
+1./(exp((-2*t*(cos(y*a)+2*cos(y*a/2).*cos(sqrt(3)*xx*a/2)+cos((sqrt(3)/3)*xx*a+sqrt(6)*z*a/3)...
+2*cos(y*a/2).*cos(sqrt(3)*xx*a/6-sqrt(6)*z*a/3))...
+U*(6-n)-u)./(k*T))+1)),0,sqrt(3)*(xx-x0),0,z0),x),x0,3*x0/2));
end
被积函数是:
1./(exp((-2*t*(cos(y*a)+2*cos(y*a/2).*cos(sqrt(3)*xx*a/2)+cos((sqrt(3)/3)*xx*a+sqrt(6)*z*a/