高等数学问题的数值解(3)--数值积分

一重积分

例:分别使用梯形公式、辛普森公式和SciPy工具库中函数quad来求积分$\int_{0}^{1}sin(\sqrt{cos(x)+x^{2}})dx$的数值解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import numpy as np
from scipy.integrate import quad

def trapezoid(f,n,a,b): # 定义梯形公式的函数
xi = np.linspace(a,b,n)
h = (b-a)/(n-1)
return h*(np.sum(f(xi))-(f(a)+f(b))/2)

def simpson(f,n,a,b): # 定义辛普森公式的函数
xi = np.linspace(a,b,2*n+1)
h = (b-a)/(2.0*n)
xe = [f(xi[i]) for i in range(len(xi)) if i%2==0]
xo = [f(xi[i]) for i in range(len(xi)) if i%2!=0]
return h*(2*np.sum(xe)+4*np.sum(xo)-f(a)-f(b))/3.0

a = 0
b = 1
n = 1000
f = lambda x:np.sin(np.sqrt(np.cos(x)+x**2))

print('梯形积分I1 = ', trapezoid(f,n,a,b))
print('辛普森积分I2 = ', simpson(f,n,a,b))
print('SciPy积分I3 = ', quad(f,a,b))

多重积分

例:计算$\int_{0}^{2}dx\int_{0}^{1}xy^{2}dy)$的数值解

1
2
3
4
5
6
import numpy as np
from scipy.integrate import dblquad

f = lambda y,x: x*y**2 # 变量的顺序不能反(先里再外)

print('I:', dblquad(f, 0,2,0,1))

例:计算$I = \iint\limits_{x^{2}+y^{2}\leq 1} e^{-\frac{x^{2}}{2}}sin(x^{2}+y)dxdy$的数值解
先把二重积分化成累次积分
$$I = \int_{-1}^{1}dx \int_{-\sqrt{1-x^{2}}}^{\sqrt{1-x^{2}}} e^{-\frac{x^{2}}{2}}sin(x^{2}+y)dy$$

1
2
3
4
5
6
7
import numpy as np
from scipy.integrate import dblquad

f = lambda y,x:np.exp(-x**2/2)*np.sin(x**2+y)
bd = lambda x:np.sqrt(1-x**2)

print('I:', dblquad(f, -1,1,lambda x:-bd(x),bd))

例:计算$\iiint\limits_{\Omega } z \sqrt{x^{2}+y^{2}+1} dxdydz$ ,其中$\Omega $为柱面$x^{2}+y^{2}-2x=0$与$z=0$,$z=6$两平面所围成的空间区域。
先把三重积分化成累次积分
$$I = \int_{0}^{2}dx\int_{-\sqrt{2x-x^{2}}}^{\sqrt{2x-x^{2}}}dy\int_{0}^{6}z\sqrt{x^{2}+y^{2}+1}dz$$

1
2
3
4
5
6
7
import numpy as np
from scipy.integrate import tplquad

f = lambda z,y,x: z*np.sqrt(x**2+y**2+1)
ybd = lambda x:np.sqrt(2*x-x**2)

print('I=', tplquad(f,0,2,lambda x:-ybd(x),ybd,0,6))
0%