【数值分析】数值微分与数值积分

发表于 | 分类于 |

本篇作为数值分析的入门教程，介绍一下数值微分与数值积分的常用算法。

概述

数值积分和数值微分是数学中用于近似求解积分和导数的计算方法。

数值积分主要用来求解定积分问题，即当被积函数复杂或无解析解时，通过数值方法近似计算积分值；数值微分则是在无法直接求得函数导数的情况下，通过数值方法近似计算导数。

这两种数值方法在工程、物理和经济学等领域中有着广泛的应用，尤其是在处理复杂系统模型和实验数据时。

本文将对常用的插值算法进行介绍。

数值微分

一阶差分

根据导数定义，可以取函数在某点及附近的函数值作为导数的近似。

最显然的可以在点的两边取：

  • 前向差分法： 前向差分法是通过函数在某点及其附近一点的值来近似该点的导数。公式为： \(f'(x) \approx \frac{f(x+h) - f(x)}{h} = G(x)\) 其中 $ h $ 是一个很小的正数。

  • 后向差分法： 后向差分法与前向差分法类似，但使用的是函数在某点及其前面一点的值。公式为： \(f'(x) \approx \frac{f(x) - f(x-h)}{h} = G(x)\) 其中 $ h $ 是一个很小的正数。

综合上面的方法，提出中心差分法：

  • 中心差分法： 中心差分法使用函数在某点两侧的值来近似该点的导数，通常比前向差分法和后向差分法更精确。公式为： \(f'(x) \approx \frac{f(x+h) - f(x-h)}{2h} = G(x)\) 其中 $ h $ 是一个很小的正数。

    进行误差分析，对$f(x)$在$x+h,x-h$处泰勒展开代入，即可求得误差

    \(\lvert f'(x)- G(X) \rvert \leq Mh^2/6\) 其中 M 为步长h邻域内三阶导数绝对值$|f’’‘(x)|$的最大值。

二阶差分

有了一阶差分法求一阶导数，可以在此基础上继续做二阶差分：

  • 二阶中心差分法： 二阶中心差分法，可以用来近似二阶导数。公式为： \(f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} = G(X)\) 其中 $ h $ 是一个很小的正数。

    进行误差分析，对$f(x)$在$x+h,x-h$处泰勒展开代入，即可求得误差

    \(\lvert f''(x)- G(X) \rvert \leq Mh^2/12\) 其中 M 为步长h邻域内四阶导数绝对值$|f’’’‘(x)|$的最大值。

数值微分的精度受到 $ h $ 的选择和函数的性质影响。

$ h $ 太大，近似误差会增大； $ h $ 太小，可能会导致数值计算中的舍入误差增大。

因此，选择合适的 $ h $ 是数值微分中一个重要的问题。

此外，数值微分还可以通过 Richardson 外推法等技术来提高精度。

数值积分

数值积分是数学中的一种方法，用于近似计算定积分的值。当被积函数的原函数无法用初等函数表示，或者被积函数仅以数据点的形式给出时，数值积分就显得尤为重要。

复合求积公式

根据牛顿莱布尼兹公式，很容易想到的做法是，把函数切分成很多段，在每一段求面积在求和。

  • 矩形法（中点法）：
    • 将积分区间分成若干个小区间，每个小区间上用矩形的面积近似代替曲线下的面积。
    • 公式：\(\int_a^b f(x) \, dx \approx \Delta x \sum_{i=1}^n f(x_i)\)
    • 其中$\Delta x = \frac{b-a}{n}$，$x_i$是第$i$个小区间的中点。
  • 复合梯形法：
    • 将积分区间分成若干个小区间，每个小区间上用梯形的面积近似代替曲线下的面积。
    • 公式：\(\int_a^b f(x) \, dx \approx \frac{\Delta x}{2} \left( f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n) \right)\)
    • 其中$\Delta x = \frac{b-a}{n}$，$x_i = a + i\Delta x$。
  • 复合辛普森法（抛物线法）：
    • 将积分区间分成若干个小区间，每个小区间上用抛物线的面积近似代替曲线下的面积。
    • 公式：\(\int_a^b f(x) \, dx \approx \frac{\Delta x}{3} \left( f(x_0) + 4\sum_{i=1,3,5,\ldots}^{n-1} f(x_i) + 2\sum_{i=2,4,6,\ldots}^{n-2} f(x_i) + f(x_n) \right)\)
    • 其中$\Delta x = \frac{b-a}{n}$，$x_i = a + i\Delta x$，且$n$为偶数。
    • 误差分析依旧根据泰勒展开带入求得，误差为 \(R_n(f) = -\frac{(b-a)}{180} \left(\frac{h}{2}\right)^4 f^{(4)}(\eta), \eta \in (a, b)\)
    • 其中，$ h = \frac{b-a}{n} $ 是每个子区间的宽度，$ f^{(4)}(\eta) $ 是函数在区间 $[a, b]$ 内某点 $\eta$ 处的四阶导数。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def simpson(f, a, b, n):
    """
    复合辛普森法则数值积分
    :param f: 被积函数
    :param a: 积分下限
    :param b: 积分上限
    :param n: 分割区间的数量，n 必须是偶数
    :return: 积分近似值
    """
    if n % 2 == 1:
        raise ValueError("n 必须是偶数")
    
    h = (b - a) / n  # 每个子区间的宽度
    x = [a + i * h for i in range(n + 1)]  # 区间端点
    s = f(x[0]) + f(x[-1])  # 首尾函数值

    # 计算中间函数值
    for i in range(1, n, 2):
        s += 4 * f(x[i])
    for i in range(2, n - 1, 2):
        s += 2 * f(x[i])

    s *= h / 3
    return s

# 示例：计算积分 ∫(x^2)dx 从 0 到 1
def integrand(x):
    return x ** 2

# 使用复合辛普森法则计算积分
result = simpson(integrand, 0, 1, 100)  # n=100，即分成100个子区间
print(f"积分结果：{result}")

龙贝格积分

龙贝格积分（Romberg Integration）是一种高效的数值积分方法，它基于梯形法则和理查森外推技术，通过逐步折半步长，递推得到更高精度的数值积分结果。

  • 基本梯形公式： \(T_0 = \frac{b-a}{2} (f(a) + f(b))\) 这是梯形公式的基本形式，其中 $a$ 和 $b$ 是积分的上下限，$f(a)$ 和 $f(b)$ 是被积函数在端点的值。

  • 递推公式： \(T_{2n} = \frac{1}{2} T_{n} + \frac{h}{2} \sum_{k=0}^{n-1} f\left(x_{k+\frac{1}{2}}\right)\) 可以发现每次 新二分的梯形面积 与 旧梯形面积 和 新加入分点的面积 的递推关系，可以简化计算。

龙贝格算法

通过改写公式，得到龙贝格算法如下：

  • 首先，使用梯形法则计算积分的近似值，记为 \(T_j^0 = [f(a) + f(b)]h/2\)，其中 $ j $ 表示二分次数。
  • 然后，利用公式 \(T_j^k = \frac{4^j}{4^j - 1} T_{j-1}^{k+1} - \frac{1}{4^j - 1} T_{j-1}^k\) 求T表，其中 $ j $ 表示外推次数。
  • 若 $ \lvert T_k^0 - T_{k-1}^0 \rvert <err$，即满足精度，终止计算， $I = T_k^0$ ；否则迭代计算k+1的T表，直到满足精度要求。

自适应积分

自适应积分是一种数值积分方法，它通过动态调整积分区间的划分来提高积分的精度和效率。

这种方法特别适用于被积函数在某些区间内变化剧烈，而在其他区间内变化平缓的情况。

自适应积分的基本思想是将积分区间不断拆分为更小的区间，直到每个小区间内的积分都满足误差要求。这样可以避免在整个积分区间上使用高精度数值积分方法，从而提高了数值积分的效率。

自适应积分在编程上使用了递归的思想： 首先使用某种数值积分方法（如辛普森法则）计算整个区间的积分，然后将区间步长折半，再次使用相同的方法计算区间的积分。

  • 如果两次计算的误差小于给定的误差限，则停止执行，返回计算结果；
  • 如果误差大于给定的误差限，则对两个子区间分别执行相同的步骤。

自适应求积方法的优点在于它可以根据积分函数的不同特点和区间的不同性质来自适应地调整分割点的位置和精度要求，从而使得计算更加高效和精确。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def adaptive_simpson(f, a, b, tol):
    """
    自适应辛普森法则数值积分
    :param f: 被积函数
    :param a: 积分下限
    :param b: 积分上限
    :param tol: 误差容限
    :return: 积分近似值
    """
    def simpson(f, a, b):
        h = (b - a) / 2
        x0, x1, x2 = a, (a + b) / 2, b
        return h / 3 * (f(x0) + 4 * f(x1) + f(x2))

    def recursive_simpson(f, a, b, tol):
        c = (a + b) / 2
        left_integral = simpson(f, a, c)
        right_integral = simpson(f, c, b)
        total_integral = left_integral + right_integral
        if abs(simpson(f, a, b) - total_integral) < 15 * tol:
            return total_integral
        else:
            return recursive_simpson(f, a, c, tol / 2) + recursive_simpson(f, c, b, tol / 2)

    return recursive_simpson(f, a, b, tol)

# 示例：计算积分 ∫(x^2)dx 从 0 到 1
def integrand(x):
    return x ** 2

# 使用自适应辛普森法则计算积分
result = adaptive_simpson(integrand, 0, 1, 1e-6)  # 误差容限为 1e-6
print(f"积分结果：{result}")

高斯积分

高斯积分个人认为是数值积分里最漂亮的方法了，不过需要读者掌握正交多项式和函数逼近等相关概念才方便叙述。 这里为方便读者理解，只对核心思想和用法做简单叙述，不加以推导和证明了。（也许未来会单独写一篇博客讲高斯积分）

高斯积分是一种高效的数值积分方法，它通过在积分区间内选择特定的点（高斯点）和相应的权重来近似积分。

高斯积分法特别适用于低算力高精度要求的情况，因为它可以使用较少的函数评估点来达到较高的精度。

高斯积分法基于正交多项式的性质，选择的高斯点是正交多项式的根。这些点和相应的权重可以使得高斯积分公式对于给定的多项式次数内的所有多项式都是精确的，即对于 $n$ 个高斯点，它可以精确积分 $2n-1$ 次以内的多项式。

高斯积分公式的一般形式为： \(\int_a^b f(x) \rho(x) \, dx \approx \sum_{i=1}^n w_i f(x_i)\) 其中 $\rho(x)$是权函数， $x_i$ 是高斯点，$w_i$ 是相应的权重。

高斯点和权重

高斯点和权重是通过正交多项式确定的。对于不同的正交多项式类型和次数，高斯点和权重也不同。

常用的正交多项式有：

  • 勒让德多项式： $\rho(x) = 1$ 有限区间积分

  • 切比雪夫多项式： $\rho(x) = \frac{1}{\sqrt{1-x^2}}$ 有限区间积分

  • 拉盖尔多项式： $\rho(x) = e^{-x}$ 无穷区间积分

对应的高斯点和权重可以推导，不过更简单的方法是查表获得。