利用梯形公式，simpson公式,Newton公式，Cotes公式计算 定积分 被积函数e^x 积分区间[0,1]

import numpy as np

# 定义被积函数
def f(x):
    return np.exp(x)

# 精确值
exact_value = np.exp(1) - np.exp(0)
print(f"精确值: {exact_value:.8f}")

# 1. 梯形公式
def trapezoidal_rule(a, b, n):
    delta_x = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    integral = delta_x / 2 * (y[0] + 2 * np.sum(y[1:-1]) + y[-1])
    return integral

# 2. 辛普森公式
def simpson_rule(a, b, n):
    if n % 2 != 0:
        raise ValueError("n 必须是偶数")
    delta_x = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    integral = delta_x / 3 * (y[0] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-2:2]) + y[-1])
    return integral

# 3. 三阶牛顿-科茨公式 (Simpson 3/8 Rule)
def newton_cotes_3(a, b, n):
    if n % 3 != 0:
        raise ValueError("n 必须是3的倍数")
    delta_x = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    integral = 3 * delta_x / 8 * (y[0] + 3 * np.sum(y[1:-1:3]) + 3 * np.sum(y[2:-1:3]) + 2 * np.sum(y[3:-1:3]) + y[-1])
    return integral

# 4. 四阶牛顿-科茨公式
def cotes_4(a, b, n):
    if n % 4 != 0:
        raise ValueError("n 必须是4的倍数")
    delta_x = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    # 系数为 [7, 32, 12, 32, 7] 的循环
    weights = np.ones(n + 1)
    weights[0::4] = 7
    weights[1::4] = 32
    weights[2::4] = 12
    weights[3::4] = 32
    weights[-1] = 7
    integral = 2 * delta_x / 45 * np.sum(y * weights)
    return integral

# 计算并输出结果
a, b = 0, 1  # 积分区间 [0, 1]


# 梯形公式 (n=2)
trap_result = trapezoidal_rule(a, b, 2)
print(f"梯形公式 (n=2): {trap_result:.8f}, 误差: {abs(trap_result - exact_value):.8f}")

# 辛普森公式 (n=2)
simp_result = simpson_rule(a, b, 2)
print(f"辛普森公式 (n=2): {simp_result:.8f}, 误差: {abs(simp_result - exact_value):.8f}")

# 三阶牛顿-科茨公式 (n=3)
nc3_result = newton_cotes_3(a, b, 3)
print(f"三阶牛顿-科茨公式 (n=3): {nc3_result:.8f}, 误差: {abs(nc3_result - exact_value):.8f}")

# 四阶科茨公式 (n=4)
cotes4_result = cotes_4(a, b, 4)
print(f"四阶科茨公式 (n=4): {cotes4_result:.8f}, 误差: {abs(cotes4_result - exact_value):.8f}")