关于使用python完成龙贝格算法
a=0 # 积分下限
b=1 # 积分上限
eps=10**-5 # 精度
T=[] # 复化梯形序列
S=[] # Simpson序列
C=[] # Cotes序列
R=[] # Romberg序列
def f(x): # 被积函数
y=4/(1+x**2)
return y
def Romberg(a,b,eps,f):
h = b - a
T.append(h * (f(a) + f(b)) / 2) # 定义T[0]
ep=eps+1
m=0
while(ep>=eps):
m=m+1
t=0 # 每一次循环都会将t初始化
for i in range(2**(m-1)-1):
t=t+f(a+(2*(i+1)-1)*h/2**m)*h/2**m
t=t+T[-1]/2 # T[-1]指的是数组T中倒数第一个元素
T.append(t)
if m>=1:
S.append((4/3)*T[-1]-(1/3)*T[-2])
if m>=2:
C.append((16/15)*S[-1]-(1/15)*S[-2])
if m>=3:
R.append((64/63)*C[-1]-(1/63)*C[-2])
if m>4:
ep=abs(10*(R[-1]-R[-2]))
Romberg(a,b,eps,f)
print(T)
print(S)
print(C)
print(R)
print("积分结果为:{:.5f}".format(R[-1]))
运行结果
[3.0, 1.5, 1.6911764705882353, 2.1358026537831245, 2.506292558631274, 2.759621739265887, 2.9189057910568406, 3.0145118477143518, 3.07021171660127, 3.1019889340333107, 3.1198359194773158, 3.1297372868164106, 3.1351765796633675, 3.1381404486809648, 3.1397444739929528, 3.140607526927665, 3.1410695722537456, 3.141315854025983, 3.141446624386673, 3.1415158242843315, 3.1415523315867686, 3.1415715389135075, 3.141581619414385, 3.14158689808331, 3.141589656627041, 3.1415910955038555]
[1.0, 1.7549019607843137, 2.2840113815147545, 2.629789193580657, 2.844064799477424, 2.9720004749871585, 3.0463805332668548, 3.0887783395635755, 3.1125813398439908, 3.125784914625317, 3.1330377425961085, 3.1369896772790193, 3.139128405020164, 3.140279149096948, 3.140895211239236, 3.141223587362439, 3.1413979479500616, 3.1414902145069035, 3.141538890916884, 3.141564500687581, 3.141577941355754, 3.141584979581344, 3.141588657639619, 3.141590576141618, 3.14159157512946]
[1.8052287581699347, 2.3192853428967837, 2.652841047718384, 2.858349839870542, 2.9805295200211406, 3.0513392038188347, 3.0916048599833568, 3.1141682065293517, 3.1266651529440725, 3.133521264460828, 3.137253139591213, 3.139270986869574, 3.1403558653687336, 3.1409362820487217, 3.1412454791039854, 3.1414095719892363, 3.141496365610693, 3.141542136010883, 3.1415662080056275, 3.141578837400299, 3.1415854487963832, 3.141588902843504, 3.141590704041751, 3.1415916417286494]
[2.3274449712257814, 2.658135582715552, 2.8616118841904172, 2.982468880023531, 3.0524631670537183, 3.092243997382793, 3.114526354887224, 3.1268635171728776, 3.13363009162776, 3.1373123757043935, 3.1393030161914526, 3.1403730856623713, 3.140945495011896, 3.141250386993751, 3.1414121766382084, 3.1414977432872235, 3.141542862525171, 3.1415665901007817, 3.141579037866881, 3.141585553739178, 3.1415889576696485, 3.1415907326321992, 3.1415916566125683]
积分结果为:3.14159
疑惑:
运行结果应该是
为什么数组T[],S[],C[],R[],中会出现那么多其他结果
题主的问题恐怕不是数组元素多而是计算错误,比如,T数组的第2个元素应该是3.1,题主计算出来却是1.5。下面是我写的一个算法,请测试。
def romberg(f, a, b, eps):
T = [(f(a)+f(b))/2]
for i in range(4):
k, m = pow(2, len(T)), 0
for j in range(1, k, 2):
m += f(j/k)
T.append(T[-1]/2 + m/k)
S = [4*T[i]/3 - T[i-1]/3 for i in range(1, len(T))]
C = [16*S[i]/15 - S[i-1]/15 for i in range(1, len(S))]
R = [64*C[i]/63 - C[i-1]/63 for i in range(1, len(C))]
while abs(R[-1] - R[-2]) > eps:
k, m = pow(2, len(T)), 0
for j in range(1, k, 2):
m += f(j/k)
T.append(T[-1]/2 + m/k)
S.append(4*T[-1]/3 - T[-2]/3)
C.append(16*S[-1]/15 - S[-2]/15)
R.append(64*C[-1]/63 - C[-2]/63)
print('T = ', T)
print('S = ', S)
print('C = ', C)
print('R = ', R)
return R[-1]
romberg(lambda x:4/(1+x*x), 0, 1, 1e-5)
T = [3.0, 3.1, 3.131176470588235, 3.138988494491089, 3.140941612041389]
S = [3.1333333333333337, 3.1415686274509804, 3.141592502458707, 3.1415926512248227]
C = [3.1421176470588232, 3.141594094125889, 3.1415926611425635]
R = [3.141585783761874, 3.1415926383967965]
3.1415926383967965
romberg(lambda x:4/(1+x*x), 0, 1, 1e-10)
T = [3.0, 3.1, 3.131176470588235, 3.138988494491089, 3.140941612041389, 3.1414298931749745, 3.141551963485656]
S = [3.1333333333333337, 3.1415686274509804, 3.141592502458707, 3.1415926512248227, 3.141592653552836, 3.1415926535892167]
C = [3.1421176470588232, 3.141594094125889, 3.1415926611425635, 3.1415926537080368, 3.1415926535916423]
R = [3.141585783761874, 3.1415926383967965, 3.1415926535900285, 3.141592653589795]
3.141592653589795
因为你一直往列表里面append啊