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
