Python的Scipy模塊中擁有大量的數(shù)值計算函數(shù)，方便我們快速進行數(shù)值計算。
Scipy中的integrate模塊提供了幾種數(shù)值積分算法，導入方式為：
from scipy import integrate
使用integrate時，需要先將要進行積分的方程定義為函數(shù)呻顽。求取一至三重積分的函數(shù)分別為：
integrate.quad(func,a,b,args,full_output)
integrate.dblquad(func,a,b,gfun,hfun,args,epsabs,epsrel)
integrate.tplquad(func,a,b,gfun,hfun,qfun,rfun,args,epsabs,epsrel)
以三重積分為例。func為運算對象函數(shù)丹墨，形式為func(z,y,x)廊遍。a,b對應變量x的積分區(qū)域，gfun,hfun對應變量y的積分區(qū)域贩挣，依次類推喉前。
注意gfun,hfun等的形式應為函數(shù)，其中gfun,hfun是自變量為x的函數(shù)王财，qfun,rfun是自變量為x,y的函數(shù)卵迂。這些函數(shù)可以使用lambda函數(shù)進行定義，形式通常為：
lambda x,y:x*y
如果是常函數(shù)绒净，則定義為：
lambda x:0
lambda x,y:1
args可選见咒，為傳遞給func的格外參數(shù)；full_output可選挂疆，非零則返回積分信息的dictionary改览。如果非零，則還會禁止顯示警告消息缤言，并將消息附加到輸出元組宝当。epsabs可選，絕對容差直接傳遞到內(nèi)部1-D正交積分胆萧。默認值為1.49e-8庆揩。epsrel可選，內(nèi)部1-D積分的相對容差鸳碧。默認值為1.49e-8盾鳞。
運算結果輸出一個元組，第一項為積分結果瞻离，第二項為絕對誤差腾仅，此外還有收斂情況等信息。
以下是一個用integrate.tplquad()計算立方體晶體的LJ勢能的例子：
#To calculate the LJ potential energy between a cubic and an atom
#Author: Lewisbase
#Date: 2019.01.01
from __future__ import division
import numpy as np
from scipy import integrate
#Cubic length width hight A,B,C
A=10
B=10
C=10
#Average sigma(s,A) epsilon(e,K) rho(r,A^-3) of the cubic
s=np.array([4.24,3.84]); e=np.array([208.4541197,69.0609]); r=0.061540355
def LJ_cubic(x,y,z):
return 4*eps*rho*((sig**2/((xb-x)**2+(yb-y)**2+(zb-z)**2))**6-(sig**2/((xb-x)**2+(yb-y)**2+(zb-z)**2))**3)
lx=50
ly=50
lz=100
dx=lx/101
dy=ly/101
dz=lz/201
x0=lx/2-A/2
x1=lx/2+A/2
y0=ly/2-B/2
y1=ly/2+B/2
z0=lz/2-C/2
z1=lz/2+C/2
for m in np.linspace(0,1,2):
for ix in np.linspace(0,100,101):
for iy in np.linspace(0,100,101):
for iz in np.linspace(0,200,201):
xb=ix*dx
yb=iy*dy
zb=iz*dz
sig=s[m]
eps=e[m]
rho=r
if abs(xb-lx/2)<(A/2) and abs(yb-ly/2)<(B/2) and abs(zb-lz/2)<(C/2) :
LJ = 1E100
else:
LJ=integrate.tplquad(LJ_cubic,
z0,
z1,
lambda y: y0,
lambda y: y1,
lambda y,x: x0,
lambda y,x: x1)
LJ=LJ[0]/298.15
print "%f %f %f %f\n"%(xb,yb,zb,LJ)
有關數(shù)值積分的方法還有高斯積分套利，龍貝格積分等推励，自己暫時還沒有用到，待到以后再做詳細學習肉迫。
參考資料
用Python作科學計算
Scipy.Integrate
Scipy.Integrate.quad
Scipy.Integrate.dblquad
Scipy.Integrate.tplquad
