Python科学计算（第2版）
上QQ阅读APP看书，第一时间看更新
上一章目录下一章

2.4.7 多项式函数

与本节内容对应的Notebook为：02-numpy/numpy-450-functions-poly.ipynb。

多项式函数是变量的整数次幂与系数的乘积之和，可以用下面的数学公式表示：

由于多项式函数只包含加法和乘法运算，因此它很容易计算，可用于计算其他数学函数的近似值。多项式函数的应用非常广泛，例如在嵌入式系统中经常会用它计算正弦、余弦等函数。在NumPy中，多项式函数的系数可以用一维数组表示，例如可以用下面的数组表示，其中a[0]是最高次的系数，a[-1]是常数项，注意x2的系数为0。

    a = np.array([1.0, 0, -2, 1])

我们可以用poly1d( )将系数转换为poly1d（一元多项式）对象，此对象可以像函数一样调用，它返回多项式函数的值：

    p = np.poly1d(a)
    print type(p)
    p(np.linspace(0, 1, 5))
    <class 'numpy.lib.polynomial.poly1d'>
    array([ 1.        , 0.515625,    0.125 , -0.078125, 0. ])

对poly1d对象进行加减乘除运算相当于对相应的多项式函数进行计算。例如：

    p + [-2, 1] # 和 p + np.poly1d([-2, 1]) 相同
    poly1d([ 1., 0., -4., 2.])

    p *
 p # 两个3次多项式相乘得到一个6次多项式
    poly1d([ 1., 0., -4., 2., 4., -4., 1.])

    p / [1, 1] # 除法返回两个多项式，分别为商式和余式
    (poly1d([ 1., -1., -1.]), poly1d([ 2.]))

由于多项式的除法不一定能正好整除，因此它返回除法所得到的商式和余式。在上面的例子中，商式为<im><??-??>，余式为2。因此将商式和被除式相乘，再加上余式就等于原来的p：

    p == np.poly1d([ 1., -1., -1.]) *
 [1,1] + 2
    True

多项式对象的deriv( )和integ( )方法分别计算多项式函数的微分和积分：

    p.deriv( )
    poly1d([ 3., 0., -2.])

    p.integ( )
    poly1d([ 0.25, 0. , -1. , 1. , 0. ])

    p.integ( ).deriv( ) == p
    True

多项式函数的根可以使用roots( )函数计算：

    r = np.roots(p)
    r
    array([-1.61803399, 1. , 0.61803399])

    p(r) # 将根带入多项式计算，得到的值近似为0
    array([ 2.33146835e-15, 1.33226763e-15, 1.11022302e-16])

而poly( )函数可以将根转换回多项式的系数：

    np.poly(r)
    array([ 1.00000000e+00, -1.66533454e-15, -2.00000000e+00,
            1.00000000e+00])

除了使用多项式对象之外，也可以直接使用NumPy提供的多项式函数对表示多项式系数的数组进行运算。可以在IPython中使用自动补全查看函数名：

    >>> np.poly # 按Tab键
    np.poly np.polyadd np.polydiv np.polyint np.polysub
    np.poly1d np.polyder np.polyfit np.polymul np.polyval

其中的polyfit( )函数可以对一组数据使用多项式函数进行拟合，找到与这组数据的误差平方和最小的多项式的系数。下面的程序用它计算-π/2～π/2区间与sin(x)函数最接近的多项式的系数：

    np.set_printoptions(suppress=True, precision=4)
    
    x = np.linspace(-np.pi / 2, np.pi / 2, 1000) ❶
    y = np.sin(x) ❷
    
    for deg in [3, 5, 7]:
     a = np.polyfit(x, y, deg) ❸
     error = np.abs(np.polyval(a, x) - y) ❹
     print "degree {}: {}".format(deg, a)
     print "max error of order %d:" % deg, np.max(error)
    degree 3: [-0.145 -0.       0.9887 0. ]
    max error of order 3: 0.00894699376707
    degree 5: [ 0.0076 -0.     -0.1658 -0.   0.9998 -0. ]
    max error of order 5: 0.000157408614187
    degree 7: [-0.0002 -0.      0.0083 0.   -0.1667 -0.    1.   0. ]
    max error of order 7: 1.52682558063e-06

❶首先通过linspace( )将-π/2～π/2区间等分为(1000-1)等份。❷计算拟合目标函数sin(x)的值。❸将表示目标函数的数组传递给polyfit( )进行拟合，第三个参数deg为多项式函数的最高阶数。polyfit( )所得到的多项式和目标函数在给定的1000个点之间的误差最小，polyfit( )返回多项式的系数数组。❹使用polyval( )计算多项式函数的值，并计算与目标函数的差的绝对值。

从程序的输出可以看到，由于正弦函数是一个奇函数，因此拟合的多项式系数中偶数次项的系数接近于0。图2-10显示了各阶多项式与正弦函数之间的误差，请注意图中Y轴为对数坐标。

图2-10 各阶多项式近似正弦函数的误差

上一章目录下一章