牛顿柯特斯公式及复合形式、龙贝格求积公式，高斯勒让德求积公式

最新推荐文章于 2024-03-30 08:56:43 发布

宏远十一冠王

最新推荐文章于 2024-03-30 08:56:43 发布

阅读量7.4k

点赞数 10

分类专栏： python 文章标签：算法

本文链接：https://blog.csdn.net/zly03/article/details/121738128

版权

python 专栏收录该内容

6 篇文章 0 订阅

订阅专栏

数值积分的研究实现

牛顿柯特斯公式

柯特斯系数

在这里插入图片描述

各阶对应公式

当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

宏远十一冠王

关注

10
点赞
踩
38

收藏

觉得还不错? 一键收藏
0
评论
牛顿柯特斯公式及复合形式、龙贝格求积公式，高斯勒让德求积公式

数值积分的研究实现牛顿柯特斯公式柯特斯系数各阶对应公式当n = 1时，对应的牛顿-柯特斯公式就是是梯形公式当n= 2时，对应的牛顿-柯特斯公式就是辛普森公式当n =4时，对应的牛顿-柯特斯公式就是柯特斯公式柯特斯系数表核心代码实现计算牛顿-柯特斯系数def calu_xishu(self,index,number,array): # 判断(-1)^n-k if (index - number)%2 ==0: res = 1
复制链接

扫一扫