怎么求导来着？别费劲了，试试这个Python的通用求导法

一时失言乱红尘 2023-01-23 13:54 10阅读 0赞

今天的CSDN问答专栏里出现了一个实战中遇到的问题：现有两组数据，一组是时间序列，一组是对应时间序列的旋转角度，请问怎样计算某一时刻的角速度呢？

显然，这是一个求导即可解决的问题。如果将时间记作 t t t，将旋转角度记作 θ \\theta θ，那么 θ \\theta θ就是 t t t的函数：

θ = f ( t ) \\theta=f(t) θ=f(t)

对上式求导，就得到角速度 v v v对时间 t t t的函数：

v = f ′ ( t ) v=f'(t) v=f′(t)

通常，我们把第1个公式叫做原函数，把第2个公式叫做导函数。有了原函数，利用六大基本初等函数的求导公式，或者应用复合函数求导法则，一般都可以得到导函数。

然而，本问题的原函数 f ( t ) f(t) f(t)在哪里呢？我们手头有的，不过是两组离散的测量数据而已。漫说 f ( t ) f(t) f(t)不见踪影，就算 f ( t ) f(t) f(t)已经存在，倘若不是基本初等函数的话，恐怕大多数程序员也早已忘记如何求导了吧？

看起来真的是山穷水尽无从下手了。不过呢，不用担心，今天推荐一个通用的求导方法，可以帮助程序员彻底摆脱求导的困惑。这个通用方法有两个步骤：第一步是通过离散数据找到原函数，第二步是不依赖于求导公式和技巧就可以得到导函数，无论原函数有多么复杂。

# 1. 寻找原函数 #

从两组离散数据中找到它们之间的内在关系，是一个回归问题，而数据插值就是对回归函数的应用。这启发了我从scipy的插值问题中寻找回归函数。为了使方法更具通用性，下面的代码用自变量 x x x代替时间 t t t，用因变量 y y y代替旋转角度 θ \\theta θ，借助scipy模块的一维插值函数interp1d，轻松得到了本文开始给出的两组数据间隐含的函数关系。

>>> import numpy as np
    >>> from scipy import interpolate
    >>> _x = np.linspace(192,199,8)
    >>> _y = np.array([7.085,10.497,14.019,17.683,21.487,25.403,29.402,33.467])
    >>> f = interpolate.interp1d(_x, _y, kind='cubic')

这段代码最后得到的f，就是旋转角度对于时间的函数，自变量的值域范围限定在原始的时间数据的变化范围内。现在我们可以用它来计算192~199范围内任意时刻的旋转角度了。

>>> f(193.3)
    array(11.54026091)
    >>> f(193.35)
    array(11.71517478)
    >>> f(193.351)
    array(11.71867619)

分别输入193.3、193.35、193.351等时刻，简单测试了一下，看上去还挺都象那么回事儿。

# 2. 求解导函数 #

原函数找是找到了，却犹抱琵琶半遮面，究竟什么模样，我们还是不得而知。没关系，下面这个通用的求导方法根本不在意原函数长什么样，只要原函数是真实存在的，就足够了。这个方法的原理，就是令自变量 x x x在某处有一个极小的增量 d x dx dx，将因变量 y y y的增量 d y dy dy与 d x dx dx之比近似视为函数在该处的导数。这里默认 d x = 1.0 × 1 0 − 10 dx=1.0\\times10^\{-10\} dx=1.0×10−10，如有必要，还可以将其设置为更小的数值。

>>> def get_derivative(f, delta=1e-10):
    		"""导函数生成器"""
    		def derivative(x):
    			"""导函数"""
    			return (f(x+delta)-f(x))/delta
    		return derivative

这不是一个闭包函数吗？熟悉闭包的读者一眼就会看出来。没错，导函数生成器就是一个闭包函数，只要告诉它原函数，它就返回对应的导函数。用刚才得到的旋转角度关于时间的原函数来测试一下，轻松得到其导函数，也就是旋转角速度关于时间的函数。

>>> fd = get_derivative(f)

用这个导函数计算一下193.5和195.8两个时刻的角速度：

>>> fd(193.5)
    3.590830143540664
    >>> fd(195.8)
    3.8961776076555026

至此，问题得解。