龙贝格函数求积

龙贝格求积公式也称为逐次分半加速法。是数值计算方法之一，用以求解数值积分。是在梯形公式、辛普森公式和柯特斯公式之间关系的基础上，构造出一种加速计算积分的方法。作为一种外推算法，在不增加计算量的前提下提高了误差的精度。

求积步骤

在这里插入图片描述

算法设计

设计思想为

梯形公式经过区间逐步分半的方法的梯形公式求积就等于辛普森公式求积

辛普森公式经过区间逐步分半的方法的梯形公式求积就等于柯特斯公式求积

柯特斯公式经过区间逐步分半的方法的梯形公式求积就等于龙贝格公式求积

逐层递推下去就得到了Romberg公式的一般形式(T数表)

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ozhCVszD-1585540820813)(C:\Users\HUAWEI\AppData\Roaming\Typora\typora-user-images\image-20200330115111105.png)]$

函数接口

T_2n(a, b, n, T_n)：通过区间逐步分半的方法的复化梯形公式求积，对应步骤2

printm™: 打印T数表

fun(x)：修改求积函数

Romberg(a, b, err_min)： #a代表上限，b代表下限,err_min代表最小误差

该函数中采用了两个列表tm,tm1来分别存储相邻两行的值，该函数使用了两重循环，第一重循环计算第一列T0，第二重循环计算每一行的T(m+1), 最后计算误差（对角线上数的差值）与预期误差比较，为防止出现死循环设置了最大迭代次数kmax

代码

import math
import numpy as np 
from numpy import *  def T_2n(a, b, n, T_n):     #计算T0(h/2)的函数, a 积分下限，b 积分上限, n是区间等分数if n<1:                print('n should larger than 1')h = (b - a)/n           #步长      sum_f = 0               #计算T(k,0)求和部分for k in range(0, n):sum_f = sum_f + fun(a + (k + 0.5)*h)T_2n = T_n/2. + sum_f*h/2.return T_2ndef printm(tm):for i in range(len(tm)):if(tm[i]!=0):print("%.5f"%(tm[i]),end=' ')print()def Romberg(a, b, err_min):            #a代表上限，b代表下限,err_min代表最小误差kmax = 6                            #控制循环次数变量tm = zeros(kmax,dtype = float)      # 第m行所有的元素tm1 = zeros(kmax,dtype = float)     #第m+1行所有的元素   tm[0] = 1/2*(b-a)*(fun(a) + fun(b))  # 初始值，梯形求积公式printm(tm)err = 1k = 0np.set_printoptions(precision = 9)while(err>err_min and k <kmax-1):  #控制循环次数和误差n = 2**k                      # n是区间等分数m = 1tm1[0] = T_2n(a, b, n, tm[0]) while(err>err_min and m <= (k+1)):  #控制循环次数tm1[m] = ((4**m)*tm1[m-1]-tm[m-1])/(4**m-1)  #对应步骤三result = tm1[m]err = abs(tm1[m]-tm[m-1])  #求误差，对应步骤四m = m+1tm = np.copy(tm1)k = k+1printm(tm)     return resultdef fun(x):           #被积函数,  x为自变量f = x**4return fprint('\nMétodo de Romberg')
print('----------------------------------')
f1 = Romberg(0, 1,1.e-12)  
print('----------------------------------')
print('result = ',f1)