计算物理基础-数值积分
在很多实际问题中，我们是计算定积分，比如求面积、体积、微分方程等。如果我们知道被积函数的原函数，定积分计算非常容易。然而，多数时候被积函数没有简单的原函数，所以定积分需要借助数值方法。另外，数值方法能通过数目较少的采样积分点，实现对积分的合理估算，比如计算体积的万能公式（通过上、下底面、中间截面）、材料计算中k空间的划分等。
（1）用库函数计算定积分
Octave提供了库函数integral实现数值积分。比如，我们要求半圆的面积，调用函数时先定义被积函数fun = @(x) sqrt(1-x.^2); 然后语法是integral(fun,-1,1)，-1和1分别对应积分区间的左端点、右端点。对于多重积分，我们这里以半球为例，x的范围是[-1,1]，y约束在圆内，所以要定义ymin=@(x) -sqrt(1 - x.^2);ymax=@(x) sqrt(1 - x.^2); 使用函数 integral2 对应二重积分，integral2 (F, -1, 1, ymin, ymax) 后面四个参数分别是x的最小、最大，y的最小、最大值。此外，附录A给出了矩形区域积分案例。
（2）计算定积分的方法对比
运用数值方法计算定积分，通常需要对积分区间进行有限划分，积分方法的优劣取决于积分点数目和精度，优秀的方法可以在较少的积分点情况下，得到较高精度的积分。比如计算体积的万能公式（通过上、下底面、中间截面），只有3个积分点，就可以作为比较好的估计值。这里以半圆面积的计算为例，介绍梯形法、Simpson方法、Gauss方法。
在梯形法中，先把自变量离散化，该面积等价于一系列梯形面积求和（有两个是三角形）。根据梯形面积公式，可知梯形法的积分权重是首尾项为0.5，其他位置是1，线性划分下积分间隔dx就是相邻的x(i+1)-x(i)，这里用的是x(2)-x(1)。在Simpson方法中，我们选择相邻的3个积分点，通过二次插值得到对应的多项式，并对该多项式进行积分，可知积分权重变成：首尾2项的系数是1/3;其它的偶数项的系数是4/3;奇数项的系数是2/3。计算体积的万能公式是Simpson方法的特例。 数值积分关键在于积分点和权重选择，可以根据待定系数法来确定。和梯形法、Simpson方法不同，Gauss方法不要求均匀取点，积分权重也是待定系数。求解系数需要解非线性方程组，在后面的视频讲解。对于多重积分，可以认为是沿着不同方向的一维积分组合而成，所以对应的权重直接相乘即可。详见附录程序B。
如上图所示，不同方法积分精度（纵轴是数值积分结果和真实值的偏差取对数）随积分点数目的变化，Simpson方法比梯形法明显改善，Gauss方法精度更高，但是积分点、权重需要查表。
（3）材料计算中k空间的划分
在材料计算中，物理量的计算需要对k空间进行有限划分，积分转变成不同位置k点对应的求和，而权重和晶格的对称性关联。其中，Monkhorst-Pack 方法可以有效减少k点的数目，以VASP为例，可以通过KPOINTS设置，在IBZKPT文件中查看。另外，对占据数的处理，有四面体方法等，详见视频讲解。参考文献：Phys. Rev. B 13,5188(1976); 49, 16223(1994); Phys. Status Solidi B 54, 469 (1972).
程序附录：
A. 库函数计算定积分
%计算半圆面积
fun = @(x) sqrt(1-x.^2);
integral(fun,-1,1)
%计算半球体积
F = @(x,y) sqrt(1-x.^2-y.^2);
ymin=@(x) -sqrt(1 - x.^2);
ymax=@(x) sqrt(1 - x.^2);
Q = integral2 (F, -1, 1, ymin, ymax)
%矩形区域的二重积分
F = @(X,Y) X.^2.*cos(X)-Y.^2.*sin(Y)+50;
Q = integral2 (F, -5, 5, -5, 5)
B. 定积分不同算法比较
%半圆面积
clear
n=5:2:15;
for i=1:length(n)
x=linspace(-1,1,n(i));
y=sqrt(1-x.^2);
p=ones(1,length(x))/3;
p(2:2:end-1)=4/3;p(3:2:end-1)=2/3;
S_data(i)=sum(y.*p)*(x(2)-x(1));
p=ones(1,length(x));
p(1)=0.5;p(end)=0.5;
T_data(i)=sum(y.*p)*(x(2)-x(1));
end
plot(n,log(abs(S_data-pi/2)),'ro-')
hold on
plot(n,log(abs(T_data-pi/2)),'bo-')
%矩形区域的二重积分
clear
n=5:2:25;
for k=1:length(n)
x=linspace(-5,5,n(k));
y=linspace(-5,5,n(k));
p=ones(1,length(x))/3;
p(2:2:end-1)=4/3;p(3:2:end-1)=2/3;
S_data(k)=0;
for i=1:length(x)
for j=1:length(y)
S_data(k)=S_data(k)+p(i)*p(j)*(x(i)^2*cos(x(i))-y(j)^2*sin(y(i))+50)*(x(2)-x(1))*(y(2)-y(1));
end
end
T_data(k)=0;
p=ones(1,length(x));
p(1)=0.5;p(end)=0.5;
for i=1:length(x)
for j=1:length(y)
T_data(k)=T_data(k)+p(i)*p(j)*(x(i)^2*cos(x(i))-y(j)^2*sin(y(i))+50)*(x(2)-x(1))*(y(2)-y(1));
end
end
end
plot(n,log(abs(S_data-4615.627270748661)),'ro-')
hold on
plot(n,log(abs(T_data-4615.627270748661)),'bo-')
xlabel('log(Num)')
ylabel('log(\Delta y)')
legend('Simpson','Trapezoid')
