This commit is contained in:
lwj
2025-06-14 09:57:46 +08:00
parent 639646c5e3
commit 39cc85200c

View File

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