? ? ? ?在很多實(shí)際問(wèn)題中，我們是計(jì)算定積分，比如求面積、體積、微分方程等。如果我們知道被積函數(shù)的原函數(shù)，定積分計(jì)算非常容易。然而，多數(shù)時(shí)候被積函數(shù)沒(méi)有簡(jiǎn)單的原函數(shù)，所以定積分需要借助數(shù)值方法。另外，數(shù)值方法能通過(guò)數(shù)目較少的采樣積分點(diǎn)，實(shí)現(xiàn)對(duì)積分的合理估算，比如計(jì)算體積的萬(wàn)能公式（通過(guò)上、下底面、中間截面）、材料計(jì)算中k空間的劃分等。

（1）用庫(kù)函數(shù)計(jì)算定積分

? ? ? ? Octave提供了庫(kù)函數(shù)integral實(shí)現(xiàn)數(shù)值積分。比如，我們要求半圓的面積，調(diào)用函數(shù)時(shí)先定義被積函數(shù)fun = @(x) sqrt(1-x.^2); 然后語(yǔ)法是integral(fun,-1,1)，-1和1分別對(duì)應(yīng)積分區(qū)間的左端點(diǎn)、右端點(diǎn)。對(duì)于多重積分，我們這里以半球?yàn)槔?，x的范圍是[-1,1]，y約束在圓內(nèi)，所以要定義ymin=@(x) -sqrt(1 - x.^2);ymax=@(x) sqrt(1 - x.^2); 使用函數(shù) integral2 對(duì)應(yīng)二重積分，integral2 (F, -1, 1, ymin, ymax) 后面四個(gè)參數(shù)分別是x的最小、最大，y的最小、最大值。此外，附錄A給出了矩形區(qū)域積分案例。

（2）計(jì)算定積分的方法對(duì)比

? ? ? 運(yùn)用數(shù)值方法計(jì)算定積分，通常需要對(duì)積分區(qū)間進(jìn)行有限劃分，積分方法的優(yōu)劣取決于積分點(diǎn)數(shù)目和精度，優(yōu)秀的方法可以在較少的積分點(diǎn)情況下，得到較高精度的積分。比如計(jì)算體積的萬(wàn)能公式（通過(guò)上、下底面、中間截面），只有3個(gè)積分點(diǎn)，就可以作為比較好的估計(jì)值。這里以半圓面積的計(jì)算為例，介紹梯形法、Simpson方法、Gauss方法。

? ? ? ?在梯形法中，先把自變量離散化，該面積等價(jià)于一系列梯形面積求和（有兩個(gè)是三角形）。根據(jù)梯形面積公式，可知梯形法的積分權(quán)重是首尾項(xiàng)為0.5，其他位置是1，線性劃分下積分間隔dx就是相鄰的x(i+1)-x(i)，這里用的是x(2)-x(1)。在Simpson方法中，我們選擇相鄰的3個(gè)積分點(diǎn)，通過(guò)二次插值得到對(duì)應(yīng)的多項(xiàng)式，并對(duì)該多項(xiàng)式進(jìn)行積分，可知積分權(quán)重變成：首尾2項(xiàng)的系數(shù)是1/3;其它的偶數(shù)項(xiàng)的系數(shù)是4/3;奇數(shù)項(xiàng)的系數(shù)是2/3。計(jì)算體積的萬(wàn)能公式是Simpson方法的特例。?數(shù)值積分關(guān)鍵在于積分點(diǎn)和權(quán)重選擇，可以根據(jù)待定系數(shù)法來(lái)確定。和梯形法、Simpson方法不同，Gauss方法不要求均勻取點(diǎn)，積分權(quán)重也是待定系數(shù)。求解系數(shù)需要解非線性方程組，在后面的視頻講解。對(duì)于多重積分，可以認(rèn)為是沿著不同方向的一維積分組合而成，所以對(duì)應(yīng)的權(quán)重直接相乘即可。詳見(jiàn)附錄程序B。

? ? ? ? 如上圖所示，不同方法積分精度（縱軸是數(shù)值積分結(jié)果和真實(shí)值的偏差取對(duì)數(shù)）隨積分點(diǎn)數(shù)目的變化，Simpson方法比梯形法明顯改善，Gauss方法精度更高，但是積分點(diǎn)、權(quán)重需要查表。

（3）材料計(jì)算中k空間的劃分

? ? ? ?在材料計(jì)算中，物理量的計(jì)算需要對(duì)k空間進(jìn)行有限劃分，積分轉(zhuǎn)變成不同位置k點(diǎn)對(duì)應(yīng)的求和，而權(quán)重和晶格的對(duì)稱性關(guān)聯(lián)。其中，Monkhorst-Pack 方法可以有效減少k點(diǎn)的數(shù)目，以VASP為例，可以通過(guò)KPOINTS設(shè)置，在IBZKPT文件中查看。另外，對(duì)占據(jù)數(shù)的處理，有四面體方法等，詳見(jiàn)視頻講解。參考文獻(xiàn)：Phys. Rev. B 13,5188(1976); 49, 16223(1994); Phys. Status Solidi B 54, 469 (1972).

程序附錄：

A. 庫(kù)函數(shù)計(jì)算定積分

%計(jì)算半圓面積

fun = @(x) sqrt(1-x.^2);

integral(fun,-1,1)

%計(jì)算半球體積

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)

%矩形區(qū)域的二重積分

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

%矩形區(qū)域的二重積分

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