利用梯形公式，simpson公式,Newton公式，Cotes公式,复化梯形公式，复化Simpson公式，复化Cotes公式 计算 定积分 被积函数x/(4+x^2) 积分区间[0,1]

import numpy as np


# 定义被积函数
def f(x):
    return x / (4 + x ** 2)


# 精确值
exact_value = 0.5 * np.log(5 / 4)
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)
    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


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


# 6. 复化辛普森公式
def composite_simpson(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


# 7. 复化四阶科茨公式
def composite_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)
    integral = 0
    for i in range(0, n, 4):
        integral += (2 * delta_x / 45) * (7 * y[i] + 32 * y[i + 1] + 12 * y[i + 2] + 32 * y[i + 3] + 7 * y[i + 4])
    return integral


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

# 简单公式
print("\n简单公式结果:")
trap_result = trapezoidal_rule(a, b, 4)
print(f"梯形公式 (n=4): {trap_result:.8f}, 误差: {abs(trap_result - exact_value):.8f}")

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

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

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

# 复化公式
print("\n复化公式结果:")
comp_trap_result = composite_trapezoidal(a, b, 12)
print(f"复化梯形公式 (n=4): {comp_trap_result:.8f}, 误差: {abs(comp_trap_result - exact_value):.8f}")

comp_simp_result = composite_simpson(a, b, 12)
print(f"复化辛普森公式 (n=4): {comp_simp_result:.8f}, 误差: {abs(comp_simp_result - exact_value):.8f}")

comp_cotes4_result = composite_cotes_4(a, b, 12)
print(f"复化四阶科茨公式 (n=4): {comp_cotes4_result:.8f}, 误差: {abs(comp_cotes4_result - exact_value):.8f}")