Python数值计算（28）——理查森外推法

1. 基础知识

理查森外推法( Richardson extrapolation)是一种提高某些数值过程精度的简单方法，在数值方法中广泛应用。
理查森外推法的基本思想是通过对原函数进行多次求导，并在每一步求导的基础上进行线性组合，得到一个新的函数，这个新的函数与原函数的差距会逐渐减小，最终趋于零。理查森外推法的优点是简单易行，适用于各种数值计算。
理查森外推法的一般步骤如下：

1. 确定原函数的导数。
2. 确定原函数的导数的导数。
3. 利用导数的导数，构造新的函数。
4. 利用原函数和新的函数的差，计算新的函数的误差。
5. 重复步骤3和4，直到误差达到要求。
6. 利用新的函数求解问题。

其递推公式为：

通过逐步精进，即可得到较高精度的导数值。

 2. 算法实现

# using Richardson extrapolation method
# to calculate the numerical differential of function

import numpy as np

def richardson(f, x, h, n):
    """
    Richardson extrapolation method to approximate the derivative of f at x.
    """
    if n<1:
        raise ValueError("n must be a positive integer.")
    Q=np.zeros((n,n))
    for i in range(n):
        # 1st column,use central difference formula
        p=h/2**i
        Q[i,0]=(f(x+p)-f(x-p))/2/p 
    
    for i in range(1,n):
        for j in range(1,i+1):
            # 2nd to i-th columns
            Q[i,j]=Q[i,j-1]+(Q[i,j-1]-Q[i-1,j-1])/(4**j-1) 
    print(Q)
    return Q[n-1,n-1]

3. 函数测试

假设需要计算函数在处的导数值，显然，因此，测试代码如下：

def f(x):
    return x*np.exp(x)

def df(x):
    return (1 + x)*np.exp(x)
x=2
h=0.2
n=4
ndf=richardson(f, x, h, n)
print(f"x={x} is {ndf} (Numerical derivative)")
print(f"x={x} is {df(x)} (Analytical derivative)")

输出结果为： 

#output
[[22.41416066  0.          0.          0.        ] 
 [22.22878688 22.16699562  0.          0.        ] 
 [22.18256486 22.16715752 22.16716831  0.        ] 
 [22.17101693 22.16716762 22.1671683  22.1671683 ]]
x=2 is 22.16716829679173 (Numerical derivative) 
x=2 is 22.16716829679195 (Analytical derivative)

可见精度已经非常高，达到了小数点后12位。