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.
(f, a, b, n)
| 110 | |
| 111 | |
| 112 | def romberg(f, a, b, n): |
| 113 | """Calculate the integral from the Romberg method. |
| 114 | |
| 115 | Args: |
| 116 | f (function): the equation f(x). |
| 117 | a (float): the initial point. |
| 118 | b (float): the final point. |
| 119 | n (int): number of intervals. |
| 120 | |
| 121 | Returns: |
| 122 | xi (float): numerical approximation of the definite integral. |
| 123 | """ |
| 124 | # Initialize the Romberg integration table |
| 125 | r = np.zeros((n, n)) |
| 126 | |
| 127 | # Compute the trapezoid rule for the first column (h = b - a) |
| 128 | h = b - a |
| 129 | r[0, 0] = 0.5 * h * (f(a) + f(b)) |
| 130 | |
| 131 | # Iterate for each level of refinement |
| 132 | for i in range(1, n): |
| 133 | h = 0.5 * h # Halve the step size |
| 134 | # Compute the composite trapezoid rule |
| 135 | sum_f = 0 |
| 136 | for j in range(1, 2**i, 2): |
| 137 | x = a + j * h |
| 138 | sum_f += f(x) |
| 139 | r[i, 0] = 0.5 * r[i - 1, 0] + h * sum_f |
| 140 | |
| 141 | # Richardson extrapolation for higher order approximations |
| 142 | for k in range(1, i + 1): |
| 143 | r[i, k] = r[i, k - 1] + \ |
| 144 | (r[i, k - 1] - r[i - 1, k - 1]) / ((4**k) - 1) |
| 145 | |
| 146 | return float(r[n - 1, n - 1]) |