72 lines
2.1 KiB
Python
72 lines
2.1 KiB
Python
import math
|
||
|
||
# 龙贝格方法 积分
|
||
def Romberg(a, b, fx, err):
|
||
table = [[],[],[],[]]
|
||
t00 = (b-a)*(fx(a)+fx(b))/2
|
||
table[0].append(t00)
|
||
print(f"t0(0)={t00}")
|
||
t01 = t10 = t11 = t20 = t21 = t30 = t31 = 0
|
||
k = 1
|
||
while True:
|
||
tmp = 0
|
||
for i in range(1, 2**(k-1)+1):
|
||
tmp += fx(a + (b-a)*(2*i-1)/(2**k))
|
||
t01 = t00/2+(b-a)*tmp/(2**k)
|
||
print(f"t0({k})={t01}")
|
||
table[0].append(t01)
|
||
|
||
if k>1:
|
||
t11 = (4*t01-t00)/3
|
||
print(f"t1({k-1})={t11}")
|
||
table[1].append(t11)
|
||
if k>2:
|
||
t21 = (16*t11-t10)/15
|
||
print(f"t2({k-2})={t21}")
|
||
table[2].append(t21)
|
||
if k>3:
|
||
t31 = (64*t21-t20)/63
|
||
print(f"t3({k-3})={t31}")
|
||
table[3].append(t31)
|
||
if abs(t31-t30) < err:
|
||
break
|
||
t30 = t31
|
||
else:
|
||
t30 = (64*t21-t20)/63
|
||
print(f"t3(0)={t30}")
|
||
table[3].append(t30)
|
||
t20 = t21
|
||
else:
|
||
t20 = (16*t11-t10)/15
|
||
print(f"t2(0)={t20}")
|
||
table[2].append(t20)
|
||
t10 = t11
|
||
else:
|
||
t10 = (4*t01-t00)/3
|
||
print(f"t1(0)={t10}")
|
||
table[1].append(t10)
|
||
t00 = t01
|
||
k += 1
|
||
print("Romberg table:")
|
||
for i in range(len(table)):
|
||
print(f"t{i}: {table[i]}")
|
||
return t31, k
|
||
|
||
|
||
# 积分原函数 ##############################################################
|
||
def fx(x):
|
||
if x == 0:
|
||
x = 1e-10 # Avoid division by zero #如果x能为0，注释掉这行##############
|
||
pass
|
||
return math.sin(x)/x #把函数改成题干的形式###################
|
||
|
||
|
||
if __name__ == "__main__":
|
||
##############################################################################################################
|
||
x_start = 0 # 积分下限
|
||
x_end = 1 # 积分上限
|
||
err = 0.5e-6 # 精度要求 P113
|
||
|
||
result, k = Romberg(x_start, x_end, fx, err)
|
||
print(f"Result: {result}, k={k}")
|