牛顿柯特斯公式及复合形式、龙贝格求积公式，高斯勒让德求积公式
数值积分的研究实现
牛顿柯特斯公式
柯特斯系数
各阶对应公式
当n = 1时，对应的牛顿-柯特斯公式就是是梯形公式
当n= 2时，对应的牛顿-柯特斯公式就是辛普森公式
当n =4时，对应的牛顿-柯特斯公式就是柯特斯公式
柯特斯系数表
核心代码实现
计算牛顿-柯特斯系数
def calu_xishu(self,index,number,array):
# 判断(-1)^n-k
if (index - number)%2 ==0:
res = 1
else:
res = -1
# 算柯特斯系数被积函数中的(t-j)多项式
for i in range(index+1):
if i!=number:
res =res*(self.x-i)
# 采用sympy库来计算积分
result = integrate(res,(self.x,0,index))
# n! 和(n-k)!
n_fac = self.jieCheng(number)
n_k_fac = self.jieCheng(index-number)
result = result/(index*n_fac*n_k_fac)
array[index-1][number] = result
计算最后的积分值
def calu_Cotes(self):
# 保存柯特斯系数，只算到n=8时，后面出现负值，一般不用
array = np.zeros((8,9),dtype = float)
for i in range(8):
for j in range(i+2):
self.calu_xishu(i+1,j,array)
# 最后的积分值
result = self.calu_final(array)
return result
def calu_final(self,array):
Jianju = (self.Big-self.Small)/self.drive
result = 0
for j in range(self.drive+1):
# 算每一个的积分值，array里面为牛顿科特斯系数
result += array[self.drive-1][j]*(self.function(self.Small+(Jianju*j)))
# 乘于（b-a)
result *=(self.Big-self.Small)
return result
复化牛顿柯特斯公式
牛顿-柯特斯公式在 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Bc27S1C5-1638719876420)(https://www.zhihu.com/equation?tex=n%5Cge+8)] 时近似程度并不好，故不能通过提高阶数的方法来提高求积精度。因此考虑把积分区间划分成固定宽度的子区间，在每个小区间用低阶的牛顿-柯特斯公式，最后再把每个小区间的近似结果加起来。这种方法称为复合求积法。下面讨论复合梯形公式、复合辛普森公式和复合柯特斯公式。
复合辛普森求积公式
牛顿-柯特斯公式在 时近似程度并不好，故不能通过提高阶数的方法来提高求积精度。因此考虑把积分区间划分成固定宽度
的子区间，在每个小区间用低阶的牛顿-柯特斯公式，最后再把每个小区间的近似结果加起来。这种方法称为复合求积法。下面讨论复
合梯形公式、复合辛普森公式和复合柯特斯公式。
核心代码实现
复化牛顿柯特斯系数和上述的牛顿柯特斯公式求积的一样
计算结果
def calu_Cotes(self):
# 计算牛顿柯特斯求积公式的系数
array = np.zeros((8,9),dtype = float)
for i in range(self.n_drgee):
for j in range(i+2):
self.calu_xishu(i+1,j,array)
# 最后的积分值
result = self.calu_final(array)
return result
def calu_final(self,array):
result =0 #结果
Jianju = (self.Big-self.Small)/self.drive/self.n_drgee
temp = self.Small
for i in range(self.drive):
for j in range(self.n_drgee+1):
# 算每一个的积分值，array里面为牛顿科特斯系数
f_val = self.function(temp)
result += array[self.n_drgee-1][j]*(f_val)
temp+=Jianju
temp -=Jianju
# 乘于h
result *=(self.Big-self.Small)/self.drive
return result
龙贝格求积公式
由复化梯形求积公式
又梯形推出辛普森的积分值
核心代码和测试数据
计算第一个T^(0)
def calu_first(self):
# 算第一个T^0
T_0 = (self.Big - self.Small)/2*(self.function(self.Small)+self.function(self.Big))
return T_0
递推部分
def final_calu(self):
result = np.zeros((10,10),dtype = float)
# print(result)
final_judge = True
i = 0
while final_judge:
if(i ==0):
res =self.calu_first()
result[i][i] = res
else:
for j in range(0,i+1):
if(j ==0):
next_val = 0
next_val = self.getSumMiddle(i)
# 递推公式
next_val +=1/2*result[i-1][j]
result[i][j] = next_val
else:
# 计算4^m/4^m-1
value_fenmu = math.pow(4,j)-1
# 递推公式
result[i][j] = (math.pow(4,j)/value_fenmu)*result[i][j-1] - (1/value_fenmu)*result[i-1][j-1]
# 判断是否已经达到精确值
if i>=1 and math.fabs(result[i][i] - result[i-1][i-1]) < self.accurancy:
print(f"最后的结果是：{result[i][i]}")
final_judge =False
# 判断是否需要加大矩阵
if (i+1)%len(result) == 0:
print("哈哈")
array = np.zeros((2*(i+1),2*(i+1)),dtype = float)
for j in range(len(result)):
for k in range(len(result[i])):
array[j][k] = result[j][k]
result = array
i+=1
def getSumMiddle(self,n):
time = self.calculate_h(n)
JianJu = (self.Big-self.Small)/(time)
# sumf(x k+1/2)的值
sum =0
# 记录已经第几个xk+1
count = 0
while count<time:
# 1/2;1/4,3/4;1/8,3/8,5/8,7/8;
scale = ((2*count)+1)/(time*2)
sum+=self.function(self.Small+scale*(self.Big-self.Small))
count+=1
sum = 1/2*JianJu*sum
return sum
高斯公式
只计算权函数p(x) = 1的情况
高斯勒让德求积公式
高斯勒让德求积公式核心代码
计算高斯勒让德公式的零点
def calu_intergration(self):
# 计算积分
result = [] #存放上一个积分计算的结果
for i in range(10):
# 进行解方程
biaoDaShi = solve(self.calu_xishu(i)) #求得零点解集
res = self.calu_n_intergration(i,biaoDaShi)
result.append(res)
if(i>=1 and i<10 and math.fabs(result[i]-result[i-1])<self.accurancy):
# print(result[i])1
break
print(result[len(result)-1])
def judgements(self,x_oringal):
# 对x不属于[-1,1]进行变换
if self.Small ==-1 and self.Big == 1:
return x_oringal
else:
X_new = (self.Big - self.Small)/2 *(x_oringal)+(self.Small+self.Big)/2
return X_new
得到AK系数
def calu_n_intergration(self,n,biaoDaShi):
# 算n=1,2,3,4,5,6,7,8,9的对应积分值
All_Ak = []
for val in biaoDaShi:
# subs()通过替换x ,令x =系数
A = 2/(1-val**2)/(diff(self.calu_xishu(n),self.x,1).subs(self.x,val))**2
All_Ak.append(A)
length = len(All_Ak)
temp =0
# 对上界和下界不为1，-1的进行处理
for index in range(length):
# 用evalf函数传入变量的值，对表达式进行求值
X_location =self.judgements(biaoDaShi[index])
temp += (self.Big-self.Small)/2*(All_Ak[index]*self.function(X_location)).evalf()
return temp
完整代码在下面链接
链接：https://pan.baidu.com/s/1FPaw_IxdLGcu0sJu4anNBA
提取码：1314
