Romberg(龙贝格)数值积分算法较高效的python实现

Kakaluotuo

已于 2022-10-25 19:20:26 修改

阅读量6.5k

点赞数 5

文章标签： python 数值分析龙贝格算法 romberg

于 2022-04-08 21:25:45 首次发布

本文链接：https://blog.csdn.net/k1052176873/article/details/124051116

版权

1. 原理与公式

2. Python代码实现

# coding=utf-8
import numpy as np
import math


# 二分法梯形公式
# func:需要积分的函数
# x_min: 积分下限
# x_max: 积分上限
# epoch: 二分次数
def compute_Tn(func, x_min=0, x_max=1, epoch=10):
    Tn_list = []
    Tn = 0
    h0 = x_max - x_min  # 积分区间的长度,即初始步长
    h = h0  # 每次迭代计算更新h步长
    x_half_list = np.array([0])  # 二分点列表

    for k in range(epoch + 1):
        if k == 0:  # 初始步长h=h0,取两个端点
            Tn = 0.5 * h * (func(x_min) + func(x_max))  # T1
            Tn_list.append({"T_%d" % 2 ** k: Tn.item(), "k": k})
            x_half_list = np.array([(x_min + x_max) / 2])  # 计算二分点
        else:
            Tn = 0.5 * Tn + 0.5 * h * np.sum(func(x_half_list))  # 上一轮的T2n = 0.5*Tn + 0.5*h*二分点处的函数值之和
            Tn_list.append({"T_%d" % 2 ** k: Tn.item(), "k": k})
            h = 0.5 * h  # 更新步长h为原来的一半
            x_half_list = np.linspace(0, 2 ** k, 2 ** k, endpoint=False)  # 计算下一轮所需的二分点，一共有2^k个点,0,1,2,...2^k-1
            x_half_list = h0 * (2 * x_half_list + 1) / (2 ** (k + 1)) + x_min  # X_(k+1/2)=a + (b-a)*(2n+1)/2^(k+1)
    return Tn_list


# 由梯形公式计算辛普森公式
def compute_Sn(Tn_list):
    Sn_list = []
    for i in range(len(Tn_list) - 1):
        Sn = list(Tn_list[i + 1].values())[0] * 4 / 3 - list(Tn_list[i].values())[0] / 3
        k = list(Tn_list[i + 1].values())[1]
        Sn_list.append({"S_%d" % 2 ** (k - 1): Sn, "k": k})
    return Sn_list


# 由辛普森公式计算柯特斯公式
def compute_Cn(Tn_list=None, Sn_list=None):
    if Sn_list is None:
        Sn_list = compute_Sn(Tn_list)
    Cn_list = []
    for i in range(len(Sn_list) - 1):
        Cn = list(Sn_list[i + 1].values())[0] * 16 / 15 - list(Sn_list[i].values())[0] / 15
        k = list(Sn_list[i + 1].values())[1]
        Cn_list.append({"C_%d" % 2 ** (k - 2): Cn, "k": k})
    return Cn_list


# 由柯特斯公式计算龙贝格公式
def compute_Rn(Tn_list=None, Cn_list=None):
    if Cn_list is None:
        Cn_list = compute_Cn(Tn_list)
    Rn_list = []
    for i in range(len(Cn_list) - 1):
        Rn = list(Cn_list[i + 1].values())[0] * 64 / 63 - list(Cn_list[i].values())[0] / 63
        k = list(Cn_list[i + 1].values())[1]
        Rn_list.append({"R_%d" % 2 ** (k - 3): Rn, "k": k})
    return Rn_list

# 例题测试用函数
def sinx_div_x(x):
    y = np.sin(x) / x
    if not isinstance(y, np.ndarray):
        y = np.array([y])
    y[np.where(np.isnan(y))] = 1
    return y



if __name__ == '__main__':
    Tn_list = compute_Tn(sinx_div_x, x_min=0, x_max=1, epoch=24)
    print(Tn_list)
    Sn_list = compute_Sn(Tn_list)
    print(Sn_list)
    Cn_list = compute_Cn(Sn_list=Sn_list)
    print(Cn_list)
    Rn_list = compute_Rn(Cn_list=Cn_list)
    print(Rn_list)

