Python数值计算（27）—— 数值微分

C囧囧

已于 2024-10-21 23:33:09 修改

阅读量585

点赞数 3

分类专栏： Python 文章标签： python numpy 差分微分

于 2024-10-21 23:19:29 首次发布

本文链接：https://blog.youkuaiyun.com/cdinten/article/details/143124070

版权

Python 专栏收录该内容

45 篇文章

订阅专栏

1. 基本知识

微积分的核心是导数，在高数中，通常是利用极限的方式来定义导数的：

当 $\Delta x$ 趋近于零时，就是该点处的导数：

$\frac{\mathrm{d} y}{\mathrm{d} x}=\lim_{\Delta x \to 0}\frac{f(x_i+\Delta x)-f(x_i))}{\Delta x}$

从几何角度看，曲线上点 $x_i$ 处的导数，就是该点处切线的斜率，导数对于函数的变化率至关重要。

2. 具有 $O(h)$ 误差的计算方法

根据前面的定义，虽然我们无法通过数值计算取得极限为零，但是当变化量足够小的时候，是可以通过差分的方法，估算处该点处的导数值的，通过对泰勒多项式的展开，以及差分选用的方式，分为了前向差分，中间差分和后向差分，以下是误差限为 $O(h)$ 的1~4阶差分公式：

#differential method for function

import numpy as np

def diff_h1(f,x,h=1e-2,order=1,method='central'):
    """
    differentiation formula of order 1 to 4 
    using the forward, central, or backward method.
    error at O(h)
    Parameters
    ----------
    f : function
        The function to be differentiated.
    x : float
        The point at which to differentiate the function.
    h : float, optional
        The step size for the finite difference method. The default is 1e-2.
    order : int, optional
        The order of the differentiation. The default is 1.
    method : str, optional
        The method to use for the finite difference method. The default is 'central'.
    
    Returns
    -------
    float
        The derivative of the function at the given point.
    
    Raises
    ------
        ValueError
        If the order is greater than 4
        or 
        if parameter h is less than 1e-3(due to floating point error).
        or 
        if the method is not 'forward', 'central', or 'backward'.
    """ 
    if order==0:
        return f(x)
    if h<1e-3:
        raise ValueError("h must be greater than 1e-3,due to floating point error")
    if order>=5:
        raise ValueError("Order must be less than 5")
    if method not in ['forward','central','backward']:
        raise ValueError("Method must be 'forward', 'central', or 'backward'")
        
    if method == 'forward':
        if order==1:
            return (f(x+h)-f(x))/h
        elif order==2:
            return (f(x+2*h)-2*f(x+h)+f(x))/(h**2)
        elif order==3:
            return (f(x+3*h)-3*f(x+2*h)+3*f(x+h)-f(x))/(h**3)
        elif order==4:
            return (f(x+4*h)-4*f(x+3*h)+6*f(x+2*h)-4*f(x+h)+f(x))/(h**4)
    
    if method == 'central':
        if order==1:
            return (f(x+h)-f(x-h))/(2*h)
        elif order==2:
            return (f(x+h)-2*f(x)+f(x-h))/(h**2)
        elif order==3:
            return (f(x+2*h)-2*f(x+h)+2*f(x-h)-f(x-2*h))/(2*h**3)
        elif order==4:
            return (f(x+2*h)-4*f(x+h)+6*f(x)-4*f(x-h)+f(x-2*h))/(h**4)
        
    if method == 'backward':
        if order==1:
            return (f(x)-f(x-h))/h
        elif order==2:
            return (f(x)-2*f(x-h)+f(x-2*h))/(h**2)
        elif order==3:
            return (f(x)-3*f(x-h)+3*f(x-2*h)-f(x-3*h))/(h**3)
        elif order==4:
            return (f(x)-4*f(x-h)+6*f(x-2*h)-4*f(x-3*h)+f(x-4*h))/h/h/h/h
            
    raise ValueError("Unknown method or order")

注意：上述公式中，计算差分的公式中，h的值不能小于0.0001，因为很有可能此时在高阶（通常是4阶）差分的时候，由于浮点数的精度带来的影响，计算结果会出现较大偏差，这并不是方法本身的问题，而是由于浮点数计算精度导致的。

3. 具有 $O(h^2)$ 误差的计算方法

通过上述各公式，可知在需要计算1阶差分的时候，需要计算两个函数值，计算4阶差分时，就需要用到五个点处的函数值，但是如果通过引入更多的点，实际上是可以构造出误差精度为 $O(h^2)$ 的1~4阶差分公式，如下所示：

