统计计算 课件 3.4 积分的计算

3.4积分的计算3.4.1使用Python计算积分3.4.2随机模拟求积分3.4.3梯形法和辛普森法3.4.1使用Python计算积分使用Python的Scipy库中的函数可以进行积分的计算。quad计算定积分，具体格式如下：quad函数返回两个值，第一个值是积分的值，第二个值是对积分值的绝对误差估计。二重积分使用函数dblquad（doublequad），具体格式如下：dblquad(func,a,

b,gfun,hfun,args=(),epsabs=1.49e-08,epsrel=1.49e-08)，被积函数func的格式为func(y,x)，先对y积分，再对x积分。a和b为x的积分限，gfun和hfun为y的积分限。gfun是y的下边界曲线，它是一个函数，接受单个浮点数参数(x)，并返回一个浮点数结果或一个浮点数，表示一个常量边界曲线。hfun是y中的上边界曲线。要是积分上下限为x或y的函数，要用lambda来表示上下限。三重积分使用函数tplquad，调用格式为tplquad(func,a,b,gfun,hfun,qfun,rfun,args=(),epsabs=1.49e-08,epsrel=1.49e-08)被积函数func的格式为func(z,y,x)，先对z积分，再对y积分，最后对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随机模拟求积分算法：（1）使用Python中的dblquad计算二重积分的解；（2）取服从区间[-1,1]的均匀分布的随机数x和y，再取服从区间[0,1]的均匀分布随机数z；（3）取Ω为长为2，宽为2，高度为1的长方体，在Ω中产生服从均匀分布的N个点；（4）若满足如下条件，计数s+1，总的计数为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计算二重积分的解；（2）取服从区间[-1,1]的均匀分布的随机数x和y，再取服从区间[0,1]的均匀分布随机数z；（3）四维超长方体。在Ω中产生服从均匀分布的N个点；（4）若满足如下条件，计数s+1，总的计数为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输出结果：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、四维圆柱体

与四维超立方体相比，把区域Ω取为四维圆柱体。只需修改（2）取Ω为的四维圆柱体。在Ω中产生服从均匀分布的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使用四维超长方体效果比四维圆柱体好一些梯形法均有下列式子准确成立对于代数精确度指的是某个求积分公式对次数不超过m的多项式均能准确成立，对m+1次多项式不成立，称该求积分公式具有m次精度。两者不相等，因此该积分公式有1次精度例7使用复合梯形公式计算积分fromegrateimportquadimportnumpyasnp#给出积分区间a,b,小区间的个数na=0b=1n=1000#将区间[a,b]n等分后取小区间的n+1个端点值,小区间的长度hX=np.linspace(a,b,n+1)h=(b-a)/n#函数ff=lambdax:np.exp(-(x**2)/2)/np.sqrt(2*np.pi)#使用梯形公式计算积分I7=(2*np.sum(f(X))-f(a)-f(b))*h/2print(I7)#使用Python的quad计算例7定积分I,r=quad(f,0,1)print(I,'使用梯形法的相对误差',abs(I7-I)/I)输出结果：0.34134472590431530.341344746068543使用梯形法的相对误差5.9072910623996923e-08辛普森法将积分区间[a,b]等分为2n份，辛普森公式使用抛物线近似代替曲线，设将该曲边梯形根据小区间分为2n个小曲边梯形，则定积分的值等于n个小曲边梯形面积的和。对于每个小曲边梯形，将其曲边近似为抛物线边，计算梯形的面积。说明对于x3成立，经验证发现对x4不成立，因此精度为3.复合辛普森公式余项为例8使用复合辛普森法计算积分fromegrateimportquadimportnumpyasnp#给出积分区间a,b,小区间的个数na=0b=1n=1000#将区间[a,b]2n等分后取小区间的n+1个端点值,小区间的长度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,'使