-
Notifications
You must be signed in to change notification settings - Fork 39
/
integration.py
146 lines (107 loc) · 3.35 KB
/
integration.py
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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
"""Methods for numerical integration."""
import numpy as np
def simpson(f, a, b, n):
"""Calculate the integral from 1/3 Simpson's Rule.
Args:
f (function): the equation f(x).
a (float): the initial point.
b (float): the final point.
n (int): number of intervals.
Returns:
xi (float): numerical approximation of the definite integral.
"""
h = (b - a) / n
sum_odd = 0
sum_even = 0
for i in range(0, n - 1):
x = a + (i + 1) * h
if (i + 1) % 2 == 0:
sum_even += f(x)
else:
sum_odd += f(x)
xi = h / 3 * (f(a) + 2 * sum_even + 4 * sum_odd + f(b))
return xi
def trapezoidal(f, a, b, n):
"""Calculate the integral from the Trapezoidal Rule.
Args:
f (function): the equation f(x).
a (float): the initial point.
b (float): the final point.
n (int): number of intervals.
Returns:
xi (float): numerical approximation of the definite integral.
"""
h = (b - a) / n
sum_x = 0
for i in range(0, n - 1):
x = a + (i + 1) * h
sum_x += f(x)
xi = h / 2 * (f(a) + 2 * sum_x + f(b))
return xi
def simpson_array(x, y):
"""Calculate the integral from 1/3 Simpson's Rule.
Args:
x (numpy.ndarray): x values.
y (numpy.ndarray): y values.
Returns:
xi (float): numerical approximation of the definite integral.
"""
if x.size != y.size:
raise ValueError("'x' and 'y' must have same size.")
h = x[1] - x[0]
n = x.size
sum_odd = 0
sum_even = 0
for i in range(1, n - 1):
if (i + 1) % 2 == 0:
sum_even += y[i]
else:
sum_odd += y[i]
xi = h / 3 * (y[0] + 2 * sum_even + 4 * sum_odd + y[n - 1])
return xi
def trapezoidal_array(x, y):
"""Calculate the integral from the Trapezoidal Rule.
Args:
x (numpy.ndarray): x values.
y (numpy.ndarray): y values.
Returns:
xi (float): numerical approximation of the definite integral.
"""
if x.size != y.size:
raise ValueError("'x' and 'y' must have same size.")
h = x[1] - x[0]
n = x.size
sum_x = 0
for i in range(1, n - 1):
sum_x += y[i]
xi = h / 2 * (y[0] + 2 * sum_x + y[n - 1])
return xi
def romberg(f, a, b, n):
"""Calculate the integral from the Romberg method.
Args:
f (function): the equation f(x).
a (float): the initial point.
b (float): the final point.
n (int): number of intervals.
Returns:
xi (float): numerical approximation of the definite integral.
"""
# Initialize the Romberg integration table
r = np.zeros((n, n))
# Compute the trapezoid rule for the first column (h = b - a)
h = b - a
r[0, 0] = 0.5 * h * (f(a) + f(b))
# Iterate for each level of refinement
for i in range(1, n):
h = 0.5 * h # Halve the step size
# Compute the composite trapezoid rule
sum_f = 0
for j in range(1, 2**i, 2):
x = a + j * h
sum_f += f(x)
r[i, 0] = 0.5 * r[i - 1, 0] + h * sum_f
# Richardson extrapolation for higher order approximations
for k in range(1, i + 1):
r[i, k] = r[i, k - 1] + \
(r[i, k - 1] - r[i - 1, k - 1]) / ((4**k) - 1)
return float(r[n - 1, n - 1])