17. SciPy求函数的导数

在SciPy里提供了很多的方法函数可以实现对某函数进行求导和求积分的操作。

SciPy的求导相对简单也容易理解。已知函数$f(x)$求其在$x_0$的导数即$f'(x_0)$。

方法一，使用scipy.misc模块下的derivative方法函数。

import numpy as np
from scipy.misc import derivative
def f(x): 
    return x**5
for x in range(1, 4):
    print derivative(f, x, dx=1e-6)

程序的执行结果：

4.999999999866223
80.00000000230045
405.00000005749826

这里的$f(x)= x^5$那么$f'(x) = 5 x^4$，当$x= 1$时$f'(x = 1) = 5$，当$x = 2$时$f'(x = 2) = 5\times 2 ^ 4 = 80$ 。 derivative方法函数还可以指定导数的阶数，下面求一下二阶导数。

import numpy as np
from scipy.misc import derivative
def f(x): 
    return x**5
for x in range(1, 4):
    print derivative(f, x, dx=1e-6, n = 2)

程序的执行结果：

19.999335520992645
160.02132952053216
540.0124791776761

这里的$f(x)= x^5$那么$f''(x) = 20 x^3$，当$x= 1$时$f''(x = 1) = 20$，当$x = 2$时$f''(x = 2) = 20\times 2 ^ 3 = 160$ 。

方法二，使用sympy模块里的diff和symbols函数。diff函数可以得到$f(x)$的导数表达式$f'(x)$,而sympy的symbols函数给出数学表达式里数学符号描述符。

from sympy import diff, symbols
t = symbols('x', real=True)
for i in range(1, 4):
    print diff(t**5, t, i)
    print diff(t**5, t, i).subs(t, i),i

程序的执行结果：

5*x**4
5 1
20*x**3
160 2
60*x**2
540 3

程序里python变量t代表数学公式里的$x$，那么diff(t ** 5, t, i)里的$t ** 5 $的意思是$f(x) = x ^ 5$这样的一个函数，t, 1两个参数的意思是对$x$求i阶导数。

当$i = 1$时，即对$f(x) = x ^ 5$一阶求导，得导数表达式$f'(x) = 5 x^4$，而diff(t**5, t, i).subs(t, i)则是求$f'(x)$当$x = i$的导数的值，即$f'(1) = 5$。

当$i = 2$时，即对$f(x) = x ^ 5$二阶求导，得导数表达式$f''(x) = 20 x^3$，而diff(t**5, t, i).subs(t, i)则是求$f'(x)$当$x = i$的导数的值，即$f'(2) = 160$。

当$i = 3$时，即对$f(x) = x ^ 5$三阶求导，得导数表达式$f^{(3)}(x) = 60 x^2$，而diff(t**5, t, i).subs(t, i)则是求$f'(x)$当$x = 3$的导数的值，即$f'(3) = 540$。

方式三，使用numpy模块里的函数也可实现求函数的导数。

import numpy as np
p = np.poly1d([1,2,0,3,0,5])
print p
for i in range(1, 4):
    print np.polyder(p,i)
    print np.polyder(p,i)(1.0)
for i in range(1, 4):
    print p.deriv(i)
    print p.deriv(i)(1.0)

程序执行结果：

   5     4     2
1 x + 2 x + 3 x + 5
   4     3
5 x + 8 x + 6 x
19.0
    3      2
20 x + 24 x + 6
50.0
    2
60 x + 48 x
108.0
   4     3
5 x + 8 x + 6 x
19.0
    3      2
20 x + 24 x + 6
50.0
    2
60 x + 48 x
108.0

用NumPy的poly1d构造$f(x)$，函数的形参是多项式的系数，最左侧的是最高次数的系数，因此p变量代表的方程为$f(x) = x^5 + 2x^4 + 3 x^2 + 5$。而NumPy的polyder和deriv函数的作用差不多都是对p所代表的多项式求导得到$f(x)$的导数的表达式。下面仅分析一下如下代码：

for i in range(1, 4):
    print np.polyder(p,i)
    print np.polyder(p,i)(1.0)

其结果为：

   4     3
5 x + 8 x + 6 x
19.0
    3      2
20 x + 24 x + 6
50.0
    2
60 x + 48 x
108.0

为了便于阅读，程序加上一些输出提示：

for i in range(1, 4):
    print "f(x)'s %d-th derivative is" %(i)
    print np.polyder(p,i)
    print "f(x)'s %d-th value when i =" % (i), i
    print np.polyder(p,i)(i)

结果为：

f(x)'s 1-th derivative is
   4     3
5 x + 8 x + 6 x
f(x)'s 1-th value when i = 1
19
f(x)'s 2-th derivative is
    3      2
20 x + 24 x + 6
f(x)'s 2-th value when i = 2
262
f(x)'s 3-th derivative is
    2
60 x + 48 x
f(x)'s 3-th value when i = 3
684

方程为$f(x) = x^5 + 2x^4 + 3 x^2 + 5$其一阶导数的表达式为$f'(x) = 5 x^4 + 8 x^3 + 6x$由语句np.polyder(p,i)其$i = 1$时得到，而$f'(1) = 5 + 8 + 6 = 19$由语句np.polyder(p,i)(i)其$i = 1$时得到。

同理当$i = 3$时，语句np.polyder(p,i)计算p所代表的表达式$f(x) = x^5 + 2x^4 + 3 x^2 + 5$,求其三阶导数的表达式为$f^{(3)}(x) = 60 x^2 + 48 x $，而语句np.polyder(p,i)(i)则是计算p所代表的函数的导数在$x = 3$的导数值，即$f^{(3)}(x) = 60 \times 3^2 + 48\times 3 = 540 + 144 = 684$。