統(tǒng)計(jì)計(jì)算 課件 3.4 積分的計(jì)算

3.4積分的計(jì)算3.4.1使用Python計(jì)算積分3.4.2隨機(jī)模擬求積分3.4.3梯形法和辛普森法3.4.1使用Python計(jì)算積分使用Python的Scipy庫中的函數(shù)可以進(jìn)行積分的計(jì)算。quad計(jì)算定積分，具體格式如下：quad函數(shù)返回兩個(gè)值，第一個(gè)值是積分的值，第二個(gè)值是對(duì)積分值的絕對(duì)誤差估計(jì)。二重積分使用函數(shù)dblquad（doublequad），具體格式如下：dblquad(func,a,

b,gfun,hfun,args=(),epsabs=1.49e-08,epsrel=1.49e-08)，被積函數(shù)func的格式為func(y,x)，先對(duì)y積分，再對(duì)x積分。a和b為x的積分限，gfun和hfun為y的積分限。gfun是y的下邊界曲線，它是一個(gè)函數(shù)，接受單個(gè)浮點(diǎn)數(shù)參數(shù)(x)，并返回一個(gè)浮點(diǎn)數(shù)結(jié)果或一個(gè)浮點(diǎn)數(shù)，表示一個(gè)常量邊界曲線。hfun是y中的上邊界曲線。要是積分上下限為x或y的函數(shù)，要用lambda來表示上下限。三重積分使用函數(shù)tplquad，調(diào)用格式為tplquad(func,a,b,gfun,hfun,qfun,rfun,args=(),epsabs=1.49e-08,epsrel=1.49e-08)被積函數(shù)func的格式為func(z,y,x)，先對(duì)z積分，再對(duì)y積分，最后對(duì)x積分。a和b為x的積分限，gfun和hfun為y的積分限，qfun和rfun為z的積分限。fromegrateimportquadimportnumpyasnpz1=lambdax:np.exp(x)I1=quad(z1,0,2)I=np.exp(2)-1積分值為6.3890560989306495，使用積分公式得到的積分的精確值為6.38905609893，fromegrateimportdblquad#例2二重積分1先y后xz2=lambday,x:x**2+y**2I2=dblquad(z2,0,1,lambdax:0,lambdax:1)print(I2)積分值為0.66667。fromegrateimportdblquadimportnumpyasnp#例3二重積分2先y后xz31=lambday,x:np.exp(-(y**2))I31=dblquad(z31,0,2,lambdax:x,lambdax:2)print(I31)#例3二重積分2先x后yz32=lambdax,y:np.exp(-(y**2))I32=dblquad(z32,0,2,lambday:0,lambday:y)print(I32,(1-np.exp(-4))/2)0.4908421805556329fromegrateimporttplquad#三重積分先z后y再xz4=lambdaz,y,x:x+y+zI4=tplquad(z4,0,1,lambdax:0,lambdax:1-x,lambdax,y:0,lambdax,y:1-x-y)print(I4)0.1253.4.2隨機(jī)模擬求積分算法：（1）使用Python中的dblquad計(jì)算二重積分的解；（2）取服從區(qū)間[-1,1]的均勻分布的隨機(jī)數(shù)x和y，再取服從區(qū)間[0,1]的均勻分布隨機(jī)數(shù)z；（3）取Ω為長為2，寬為2，高度為1的長方體，在Ω中產(chǎn)生服從均勻分布的N個(gè)點(diǎn)；（4）若滿足如下條件，計(jì)數(shù)s+1，總的計(jì)數(shù)為sfromegrateimportdblquadimportnumpyasnp#例5二重積分2先y后xz5=lambday,x:np.cos(np.sqrt(x**2+y**2))I5,r5=dblquad(z5,-1,1,lambdax:-np.sqrt(1-x**2),lambdax:np.sqrt(1-x**2))print(I5)N=20000s=0foriinrange(N):x=np.random.uniform(-1,1)y=np.random.uniform(-1,1)z=np.random.uniform(0,1)ifx**2+y**2<=1andz<=np.cos(np.sqrt(x**2+y**2))andz>=0:s+=1I=4*s/Nprint(I,2*np.pi*(np.sin(1)+np.cos(1)-1))2.3987523306493152.39742.3987523306492724算法：（1）使用Python中的dblquad計(jì)算二重積分的解；（2）取服從區(qū)間[-1,1]的均勻分布的隨機(jī)數(shù)x和y，再取服從區(qū)間[0,1]的均勻分布隨機(jī)數(shù)z；（3）四維超長方體。在Ω中產(chǎn)生服從均勻分布的N個(gè)點(diǎn)；（4）若滿足如下條件，計(jì)數(shù)s+1，總的計(jì)數(shù)為sfromegrateimporttplquadimportnumpyasnp#例6三重積分先z再y后x,使用四維長方體z6=lambdaz,y,x:2*np.sqrt(x**2+y**2)I61,r6=tplquad(z6,0,1,lambdax:0,lambdax:np.sqrt(1-x**2),lambdax,y:0,lambdax,y:np.sqrt(1-x**2-y**2))print(I61)

