56 lines
1.1 KiB
Python
56 lines
1.1 KiB
Python
import math
|
|
|
|
def fx(x):
|
|
if x == 0:
|
|
x = 1e-10 # Avoid division by zero
|
|
return math.sin(x)/x
|
|
|
|
|
|
|
|
def Romberg(a, b, err):
|
|
t00 = (b-a)*(fx(a)+fx(b))/2
|
|
t01 = 0
|
|
t10 = 0
|
|
t11 = 0
|
|
t20 = 0
|
|
t21 = 0
|
|
t30 = 0
|
|
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)
|
|
|
|
if k>1:
|
|
t11 = (4*t01-t00)/3
|
|
|
|
if k>2:
|
|
t21 = (16*t11-t10)/15
|
|
|
|
if k>3:
|
|
t31 = (64*t21-t20)/63
|
|
if abs(t31-t30) < err:
|
|
break
|
|
t30 = t31
|
|
else:
|
|
t30 = (64*t21-t20)/63
|
|
t20 = t21
|
|
else:
|
|
t20 = (16*t11-t10)/15
|
|
t10 = t11
|
|
else:
|
|
t10 = (4*t01-t00)/3
|
|
t00 = t01
|
|
k += 1
|
|
return t31, k
|
|
|
|
|
|
if __name__ == "__main__":
|
|
a = 0
|
|
b = 1
|
|
err = 0.5e-6
|
|
result, k = Romberg(a, b, err)
|
|
print(f"Result: {result}, k={k}")
|