From 39cc85200c0d59d24be1a8847b148ae6d6fcbbda Mon Sep 17 00:00:00 2001 From: 123 <629825095@qq.com> Date: Sat, 14 Jun 2025 09:57:46 +0800 Subject: [PATCH] q11 --- 120-5.py | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/120-5.py b/120-5.py index b4753b1..fd1dad8 100644 --- a/120-5.py +++ b/120-5.py @@ -2,12 +2,16 @@ import math def fx(x): if x == 0: - x = 1e-10 # Avoid division by zero - return math.sin(x)/x + x = 1e-10 # Avoid division by zero #如果x能为0,注释掉这行############## + pass + return math.sin(x)/x #把函数改成题干的形式################### # 龙贝格方法 积分 def Romberg(a, b, 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: @@ -15,34 +19,49 @@ def Romberg(a, b, err): 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 - +# 把范围ab与精度换成题干的############ if __name__ == "__main__": a = 0 b = 1 - err = 0.5e-6 + err = 0.5e-6 #把精度要求改成题干的################## result, k = Romberg(a, b, err) print(f"Result: {result}, k={k}")