Untitled

import numpy as np

def f(x):
    return 1 + np.exp(-x) * np.sin(4 * x)

# Given analyitical solution
exact_integral = (21 * np.exp(1) - 4 * np.cos(4) - np.sin(4)) / (17 * np.exp(1))

# Integral limits
a = 0
b = 1

# Num of intervals
n = 6  # multiple of 2 and 3
h = (b - a) / n


x = np.linspace(a, b, n+1)
y = f(x)

# Trap rule
trapezoidal_integral = h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])

# Simpson's 1/3 rule
simpson_13_integral = h/3 * (y[0] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-2:2]) + y[-1])

# Simpson's 3/8 rule
simpson_38_integral = 3 * h / 8 * (y[0] + 3 * np.sum(y[1:-1:3]) + 3 * np.sum(y[2:-1:3]) + 2 * np.sum(y[3:-2:3]) + y[-1])

# Errors
trapezoidal_error = np.abs(trapezoidal_integral - exact_integral)
simpson_13_error = np.abs(simpson_13_integral - exact_integral) 
simpson_38_error = np.abs(simpson_38_integral - exact_integral) 

# Print results
print(f"Exact Integral: {exact_integral:.16f}")
print(f"Trapezoidal Rule: {trapezoidal_integral:.16f}, Error: {trapezoidal_error:.16f}")
print(f"Simpson's 1/3 Rule: {simpson_13_integral:.16f}, Error: {simpson_13_error:.16f}")
print(f"Simpson's 3/8 Rule: {simpson_38_integral:.16f}, Error: {simpson_38_error:.16f}")