{'T_1': 0.9207354924039483, 'k': 0}
{'T_2': 0.9397932848061772, 'k': 1}
{'T_4': 0.9445135216653896, 'k': 2}
{'T_8': 0.9456908635827013, 'k': 3}
{'T_16': 0.945985029934386, 'k': 4}
{'T_32': 0.9460585609627681, 'k': 5}
{'T_64': 0.9460769430600631, 'k': 6}
{'T_128': 0.946081538543152, 'k': 7}
{'T_256': 0.946082687411347, 'k': 8}
{'T_512': 0.9460829746282348, 'k': 9}
{'T_1024': 0.9460830464324467, 'k': 10}
{'T_2048': 0.946083064383499, 'k': 11}
{'T_4096': 0.946083068871262, 'k': 12}
{'T_8192': 0.9460830699932028, 'k': 13}
{'T_16384': 0.946083070273688, 'k': 14}
{'T_32768': 0.9460830703438092, 'k': 15}
{'T_65536': 0.9460830703613395, 'k': 16}
{'T_131072': 0.9460830703657221, 'k': 17}
{'T_262144': 0.9460830703668177, 'k': 18}
{'T_524288': 0.9460830703670917, 'k': 19}
{'T_1048576': 0.9460830703671603, 'k': 20}
{'T_2097152': 0.9460830703671774, 'k': 21}
{'T_4194304': 0.9460830703671814, 'k': 22}
{'T_8388608': 0.9460830703671824, 'k': 23}
{'T_16777216': 0.946083070367183, 'k': 24}
{'S_1': 0.9461458822735869, 'k': 1}
{'S_2': 0.9460869339517937, 'k': 2}
{'S_4': 0.9460833108884718, 'k': 3}
{'S_8': 0.9460830853849478, 'k': 4}
{'S_16': 0.9460830713055621, 'k': 5}
{'S_32': 0.9460830704258281, 'k': 6}
{'S_64': 0.9460830703708483, 'k': 7}
{'S_128': 0.9460830703674121, 'k': 8}
{'S_256': 0.9460830703671974, 'k': 9}
{'S_512': 0.9460830703671841, 'k': 10}
{'S_1024': 0.9460830703671832, 'k': 11}
{'S_2048': 0.946083070367183, 'k': 12}
{'S_4096': 0.946083070367183, 'k': 13}
{'S_8192': 0.9460830703671832, 'k': 14}
{'S_16384': 0.946083070367183, 'k': 15}
{'S_32768': 0.946083070367183, 'k': 16}
{'S_65536': 0.946083070367183, 'k': 17}
{'S_131072': 0.946083070367183, 'k': 18}
{'S_262144': 0.946083070367183, 'k': 19}
{'S_524288': 0.9460830703671832, 'k': 20}
{'S_1048576': 0.9460830703671832, 'k': 21}
{'S_2097152': 0.9460830703671828, 'k': 22}
{'S_4194304': 0.9460830703671828, 'k': 23}
{'S_8388608': 0.9460830703671831, 'k': 24}
{'C_1': 0.9460830040636741, 'k': 2}
{'C_2': 0.946083069350917, 'k': 3}
{'C_4': 0.9460830703513795, 'k': 4}
{'C_8': 0.9460830703669365, 'k': 5}
{'C_16': 0.9460830703671791, 'k': 6}
{'C_32': 0.946083070367183, 'k': 7}
{'C_64': 0.946083070367183, 'k': 8}
{'C_128': 0.9460830703671831, 'k': 9}
{'C_256': 0.9460830703671832, 'k': 10}
{'C_512': 0.9460830703671832, 'k': 11}
{'C_1024': 0.946083070367183, 'k': 12}
{'C_2048': 0.9460830703671831, 'k': 13}
{'C_4096': 0.9460830703671833, 'k': 14}
{'C_8192': 0.946083070367183, 'k': 15}
{'C_16384': 0.9460830703671831, 'k': 16}
{'C_32768': 0.9460830703671831, 'k': 17}
{'C_65536': 0.9460830703671831, 'k': 18}
{'C_131072': 0.9460830703671831, 'k': 19}
{'C_262144': 0.9460830703671833, 'k': 20}
{'C_524288': 0.9460830703671832, 'k': 21}
{'C_1048576': 0.9460830703671828, 'k': 22}
{'C_2097152': 0.9460830703671829, 'k': 23}
{'C_4194304': 0.9460830703671831, 'k': 24}
{'R_1': 0.9460830703872224, 'k': 3}
{'R_2': 0.9460830703672599, 'k': 4}
{'R_4': 0.9460830703671834, 'k': 5}
{'R_8': 0.946083070367183, 'k': 6}
{'R_16': 0.946083070367183, 'k': 7}
{'R_32': 0.946083070367183, 'k': 8}
{'R_64': 0.9460830703671831, 'k': 9}
{'R_128': 0.9460830703671832, 'k': 10}
{'R_256': 0.9460830703671832, 'k': 11}
{'R_512': 0.946083070367183, 'k': 12}
{'R_1024': 0.9460830703671831, 'k': 13}
{'R_2048': 0.9460830703671833, 'k': 14}
{'R_4096': 0.946083070367183, 'k': 15}
{'R_8192': 0.9460830703671831, 'k': 16}
{'R_16384': 0.9460830703671831, 'k': 17}
{'R_32768': 0.9460830703671831, 'k': 18}
{'R_65536': 0.9460830703671831, 'k': 19}
{'R_131072': 0.9460830703671833, 'k': 20}
{'R_262144': 0.9460830703671832, 'k': 21}
{'R_524288': 0.9460830703671828, 'k': 22}
{'R_1048576': 0.9460830703671829, 'k': 23}
{'R_2097152': 0.9460830703671831, 'k': 24}

Process finished with exit code 0