python 实现simpson approx辛普森算法

simpson approx辛普森算法介绍

辛普森算法（Simpson’s rule）是一种用于数值积分的近似方法，它通过用抛物线来近似函数曲线上的小段来估计定积分的值。这种方法特别适用于那些难以找到精确解析解或需要快速数值解的积分问题。

基本原理

辛普森算法基于以下观察：如果函数𝑓(𝑥)在区间[a,b]上是连续的，并且我们可以找到该区间上三个点（通常是区间的两个端点和一个中点）的函数值，那么我们可以使用这三个点来构造一个二次多项式（抛物线），并用这个抛物线的积分来近似原函数的积分。

公式

对于单个区间[a,b]，辛普森公式为：

$$[{\int }_a^{b}f(x),dx\approx \frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]]$$

这个公式给出了一个基于区间两端和中点函数值的积分近似。

多区间扩展

对于更长的积分区间或需要更高精度的场景，可以将积分区间[a,b]分割成多个小区间，并在每个小区间上应用辛普森公式，然后将所有小区间的结果相加。如果区间被等分为𝑛个小区间（𝑛通常为偶数，以便在每个小区间上都能应用辛普森公式），则整个积分的近似值为：
$$[{\int }_a^{b}f(x),dx\approx \frac{h}{3}\left[f(a)+2\sum _{i=1}^{n/2-1}f(a+(2i-1)h)+4\sum _{i=1}^{n/2}f(a+2ih)+f(b)\right]]$$

其中，$h=\frac{b-a}{n}$是每个小区间的宽度。

自适应辛普森算法

自适应辛普森算法是一种更高级的数值积分方法，它根据当前的误差估计自动调整积分区间的分割，以在保持计算效率的同时提高精度。这种方法通常用于需要高精度解的情况，并且当被积函数在积分区间内变化较大时特别有用。

Python 实现示例

下面是一个简单的Python函数，实现了辛普森算法来计算单个区间的积分近似值：

def simpson(f, a, b, n):
    """
    使用辛普森算法计算定积分的近似值。
    
    参数:
    f -- 被积函数
    a -- 积分下限
    b -- 积分上限
    n -- 区间分割数（应为偶数）
    
    返回:
    积分的近似值
    """
    if n % 2 != 0:
        raise ValueError("n must be even")
    
    h = (b - a) / n
    x = a
    s = f(a) + f(b)
    
    for i in range(1, n, 2):
        s += 4 * f(x + i * h)
        
    for i in range(2, n-1, 2):
        s += 2 * f(x + i * h)
    
    return s * h / 3

# 示例：计算积分 ∫(0, 1) x^2 dx
import math

def f(x):
    return x**2

approx_integral = simpson(f, 0, 1, 10)  # 使用10个小区间
print(f"Approximate integral: {approx_integral}")
print(f"Exact integral: {1/3}")  # 精确值为1/3

请注意，这个实现是基本的辛普森算法，并没有包括自适应部分。对于需要自适应精度的场景，您可能需要实现更复杂的逻辑来动态调整n的值。

simpson approx辛普森算法python实现样例

下面是用Python实现辛普森算法（Simpson’s rule）的代码示例：

def simpsons_rule(f, a, b, n):
    """计算函数 f 在区间 [a, b] 上的定积分近似值，n为偶数个节点数"""
    h = (b - a) / n
    x = a
    integral = f(x)  # 累加第一个节点的函数值
    for i in range(1, n):
        x += h
        if i % 2 == 0:
            integral += 2 * f(x)  # 偶数节点的函数值乘以2
        else:
            integral += 4 * f(x)  # 奇数节点的函数值乘以4
    integral += f(b)  # 累加最后一个节点的函数值
    integral *= h / 3
    return integral

上述代码中，f 是被积函数，a 和 b 是积分区间的左右边界，n 是偶数个节点数。

使用示例：

def f(x):
    return x**2  # 示例函数为 f(x) = x^2

a = 0  # 积分区间左边界
b = 1  # 积分区间右边界
n = 10  # 节点数（偶数）

integral = simpsons_rule(f, a, b, n)
print("积分近似值:", integral)

输出结果：

积分近似值: 0.3333333333333333

这段代码实现了辛普森算法，可以计算给定的函数 f 在给定的区间 [a, b] 上的定积分近似值。算法的基本思想是将积分区间划分为若干个子区间，然后针对每个子区间利用辛普森公式进行积分，最后将所有子区间的积分结果累加得到最终的积分近似值。