Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import math
def TrapezeStep(function, x_min, x_max, k, old_int):
if k == 1:
return (function(x_min)+function(x_max))*(x_max-x_min)/2
else:
n_steps = 2**(k-2)
h = (x_max - x_min)/n_steps
ii = 0
x = x_min + h/2
for i in range(n_steps):
ii += function(x)
x += h
return old_int/2 + h*ii/2
def TrapezeIntegrator(function, x_min, x_max, precision=1e-9, step_min = 3, step_max = 25, silent = True):
previous_integral = 0
integral = 0
error = precision+1
for step in range(1,step_max):
integral = TrapezeStep(function, x_min, x_max, step, previous_integral)
error = abs(previous_integral - integral)
previous_integral = integral
if not silent:
print("Step {} : precision of {}".format(step, error))
if error < precision and step > step_min:
break
return integral
def RombergIntegrator(function, x_min, x_max, precision=1e-9, step_min = 3, step_max = 20, silent = True):
r = [TrapezeStep(function, x_min, x_max, 1, 0)]
for k in range(2,step_max):
new_r = [ TrapezeStep(function, x_min, x_max, k, r[0]) ]
for m in range(1,k):
new_r.append(1/(4**m - 1) * (4**m * new_r[m-1] - r[m-1]))
error = abs(r[-1] - new_r[-1])
if not silent:
print("Step {} : precision of {}".format(k, error))
r = new_r
if error < precision and k > step_min:
break
return r[-1]
#self.integrals.append(self.r[-1])
if __name__ == "__main__":
def my_func(x):
return math.exp(-0.5*x**2)
ii = RombergIntegrator(my_func, -5, 5 , silent = False)
print (ii)