Numerical Methods (Spring 2025)¶

Homework 3¶

Given the following function:

$$ f(x) = |sin(100x)+0.5| $$

  1. Plot the function $f(x)$ in $[0, 0.05]$ and highlight the area under the curve.

  2. Compute $\int_{a}^{b} f(x) dx$, where $a = 0$ and $b = 0.05$, by using Simpson 1/3 and Simpson 3/8 rules (write these two new functions). Use the minimum number of intervals needed by each of the two methods. Discuss their estimation errors. The estimation error is defined as the difference between the integral obtained through these methods and the value returned by scipy.integrate.quad().

  3. Perform the previous computation by splitting the interval into 8 regions. Choose the most appropriate method between Simpson 1/3 and Simpson 3/8 rules and discuss your choice.

  4. Solve the problem discussed in Exercise 2 using Romberg integration with the Trapezoid rule (see notebook section 3 for this function or romb.py). Use max_iter = 3, 4, 5 and discuss the results (explain the function of the max_iter parameter). Compare with the previous exercise.

  5. Compute $\int_{a}^{b} f(x) dx$, where $a = 0$ and $b = 0.05$, by splitting the interval into 13 regions. Choose the most appropriate method amongst the ones you previously used and discuss your choice.

Libraries¶

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import  quad

Ex. 0¶

In [2]:
# Define function
def f(x):
    return np.abs(np.sin(100*x)+0.5)

# Define integration limits
a, b = 0, 0.05

# Generate x values for plotting
x_vals = np.linspace(a, b, 1000)
y_vals = f(x_vals)

# Plot the function
plt.figure(figsize=(8, 5))
plt.plot(x_vals, y_vals, color='b')
#plto vertical lines
plt.axvline(x=a, color='r', linestyle='--', label='a')
plt.axvline(x=b, color='r', linestyle='--', label='b')
# Fill the area under the curve
plt.fill_between(x_vals, y_vals, color='b', alpha=0.1)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Plot of f(x)')
plt.legend()
plt.grid()
plt.show()
No description has been provided for this image

Ex. 1¶

In [3]:
interval = [a, b]
interval
Out[3]:
[0, 0.05]
In [4]:
def compSimp13(f, interval, n):
    
    a = interval[0]  # store the beginning of the interval
    b = interval[1]  # store the end of the interval

    # the % operator returns the remainder 
    if n % 2 == 1 :
        print('Value error: please only use an even number for n')
        return float('nan')
    else:
        h = (b - a) / n  # interval width
        s = f(a) + f(b)  # Get the values of f(a) and f(b)
        
        for i in range(1, n, 2):
            s += 4 * f(a + h * i)

        for i in range(2, n - 1, 2):
            s += 2 * f(a + h * i)

        intF = s * h / 3  # Calculate the integral   
        return intF
In [5]:
def compSimp38(f, interval, n):
    
    a, b = interval
    if n % 3 != 0:
        print('Value error: please use a multiple of 3 for n')
        return float('nan')
    
    h = (b - a) / n
    s = f(a) + f(b)
    
    for i in range(1, n, 3):
        s += 3 * f(a + h * i) + 3 * f(a + h * (i + 1))
    
    for i in range(3, n, 3):
        s += 2 * f(a + h * i)
    
    intF = 3 * h * s / 8
    return intF
In [6]:
simpson_13_result = compSimp13(f, interval, 2) # even number of intervals
simpson_38_result = compSimp38(f, interval, 3) # multiple of 3 intervals
quad_result, quad_error = quad(f, interval[0], interval[1])
In [7]:
print(f'The result of the composite Simpson 1/3 rule is: {simpson_13_result}, with error: {abs(quad_result - simpson_13_result)}')
print(f'The result of the composite Simpson 3/8 rule is: {simpson_38_result}, with error: {abs(quad_result - simpson_38_result)}')
print(f'The result of the quad function is: {quad_result}, with error: {quad_error}')

The result of the composite Simpson 1/3 rule is: 0.044606773758991365, with error: 0.002797729205530286
The result of the composite Simpson 3/8 rule is: 0.03983402662057486, with error: 0.001975017932886222
The result of the quad function is: 0.04180904455346108, with error: 8.359068311025586e-10

The Simpson's 3/8 rule provided a slightly better estimate than the 1/3 rule, with an error approximately 1.5 times smaller.

However, both methods exhibited significantly higher errors compared to the reference scipy.quad function. A clear distinction between the two methods could not be identified, as they share the same theoretical error estimation given by $ O(h^4) $. Consequently, when $ h $ is halved (i.e., the number of intervals is doubled), the estimation error is reduced by a factor of 16 in both cases.

Ex. 2¶

In [8]:
simpson_13_result = compSimp13(f, interval, 8) # even number of intervals
print(f'The result of the composite Simpson 1/3 rule is: {simpson_13_result}, with error: {abs(quad_result - simpson_13_result)}')

The result of the composite Simpson 1/3 rule is: 0.04207197707372662, with error: 0.0002629325202655394