N=100000s=0輸出結(jié)果：0.61685027506797570.616540.6168502750680849foriinrange(N):x=np.random.uniform(0,1)y=np.random.uniform(0,1)z=np.random.uniform(0,1)g=np.random.uniform(0,2)ifx**2+y**2<=1and0<=z<=np.sqrt(1-x**2-y**2)andg<=2*np.sqrt(x**2+y**2):s+=1I62=2*s/Nprint(I62,np.pi**2/16)2、四維圓柱體

與四維超立方體相比，把區(qū)域Ω取為四維圓柱體。只需修改（2）取Ω為的四維圓柱體。在Ω中產(chǎn)生服從均勻分布的N個(gè)點(diǎn)，四維圓柱體的體積為Π/2fromegrateimportdblquad,tplquadimportnumpyasnpN=20000s=0foriinrange(N):x=np.random.uniform(0,1)y=np.random.uniform(0,np.sqrt(1-x**2))z=np.random.uniform(0,1)g=np.random.uniform(0,2)ifx**2+y**2<=1and0<=z<=np.sqrt(1-x**2-y**2)andg<=2*np.sqrt(x**2+y**2):s+=1I63=(np.pi/2)*s/Nprint(I63)0.6024946391054505使用四維超長方體效果比四維圓柱體好一些梯形法均有下列式子準(zhǔn)確成立對(duì)于代數(shù)精確度指的是某個(gè)求積分公式對(duì)次數(shù)不超過m的多項(xiàng)式均能準(zhǔn)確成立，對(duì)m+1次多項(xiàng)式不成立，稱該求積分公式具有m次精度。兩者不相等，因此該積分公式有1次精度例7使用復(fù)合梯形公式計(jì)算積分fromegrateimportquadimportnumpyasnp#給出積分區(qū)間a,b,小區(qū)間的個(gè)數(shù)na=0b=1n=1000#將區(qū)間[a,b]n等分后取小區(qū)間的n+1個(gè)端點(diǎn)值,小區(qū)間的長度hX=np.linspace(a,b,n+1)h=(b-a)/n#函數(shù)ff=lambdax:np.exp(-(x**2)/2)/np.sqrt(2*np.pi)#使用梯形公式計(jì)算積分I7=(2*np.sum(f(X))-f(a)-f(b))*h/2print(I7)#使用Python的quad計(jì)算例7定積分I,r=quad(f,0,1)print(I,'使用梯形法的相對(duì)誤差',abs(I7-I)/I)輸出結(jié)果：0.34134472590431530.341344746068543使用梯形法的相對(duì)誤差5.9072910623996923e-08辛普森法將積分區(qū)間[a,b]等分為2n份，辛普森公式使用拋物線近似代替曲線，設(shè)將該曲邊梯形根據(jù)小區(qū)間分為2n個(gè)小曲邊梯形，則定積分的值等于n個(gè)小曲邊梯形面積的和。對(duì)于每個(gè)小曲邊梯形，將其曲邊近似為拋物線邊，計(jì)算梯形的面積。說明對(duì)于x3成立，經(jīng)驗(yàn)證發(fā)現(xiàn)對(duì)x4不成立，因此精度為3.復(fù)合辛普森公式余項(xiàng)為例8使用復(fù)合辛普森法計(jì)算積分fromegrateimportquadimportnumpyasnp#給出積分區(qū)間a,b,小區(qū)間的個(gè)數(shù)na=0b=1n=1000#將區(qū)間[a,b]2n等分后取小區(qū)間的n+1個(gè)端點(diǎn)值,小區(qū)間的長度hX=np.linspace(a,b,n+1)h=(b-a)/nf=lambdax:np.exp(-(x**2)/2)/np.sqrt(2*np.pi)x1=[]x2=[]foriinrange(n+1):ifi%2==1:t=f(X[i])x1.append(t)else:t=f(X[i])x2.append(t)I8=(4*np.sum(x1)+2*np.sum(x2)-f(a)-f(b))*h/3print(I8)I,r=quad(f,0,1)print(I,'使