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