def diff_h2(f,x,h=1e-2,order=1,method='central'):
    """
    differentiation formula of order 1 to 4 
    using the forward, central, or backward method.
    error at O(h^2)
    Parameters
    ----------
    f : function
        The function to be differentiated.
    x : float
        The point at which to differentiate the function.
    h : float, optional
        The step size for the finite difference method. The default is 1e-2.
    order : int, optional
        The order of the differentiation. The default is 1.
    method : str, optional
        The method to use for the finite difference method. The default is 'central'.
    
    Returns
    -------
    float
        The derivative of the function at the given point.
    
    Raises
    ------
    ValueError
        If the order is greater than 4 
        or 
        if parameter h is less than 1e-3(due to floating point error).
        or 
        if the method is not 'forward', 'central', or 'backward'.
    
    """ 
    if order==0:
        return f(x)
    if h<1e-3:
        raise ValueError("h must be greater than 1e-3,due to floating point error")
    if order>=5:
        raise ValueError("Order must be less than 5")
    if method not in ['forward','central','backward']:
        raise ValueError("Method must be 'forward', 'central', or 'backward'")
        
    if method == 'forward':
        if order==1:
            return (-f(x+2*h)+4*f(x+h)-3*f(x))/(2*h)
        elif order==2:
            return -(f(x+3*h)+4*f(x+2*h)-5*f(x+h)+2*f(x))/(h**2)
        elif order==3:
            return (-3*f(x+4*h)+14*f(x+3*h)-24*f(x+2*h)+18*f(x+h)-5*f(x))/(2*h**3)
        elif order==4:
            return (-2*f(x+5*h)+11*f(x+4*h)-24*f(x+3*h)+26*f(x+2*h)-14*f(x+h)+3*f(x))/(h**4)
        
    if method == 'central':
        if order==1:
            return (-f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))/(12*h)
        elif order==2:
            return (-f(x+2*h)+16*f(x+h)-30*f(x)+16*f(x-h)-f(x-2*h))/(12*h**2)
        elif order==3:
            return (-f(x+3*h)+8*f(x+2*h)-13*f(x+h)+13*f(x-h)-8*f(x-2*h)+f(x-3*h))/(8*h**3)
        elif order==4:
            return (-f(x+3*h)+12*f(x+2*h)-39*f(x+h)+56*f(x)-39*f(x-h)+12*f(x-2*h)-f(x-3*h))/(6*h**4)
            
    if method == 'backward':
        if order==1:
            return (3*f(x)-4*f(x-h)+f(x-2*h))/(2*h)
        elif order==2:
            return (2*f(x)-5*f(x-h)+4*f(x-2*h)-f(x-3*h))/(h**2)
        elif order==3:
            return (5*f(x)-18*f(x-h)+24*f(x-2*h)-14*f(x-3*h)+3*f(x-4*h))/(2*h**3)
        elif order==4:
            return (3*f(x)-14*f(x-h)+26*f(x-2*h)-24*f(x-3*h)+11*f(x-4*h)-2*f(x-5*h))/h/h/h/h
            
    raise ValueError("Unknown method or order")

4. 测试

测试如下：

假设有原函数：

$f(x)=e^{-x}sinx$

其4阶导函数为：

$f(x)=-4e^{-x}sinx$

在此处需要提出注意，就是在很多时候我们并不一定能够获取到实际的解析函数关系式，但是又需要计算采样点附近的导数值，以上给出一个具体的具有解析形式的函数，主要是为了演示差分逼近导数的算法。

以下测试代码输出真实导数值，以及通过差分方式给出的近似结果如下：

x=0.5
print(df4(x))
h=1e-3
print(diff_h2(f,x,h,4,'central'))
# output:
#-1.1631451528507675
#-1.1632546777680846

可以看到两者已经非常接近，其差为0.0001095249173171，似乎并没有达到号称的 $O(h^2)$ 误差范围，这个问题还是由于分母太小的缘故。

如果选择 $h=0.1$ ：

x=0.5
print(df4(x))
h=1e-1
print(diff_h2(f,x,h,4,'central'))
# output:
#-1.1631451528507675
#-1.1631586495957265

这次两者的差为0.000013496744959，可以看到实际上已经达到了0.00001的误差水平，不仅比前面的误差还小，甚至还比 $O(h^2)$ 误差范围还小很多。

所以，在数值计算的时候，虽然公式可以保证，较小的步长可以产生更加精确的结果，但是由于在实际计算时，由于浮点数计算的限制，以及分母太小，容易造成分子的误差被放大，最终结果是更小的步长并不一定会产生更小的误差。

Python数值计算（27）—— 数值微分

1. 基本知识

2. 具有误差的计算方法

3. 具有误差的计算方法

4. 测试

2. 具有 $O(h)$ 误差的计算方法

3. 具有 $O(h^2)$ 误差的计算方法