数值积分

应用背景：微积分的基本计算公式为经典的牛顿莱布尼茨公式，但是作为计算机来说，计算原函数等操作过于困难，因此需要其他方法进行优化，代替牛顿莱布尼茨公式进行计算，本文将分享两种求积公式及其代码实现。

使用语言： Python-Numpy

使用两种方法计算$f(x)=sin(x)$在[1,2]上的积分

复化梯形的递推公式

由梯形公式改进，将区间[a, b]等分为n个小区间 $[x_i, x_{i+1}]$

$\int_{a}^{b}f(x)dx=\sum_{i=0}^{n-1}\int_{x_i}^{x_{i+1}}f(x)dx\approx \frac{h_i}{2}[f(x_i)+f(x_{i+1})]$

根据梯形公式化简为区间端点值组合

将求和公式展开，分离出左右端点

得到最终的化简结果

$T_n = \frac{h}{2} [f(a)+2\sum_{i=1}^{n-1}f(x_i)+f(b)]$

但是该方法的步长太小的话会导致计算次数太大，步长太大又难以保证效率

采用区间不断对分的方法，取$n = 2^k$，反复使用复合求积公式

$T^{(k)}=\frac{1}{2}T^{(k-1)}+\frac{h_{k-1}}{2}\sum_{i=0}^{2^{k-1}-1}f(a+ih_{k-1} + 0.5h_{k-1})$

其中$h_{k-1} = \frac{b-a}{2^{k-1}}$

Code

 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


import numpy as np

def f(x):
    return np.sin(x)

def Calc():
    eps = 1e-7
    Flag = True
    l = 1.0
    r = 2.0
    h = r - l
    t1 = 1.0*(h/2)*(1+f(r))
    t2 = 0.0
    while Flag:
        sum = 0
        x = l + h/2
        while x < r:
            sum += f(x)
            x += h
        t2 = t1/2 + h*sum/2
        h /= 2.0
        if abs(t2 - t1) < abs(eps):
            Flag = False
        t1 = t2
    res = t2
    return res

answer = Calc()
print(answer, 6)

0.9564492180077249

Romberg 算法

基于梯形递推公式，将上一级递推公式的结果进行线性组合，得到Simpson公式，再对Simpson公式结果进行线性组合得出Cotes公式结果，最后对Cotes公式线性组合得出Romberg算法的递推公式

该方法为Richardson外推法

$T_{m}^{(k)}=\frac{4^m T_{m-1}^{(k+1)} - T_{m-1}^{(k)}}{4^m - 1}$

k为Romberg算法的阶数，Romberg算法具有收敛性，一般情况下我们取四阶 k = 3进行计算就可以满足精度要求

 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


import numpy as np

def f(x):
    return np.sin(x)

# 计算梯形公式的结果
def trapezcomp(l, r, k):
    h = (r - l)/k
    x = l
    Sum = f(x)
    for i in range(1, k):
        x += h
        Sum += 2*f(x)
    return (Sum + f(r))*h*0.5

def romberg(x, y, n):
    '''
    x : 积分下限
    y : 积分上限
    n : Romberg算法求解的阶数
    '''
    I = np.zeros((n, n))
    for i in range(0, n):
        # 获得梯形公式的结果
        I[i, 0] = trapezcomp(x, y, 2**i)
        for j in range(0, i):
            # 使用加速递推公式计算
            I[i, j+1] = (4**(j+1)*I[i, j] - I[i-1, j])/(4**(j+1)-1)
    # 返回4阶Romberg算法的递推值       
    return I[n-1, n-1]

answer = romberg(1.0, 2.0, 4)
print(answer)

运行结果： 0.9564491426149814

数值积分方法

数值积分

复化梯形的递推公式

Romberg 算法