In this case, with ( n = 8 ), Simpson's 1/3 rule should be used since the number of intervals is even.

The error is about 10 times smaller than the one of the previous excercise, with respect to the reference integral. Naturally, the higher number of intervals improved the estimation of the integral.

Ex. 3¶

In [9]:
def romberg(f, a, b, max_iter):
    """Calculate the integral from the Romberg method.

    Args:
        f (function): the equation f(x).
        a (float): the initial point.
        b (float): the final point.
        max_iter (int).

    Returns:
        xi (float): numerical approximation of the definite integral.
    """
    # Initialize the Romberg integration table
    r = np.zeros((max_iter, max_iter))

    # Compute the trapezoid rule for the first column (h = b - a)
    h = b - a
    r[0, 0] = 0.5 * h * (f(a) + f(b))

    for i in range(1, max_iter):
        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 r
In [10]:
max_iter = 3
romb_table = romberg(f, a, b, max_iter)
result = romb_table[max_iter - 1, max_iter - 1]
print(f'The result of the Romberg integration is: {result}, with error: {abs(quad_result - result)}')

The result of the Romberg integration is: 0.03808422537279067, with error: 0.0037248191806704065
In [11]:
max_iter = 4
romb_table = romberg(f, a, b, max_iter)
result = romb_table[max_iter - 1, max_iter - 1]
print(f'The result of the Romberg integration is: {result}, with error: {abs(quad_result - result)}')

The result of the Romberg integration is: 0.04237773600582095, with error: 0.000568691452359868
In [12]:
max_iter = 5
romb_table = romberg(f, a, b, max_iter)
result = romb_table[max_iter - 1, max_iter - 1]
print(f'The result of the Romberg integration is: {result}, with error: {abs(quad_result - result)}')

The result of the Romberg integration is: 0.041878217769367745, with error: 6.917321590666631e-05

The max_iter parameter serves as a stopping criterion in the Romberg function, determining the final size of Romberg’s table. When max_iter = k, the function uses up to $2^k$ intervals. After $k$ extrapolations, the final result has an error of $O(h^{2k})$. To be more precise, for i in range(1, max_iter): is going to stop at max_iter - 1.

This method outperformed the one presented in the previous exercise (Simpson’s rule) for $maxiter \geq 4$. However, when $maxiter = 3$, the error was $0.0037$, which is greater than $0.00026$.

Initially, Romberg integration may perform worse because it relies on successive Richardson extrapolations, which require sufficient refinement to improve accuracy. When $k$ is small, the method has not yet reached a high enough order of accuracy, leading to a larger error compared to Simpson’s rule. However, as $k$ increases, the higher-order convergence of Romberg integration allows it to surpass Simpson’s rule in precision.

Ex. 4¶

In [13]:
interval
Out[13]:
[0, 0.05]
In [14]:
h = (interval[1] - interval[0]) / 13
h
Out[14]:
0.0038461538461538464
In [15]:
simpson_13_result = compSimp13(f, [interval[0], h*10], 10)
print(f'The result of the Simpson 1/3 rule is: {simpson_13_result}')

The result of the Simpson 1/3 rule is: 0.037230624763494945
In [16]:
simpson_38_result = compSimp38(f, [h*10, interval[1]], 3)
print(f'The result of the Simpson 3/8 rule is: {simpson_38_result}')

The result of the Simpson 3/8 rule is: 0.004689313722365598
In [17]:
total1 = simpson_13_result + simpson_38_result
print(f'The total result of the composite Simpson rule is: {total1}, with error: {abs(quad_result - total1)}')

The total result of the composite Simpson rule is: 0.04191993848586054, with error: 0.00011089393239946399
In [ ]:
romberg_res = romberg(f, interval[0], h*8, 4)[-1, -1] # it is goinf to use 2^3 = 8 intervals over the first 8 intervals of the original range
#notice that we need to pass 4 as max_iter to go up to k=3
print(f'The result of the Romberg integration is: {romberg_res}')

The result of the Romberg integration is: 0.03536375805377784
In [19]:
simpson_13_result_2 = compSimp13(f, [h*8, interval[1]], 4)
print(f'The result of the Simpson 1/3 rule is: {simpson_13_result_2}')

The result of the Simpson 1/3 rule is: 0.0062421951988889
In [20]:
total1 = simpson_13_result_2 + romberg_res
print(f'The total result of the Simpson 1/3  + Romberg rule is: {total1}, with error: {abs(quad_result - total1)}')

The total result of the Simpson 1/3  + Romberg rule is: 0.041605953252666736, with error: 0.0002030913007943433

Romberg's method requires the number of intervals to be powers of 2, Simpson 1/3 requires an even number of intervals and Simpson 3/8 requires multiple of 3. In this case, 13 intervals are required, thus a mixed approach shall be used. In general, different options can be used. For example:

  • Simpson 1/3 over the first 10 intervals + Simpson 3/8 on the remaining ones.

  • Romberg over the first 8 intervals ($2^3$) and Simpson 1/3 over the remaining 4 intervals.