# 12-04：Gauss求积公式

``````import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
``````

# Part 1:Gauss-Legendre求积公式

``````def f(x):
return ((10/x)**2)*np.sin(10/x)
``````
``````# 真实值
from scipy import integrate
def f(x):
return ((10/x)*8*2)*np.sin(10/x)
v, err = integrate.quad(f, 1, 3)
print (v,err)
``````
``````-1.426024756346262 1.674207356451877e-08
``````

#### 两点

``````a,b = 1,3
A0 = (b-a)/2

t = np.array([-0.5773502692,0.5773502692])
A = np.array([1,1])

X = (a+b)/2+((b-a)/2)*t

F = np.sum(A*A0*f(X))
print(F)
``````
``````23.399780052805777
``````
``````def Two(a,b):
t = np.array([-0.5773502692,0.5773502692])
A = np.array([1,1])
X = (a+b)/2+((b-a)/2)*t
F = (b-a)*np.sum(A*f(X))/2
return F
a,b = 1,3
Two(a,b)
``````
``````23.399780052805777
``````

#### 三点

``````t = np.array([-0.7745966692,0,0.7745966692])
A = np.array([5/9,8/9,5/9])

X = (a+b)/2+((b-a)/2)*t

F = np.sum(A*A0*f(X))
print(F)
``````
``````10.74235460649307
``````
``````def Three(a,b):
t = np.array([-0.7745966692,0,0.7745966692])
A = np.array([5/9,8/9,5/9])
X = (a+b)/2+((b-a)/2)*t
F = (b-a)*np.sum(A*f(X))/2
return F
a,b = 1,3
Three(a,b)
``````
``````10.74235460649307
``````

#### 四点

``````t = np.array([-0.8611363116,-0.3399810436,0.3399810436,0.8611363116])
A = np.array([0.3478548451,0.6521451549,0.6521451549,0.3478548451])

X = (a+b)/2+((b-a)/2)*t

F = np.sum(A*A0*f(X))
print(F)
``````
``````-2.212858309779213
``````

#### 五点

``````t = np.array([-0.9061798459,-0.5384693101,0,0.5384693101,0.9061798459])
A = np.array([0.2369268851,0.4786286705,0.568888889,0.4786286705,0.2369268851])

X = (a+b)/2+((b-a)/2)*t

F = np.sum(A*A0*f(X))
print(F)
``````
``````-2.379448453092436
``````

#### 七点

``````t = np.array([-0.9602898565,-0.79666647774,-0.5255324099,-0.1834346425,0.1834346425,0.5255324099,0.79666647774,0.9602898565])
A = np.array([0.1012285361,0.2223810345,0.3137066459,0.3626837834,0.3626837834,0.3137066459,0.2223810345,0.1012285361])

X = (a+b)/2+((b-a)/2)*t

F = np.sum(A*A0*f(X))
print(F)
``````
``````-1.4206298545026215
``````
``````def Seven(a,b):
t = np.array([-0.9602898565,-0.79666647774,-0.5255324099,-0.1834346425,0.1834346425,0.5255324099,0.79666647774,0.9602898565])
A = np.array([0.1012285361,0.2223810345,0.3137066459,0.3626837834,0.3626837834,0.3137066459,0.2223810345,0.1012285361])
X = (a+b)/2+((b-a)/2)*t
F = (b-a)*np.sum(A*f(X))/2
return F
a,b = 1,3
Seven(a,b)
``````
``````-1.4206298545026215
``````

# Part 2 :复合Gauss-Legendre求积公式

#### 分段二阶

``````a,b,n = 1,3,100
x = np.linspace(a,b,n)
re = []
for i in range(len(x)-1):
xi = x[i]
xii = x[i+1]

t = np.array([-0.5773502692,0.5773502692])
A = np.array([1,1])

X = (xi+xii)/2+((xii-xi)/2)*t
F = (xii-xi)*np.sum(A*f(X))/2
re.append(F)

print(np.sum(re))
``````
``````-1.426029291763135
``````
``````def TwoPlus(a,b,n):
x = np.linspace(a,b,n)
re = []
for i in range(len(x)-1):
xi = x[i]
xii = x[i+1]
F = Two(xi,xii)
re.append(F)
return np.sum(re)
a,b,n = 1,3,100
TwoPlus(a,b,n)
``````
``````-1.426029291763135
``````

#### 分段三阶

``````a,b,n = 1,3,10
x = np.linspace(a,b,n)
re = []
for i in range(len(x)-1):
xi = x[i]
xii = x[i+1]

t = np.array([-0.7745966692,0,0.7745966692])
A = np.array([5/9,8/9,5/9])

X = (xi+xii)/2+((xii-xi)/2)*t
F = (xii-xi)*np.sum(A*f(X))/2
re.append(F)

print(np.sum(re))
``````
``````-1.4268535580329276
``````
``````def ThreePlus(a,b,n):
x = np.linspace(a,b,n)
re = []
for i in range(len(x)-1):
xi = x[i]
xii = x[i+1]
F = Three(xi,xii)
re.append(F)
return np.sum(re)
a,b,n = 1,3,20
ThreePlus(a,b,n)
``````
``````-1.426030077865446
``````

#### 分段七阶

``````def SevenPlus(a,b,n):
x = np.linspace(a,b,n)
re = []
for i in range(len(x)-1):
xi = x[i]
xii = x[i+1]
F = Seven(xi,xii)
re.append(F)
return np.sum(re)
a,b,n = 1,3,5
SevenPlus(a,b,n)
``````
``````-1.4260247638153247
``````

#### 复合三阶随分段数变化的积分值变化

``````n = np.arange(5,55,2)
a,b = 1,3
result = []
for i in n:
res = ThreePlus(a,b,i)
result.append(res)
result = np.array(result)

plt.figure(figsize=(8,6))
plt.plot(n,result-v)
plt.scatter(n,result-v)
plt.xlabel('分段数')
plt.ylabel('误差')
plt.title('复合三阶Gauss-Legendre求积误差与分段数关系')
plt.grid()
``````
image.png