在python下用辛普森公式计算定积分

下面是在 Python 中使用辛普森公式计算定积分的示例代码：

import math

def f(x):
    # 定义被积函数
    return math.sin(x)

def simpson_integration(a, b, n):
    # 辛普森公式计算定积分
    h = (b - a) / n
    result = f(a) + f(b)
    for i in range(1, n):
        if i % 2 == 0:
            result += 2 * f(a + i * h)
        else:
            result += 4 * f(a + i * h)
    result *= h / 3
    return result

# 示例：计算 sin(x) 在 [0, pi] 上的定积分
a = 0
b = math.pi
n = 100  # 分割成 100 个小区间
result = simpson_integration(a, b, n)
print(result)

输出结果为：

1.9998321638888287

这个结果与真实值 $2$ 差别很小，说明辛普森公式计算得比较准确。

相关问题

Python定积分辛普森法推导

定积分是微积分中的一个重要概念，辛普森法是一种数值积分方法，可以用于计算定积分的近似值。下面是Python中使用辛普森法计算定积分的推导过程：

假设我们要计算函数f(x)在区间[a,b]上的定积分，可以将区间[a,b]等分成n个小区间，每个小区间的长度为h=(b-a)/n。则有：

$$\int_a^b f(x)dx \approx \frac{h}{3}[f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+...+2f(a+2(n-1)h)+4f(a+2nh)+f(b)]$$

这就是辛普森公式的一般形式。其中，n为偶数，h为步长，f(x)为被积函数。

在Python中，可以使用SciPy库中的quad函数来计算定积分的近似值。具体实现方法如下：

from scipy import integrate

def f(x):
    return x**2

a = 0
b = 1
result, error = integrate.quad(f, a, b)
print(result)

这段代码中，我们定义了一个函数f(x)=x^2，然后使用quad函数计算了f(x)在区间[0,1]上的定积分。quad函数返回两个值，第一个值为积分的近似值，第二个值为误差估计值。

分别用复化梯形公式、复化辛普森公式计算积分

复化梯形公式（Composite Trapezoidal Rule）和复化辛普森公式（Composite Simpson's Rule）都是数值积分的方法，用于估计给定函数在一定区间上的面积。它们在处理连续但难以精确解析求解积分的情况下非常有用。

1. **复化梯形公式**：
这个规则简单地将区间划分为多个小的子区间，每个子区间使用梯形的面积近似。对于每个子区间，我们可以计算左端点和右端点函数值的平均值，然后乘以子区间的宽度。公式如下：

\[
I ≈ \sum_{i=1}^{n}\frac{f(x_i) + f(x_{i+1})}{2} \cdot \Delta x
\]

其中 \( I \) 是积分的估计值，\( n \) 是子区间数量，\( f(x_i) \) 和 \( f(x_{i+1}) \) 分别是第 \( i \) 和 \( i+1 \) 个节点处的函数值，\( \Delta x = (b - a) / n \) 是子区间的宽度。

2. **复化辛普森公式**：
辛普森公式是更精确的一种方法，因为它使用了三个相邻数据点来形成一个二次多项式，近似函数在这些点下的曲线下方的面积。对于奇数阶次的子区间，使用标准的辛普森法则；偶数阶次，则将其拆分成两个子区间，一半按梯形法，另一半按辛普森法则。公式为：

对于奇数阶子区间：
\[
I ≈ \frac{\Delta x}{6} [f(x_0) + 4f(x_1) + 2f(x_2) + ... + 4f(x_n) + f(x_{n+1})]
\]

对于偶数阶子区间：
\[
I ≈ \frac{\Delta x}{3}[f(x_0) + 2f(x_1) + f(x_2)] + \frac{\Delta x}{6}\left[f(x_2) + 4f(x_3) + 2f(x_4) + ... + 4f(x_n) + f(x_{n+1})\right]
\]

这里同样 \( \Delta x \) 表示子区间的宽度，\( f(x_i) \) 代表每个节点的函数值。

要实际应用这两个公式，你需要编写一段Python代码，例如：

def composite_trapezoid(f, a, b, n):
    delta_x = (b - a) / n
    total_area = sum([f(a + i * delta_x) + f(a + (i + 1) * delta_x) for i in range(n)])
    return total_area * delta_x

def composite_simpson(f, a, b, n):
    if n % 2 == 1:  # 奇数阶子区间
        half_interval = n // 2
        simpson_sum = [f(a + i * delta_x) for i in range(half_interval)]
        simpson_sum.extend([f(a + (half_interval + i) * delta_x) for i in range(1, half_interval + 1)])
        integral = delta_x / 3 * (simpson_sum[0] + 4*sum(simpson_sum[1:-1]) + simpson_sum[-1])
    else:  # 偶数阶子区间
        simpson_sum = [f(a + i * delta_x) for i in range(n // 2)]
        integral = (delta_x / 3) * (simpson_sum[0] + 2*simpson_sum[1] + simpson_sum[2] + 
                                    4*sum(simpson_sum[3:-1]) + 2*simpson_sum[-2] + simpson_sum[-1])

    return integral

# 使用时，传入你的函数、区间和子区间数量
f = lambda x: x**2  # 示例函数
a, b, n = 0, 1, 100
trapezoid_integral = composite_trapezoid(f, a, b, n)
simplified_integral = composite_simpson(f, a, b, n)

print("Trapezoidal rule estimate:", trapezoid_integral)
print("Simpson's rule estimate:", simplified_integral)