龙贝格求积分算法例题_龙贝格求积公式（Python实现）

最新推荐文章于 2021-12-14 10:34:27 发布

慕容隽

最新推荐文章于 2021-12-14 10:34:27 发布

阅读量2.1k

点赞数

文章标签：龙贝格求积分算法例题

本文链接：https://blog.csdn.net/weixin_29716207/article/details/112431252

版权

本文介绍了龙贝格求积分算法，并通过Python代码展示了如何运用该算法进行数值积分计算，详细解释了算法过程和代码实现。

摘要由CSDN通过智能技术生成

#Author：glm233
#这个程序就是一个黑箱接口，把要进行龙贝格求积的函数放在func函数里，可以自行修改，然后范围就是在Romberg里改

import math

'''
给定一个函数，如：f(x)= x^(3/2），和积分上下限a,b，用机械求积Romberg公式求积分。

'''
import numpy as np


def func(x):
    return x**2+math.sin(x)/x

class Romberg:
    def __init__(self, integ_dowlimit, integ_uplimit):
        '''
        初始化积分上限integ_uplimit和积分下限integ_dowlimit
        输入一个函数，输出函数在积分上下限的积分

        '''
        self.integ_uplimit = integ_uplimit
        self.integ_dowlimit = integ_dowlimit



    def calc(self):
        '''
        计算Richardson外推算法的四个序列

        '''
        t_seq1 = np.zeros(5, 'f')
        s_seq2 = np.zeros(4, 'f')
        c_seq3 = np.zeros(3, 'f')
        r_seq4 = np.zeros(2, 'f')
        # 循环生成hm间距序列
        hm = [(self.integ_uplimit - self.integ_dowlimit) / (2 ** i) for i in range(0,5)]
        print(hm)
        # 循环生成t_seq1
        fa = func(self.integ_dowlimit)
        fb = func(self.integ_uplimit)

        t0 = (1 / 2) * (self.integ_uplimit - self.integ_dowlimit) * (fa+fb)
        t_seq1[0] = t0

        for i in range(1, 5):
            sum = 0
            # 多出来的点的累加和
            for each in range(1, 2**i,2):
                sum =sum + hm[i]*func( self.integ_dowlimit+each * hm[i])#计算两项值
            temp1 = 1 / 2 * t_seq1[i - 1]
            temp2 =sum
            temp =  temp1 + temp2
            # 求t_seql的1-4位
            t_seq1[i] = temp
        print('T序列：'+ str(list(t_seq1)))
        # 循环生成s_seq2
        s_seq2 = [round((4 * t_seq1[i + 1] - t_seq1[i]) / 3,6) for i in range(0, 4)]
        print('S序列：' + str(list(s_seq2)))
        # 循环生成c_seq3
        c_seq3 = [round((4 ** 2 * s_seq2[i + 1] - s_seq2[i]) / (4 ** 2 - 1),6) for i in range(0, 3)]
        print('C序列：' + str(list(c_seq3)))
        # 循环生成r_seq4
        r_seq4 = [round((4 ** 3 * c_seq3[i + 1] - c_seq3[i]) / (4 ** 3 - 1),6) for i in range(0, 2)]
        print('R序列：' + str(list(r_seq4)))
        r_seq5 = [round((4 ** 4 * r_seq4[i + 1] - r_seq4[i]) / (4 ** 4 - 1),6) for i in range(0, 1)]
        print('A序列：' + str(list(r_seq5)))
        return 'end'


rom = Romberg(0.3, 0.8)
print(rom.calc())