理论知识:
勒让德多项式及性质
第四节 高斯(Gauss)求积公式
高斯-勒让德求积公式及Matlab实现
目标函数求积区间为[-1,1]时
代码:
#高斯-勒让德求积公式
from sympy import *
from scipy.special import perm,comb #排列,组合
n = 2 #n次多项式正交,n越大精度越高(n=0,1,2,...)
x = symbols("x")
#勒让德多项式
def L(n):
df = diff((x ** 2 - 1) ** (n + 1), x, n + 1)
# Python内置阶乘函数factorial
# L = factorial(n+1) / factorial(2*n+2) * df
L = 1 /2**(n+1)/factorial(n+1) * df
return L
#高斯点x求取
def Gauss_point(n):
x_k_list = solve(L(n)) #求得零点解集
return x_k_list
#求积系数A
def Quadrature_coefficient(x_k_list):
A_list = []
for x_k in x_k_list:
A = 2/(1-x_k**2)/(diff(L(n),x,1).subs(x,x_k))**2
A_list.append(A)
return A_list
#需要求积的目标函数
def f(x):
f = sqrt(x+1.5)
return f
result = 0
x_k_list = Gauss_point(n)
A_list = Quadrature_coefficient(Gauss_point(n))
sum = len(A_list)
for i in range(sum):
result += (A_list[i]*f(x_k_list[i])).evalf()
print(result)
结果:
2.39970807094290
目标函数求积区间不是[-1,1]时
代码:
#高斯-勒让德求积公式
from sympy import *
from scipy.special import perm,comb #排列,组合
x,t = symbols("x,t")
#积分区间
a = 1
b = 2
#需要求积的目标函数
def f(x):
f = 1/x
return f
n = 2 #n次多项式正交,n越大精度越高(n=0,1,2,...)
#勒让德多项式
def L(n):
df = diff((x ** 2 - 1) ** (n + 1), x, n + 1)
# Python内置阶乘函数factorial
L = 1 /2**(n+1)/factorial(n+1) * df
return L
#高斯点x求取
def Gauss_point(n):
x_k_list = solve(L(n)) #求得零点解集
return x_k_list
#求积系数A
def Quadrature_coefficient(x_k_list):
A_list = []
for x_k in x_k_list:
A = 2/(1-x_k**2)/(diff(L(n),x,1).subs(x,x_k))**2
A_list.append(A)
return A_list
result = 0
x_k_list = Gauss_point(n)
A_list = Quadrature_coefficient(Gauss_point(n))
sum = len(A_list)
#区间变换
if a == -1 and b == 1:
for i in range(sum):
result += (A_list[i] * f(x_k_list[i])).evalf()
print(result)
#将求求粉公式中的区间(a,b)转换为[-1,1]
else:
for i in range(sum):
X = (b-a)/2 * x_k_list[i] + (a+b)/2 #区间变换
result += (b-a)/2 * (A_list[i] * f(X)).evalf()
print(result)
结果:
0.693121693121693