数值积分的研究实现

牛顿柯特斯公式

柯特斯系数

各阶对应公式

当n = 1时，对应的牛顿-柯特斯公式就是是梯形公式

当n= 2时，对应的牛顿-柯特斯公式就是辛普森公式

当n =4时，对应的牛顿-柯特斯公式就是柯特斯公式

柯特斯系数表

核心代码实现

计算牛顿-柯特斯系数

def calu_xishu(self,index,number,array):# 判断(-1)^n-kif (index - number)%2 ==0:res = 1else: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 resultdef calu_final(self,array):Jianju = (self.Big-self.Small)/self.driveresult = 0for 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 resultdef calu_final(self,array):result =0   #结果Jianju = (self.Big-self.Small)/self.drive/self.n_drgeetemp = self.Smallfor 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+=Jianjutemp -=Jianju# 乘于hresult *=(self.Big-self.Small)/self.drivereturn result

龙贝格求积公式

由复化梯形求积公式

又梯形推出辛普森的积分值

核心代码和测试数据

计算第一个T^(0)

def calu_first(self):# 算第一个T^0T_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 = Truei = 0while final_judge:if(i ==0):res =self.calu_first()result[i][i] = reselse:for j in range(0,i+1):if(j ==0):next_val = 0next_val = self.getSumMiddle(i)# 递推公式next_val +=1/2*result[i-1][j]result[i][j] = next_valelse:# 计算4^m/4^m-1value_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 = arrayi+=1def getSumMiddle(self,n):time = self.calculate_h(n)JianJu = (self.Big-self.Small)/(time)# sumf(x k+1/2)的值sum =0# 记录已经第几个xk+1count = 0while 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+=1sum = 1/2*JianJu*sumreturn 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])1breakprint(result[len(result)-1])def judgements(self,x_oringal):# 对x不属于[-1,1]进行变换if self.Small ==-1 and self.Big == 1:return x_oringalelse:X_new = (self.Big - self.Small)/2 *(x_oringal)+(self.Small+self.Big)/2return 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))**2All_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