Add "Romberg method"

Author: cfgnunes

README.md

@@ -71,6 +71,7 @@ python3 main.py
 
 - Composite Trapezoidal method
 - Composite 1/3 Simpson's method
+- Romberg method
 
 ### Initial-value problems for ordinary differential equations
 

integration.py

@@ -1,5 +1,7 @@
 """Methods for numerical integration."""
 
+import numpy as np
+
 
 def simpson(f, a, b, n):
     """Calculate the integral from 1/3 Simpson's Rule.
@@ -11,7 +13,7 @@ def simpson(f, a, b, n):
         n (int): number of intervals.
 
     Returns:
-        xi (float): integral value.
+        xi (float): numerical approximation of the definite integral.
     """
     h = (b - a) / n
 
@@ -39,7 +41,7 @@ def trapezoidal(f, a, b, n):
         n (int): number of intervals.
 
     Returns:
-        xi (float): integral value.
+        xi (float): numerical approximation of the definite integral.
     """
     h = (b - a) / n
 
@@ -61,7 +63,7 @@ def simpson_array(x, y):
         y (numpy.ndarray): y values.
 
     Returns:
-        xi (float): integral value.
+        xi (float): numerical approximation of the definite integral.
     """
     if y.size != y.size:
         raise ValueError("'x' and 'y' must have same size.")
@@ -90,7 +92,7 @@ def trapezoidal_array(x, y):
         y (numpy.ndarray): y values.
 
     Returns:
-        xi (float): integral value.
+        xi (float): numerical approximation of the definite integral.
     """
     if y.size != y.size:
         raise ValueError("'x' and 'y' must have same size.")
@@ -105,3 +107,40 @@ def trapezoidal_array(x, y):
 
     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])]

main.py

@@ -468,6 +468,28 @@ def f(x):
     print(f"xi = {xi:.5f}")
 
 
+@print_docstring
+def example_romberg():
+    """Run an example 'Integration: Romberg method'."""
+    def f(x):
+        return x ** 2 * math.log(x ** 2 + 1)
+
+    a = 0.0
+    b = 2.0
+    h = 0.25
+    n = int((b - a) / h)
+
+    print("Inputs:")
+    print(f"a = {a}")
+    print(f"b = {b}")
+    print(f"n = {n}")
+
+    [xi] = integration.romberg(f, a, b, n)
+
+    print("Output:")
+    print(f"xi = {xi:.5f}")
+
+
 @print_docstring
 def example_ode_euler():
     """Run an example 'ODE: Euler'."""
@@ -739,6 +761,7 @@ def main():
     example_trapezoidal()
     example_simpson_array()
     example_simpson()
+    example_romberg()
 
     # Initial-value problems for ordinary differential equations
     example_ode_euler()