2.4.8 多项式函数类

numpy.polynomial模块中提供了更丰富的多项式函数类，例如Polynomial、Chebyshev、Legendre等。它们和前面介绍的numpy.poly1d相反，多项式各项的系数按照幂从小到大的顺序排列，下面使用Polynomial类表示多项式x3 -2x＋1，并计算x=2处的值：

    from numpy.polynomial import Polynomial, Chebyshev
    p = Polynomial([1, -2, 0, 1])
    print p(2.0)
    5.0

Polynomial对象提供了众多的方法对多项式进行操作，例如deriv( )计算导函数：

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

切比雪夫多项式是一个正交多项式序列Ti (x)，一个n次多项式可以表示为多个切比雪夫多项式的加权和。在NumPy中，使用Chebyshev类表示由切比雪夫多项式组成的多项式p(x)：

Ti (x)多项式可以通过Chebyshev.basis(i)获得，图2-11显示了0到4次切比雪夫多项式。通过多项式类的convert( )方法可以在不同类型的多项式之间相互转换，转换的目标类型由kind参数指定。例如下面将T4 (x)转换成Polynomial类。由结果可知：。

图2-11 0到4次切比雪夫多项式

    Chebyshev.basis(4).convert(kind=Polynomial)
    Polynomial([ 1., 0., -8., 0., 8.], [-1., 1.], [-1., 1.])

切比雪夫多项式的根被称为切比雪夫节点，可以用于多项式插值。相应的插值多项式能最大限度地降低龙格现象，并且提供多项式在连续函数的最佳一致逼近。下面以函数插值为例演示切比雪夫节点与龙格现象。

❶在[-1,1]区间上等距离取n个取样点。❷使用n阶切比雪夫多项式的根作为取样点。❸使用两种取样点分别对f(x)进行多项式插值，即计算一个多项式经过所有的插值点。图2-12显示了两种插值点所得到的插值多项式，由下左图可知等距离插值多项式在两端有非常大的振荡，这种现象被称为龙格现象，n越大振荡也越大；而下右图采用切比雪夫节点作为插值点，插值多项式的振荡明显减小，并且n越大振荡越小。

图2-12 等距离插值点（左）、切比雪夫插值点（右）

插值与拟合

所谓多项式插值就是找到一个多项式经过所有的插值点。一个n阶多项式有n+1个系数，因此可以通过解方程求解经过n+1个插值点的n阶多项式的系数。fit( )方法虽然计算与目标点拟合的多项式系数，但是当使用n阶多项式拟合n+1的目标点时，多项式将经过所有目标点，因此其结果与多项式插值相同。

    def f(x):
     return 1.0 / ( 1 + 25 *
 x**
2)
    
    n = 11
    x1 = np.linspace(-1, 1, n) ❶
    x2 = Chebyshev.basis(n).roots( ) ❷
    xd = np.linspace(-1, 1, 200)
    
    c1 = Chebyshev.fit(x1, f(x1), n - 1, domain=[-1, 1]) ❸
    c2 = Chebyshev.fit(x2, f(x2), n - 1, domain=[-1, 1])
    
    print u"插值多项式的最大误差：",
    print u"等距离取样点：", abs(c1(xd) - f(xd)).max( ),
    print u"切比雪夫节点：", abs(c2(xd) - f(xd)).max( )
    插值多项式的最大误差：等距离取样点： 1.91556933029 切比雪夫节点： 0.109149825014

在使用多项式逼近函数时，使用切比雪夫多项式进行插值的误差比一般多项式要小许多。在下面的例子中，对g(x)在100个切比雪夫节点之上分别使用Polynomial和Chebyshev进行插值，结果如图2-13所示。在使用Polynomial.fit( )插值时，产生了RankWarning: The fit may be poorly conditioned警告，因此其结果多项式未能经过所有插值点。

    def g(x):
        x = (x - 1) *
 5
        return np.sin(x**
2) + np.sin(x)**
2
    
    n = 100
    x = Chebyshev.basis(n).roots( )
    xd = np.linspace(-1, 1, 1000)
    
    p_g = Polynomial.fit(x, g(x), n - 1, domain=[-1, 1])
    c_g = Chebyshev.fit(x, g(x), n - 1, domain=[-1, 1])
    
    print "Max Polynomial Error:", abs(g(xd) - p_g(xd)).max( )
    print "Max Chebyshev Error:", abs(g(xd) - c_g(xd)).max( )
    Max Polynomial Error: 1.19560558744
    Max Chebyshev Error: 6.47575726376e-09

trim( )方法可以降低多项式的次数，将尾部绝对值小于参数tol的高次系数截断。下面使用trim( )方法获取c中前68个系数，得到一个新的Chebyshev对象c_trimed，其最大误差上升到0.09左右。

    c_trimed = c_g.trim(tol=0.05)
    print "degree:", c_trimed.degree( )
    print "error:", abs(g(xd) - c_trimed(xd)).max( )
    degree: 68
    error: 0.0912094835458

下面用同样的方法对函数h(x)进行19阶的切比雪夫多项式插值，得到插值多项式c_h：

    def h(x):
        x = 5 *
 x
        return np.exp(-x**
2 / 10)
    
    n = 20
    x = Chebyshev.basis(n).roots( )
    c_h = Chebyshev.fit(x, h(x), n - 1, domain=[-1, 1])
    
    print "Max Chebyshev Error:", abs(h(xd) - c_h(xd)).max( )
    Max Chebyshev Error: 1.66544267266e-09

多项式类支持四则运算，下面将c_g和c_h相减得到c_diff，并调用其roots( )计算其所有根。然后找出其中所有的实数根real_roots，它们就是g(x)与h(x)交点的横坐标。图2-14显示了这两条函数曲线以及通过插值多项式计算的交点：

图2-14 使用Chebyshev插值计算两条曲线在[-1，1]之间的所有交点

    c_diff = c_g - c_h
    roots = c_diff.roots( )
    real_roots = roots[roots.imag == 0].real
    print np.allclose(c_diff(real_roots), 0)
    True

切比雪夫多项式在区间[-1,1]上为正交多项式，因此只有在该区间才能对目标函数正确插值。为了对任何区域的目标函数进行插值，需要对自变量的区间进行缩放和平移变换。可以通过domain参数指定拟合点的区间。在下面的例子中，对g2(x)在区间[-10,0]之内使用切比雪夫多项式进行插值。❶为了产生目标区间的切比雪夫节点，在通过basis( )方法创建Tn (x)时，通过domain参数指定目标区间。❷在调用fit( )方法进行拟合时，通过domain参数指定同样的区间。❸最后输出拟合得到的c_g2多项式在[-10,0]中与目标函数的最大误差。

    def g2(x):
        return np.sin(x**
2) + np.sin(x)**
2
    
    n = 100
    x = Chebyshev.basis(n, domain=[-10, 0]).roots( ) ❶
    xd = np.linspace(-10, 0, 1000)
    
    c_g2 = Chebyshev.fit(x, g2(x), n - 1, domain=[-10, 0]) ❷
    
    print "Max Chebyshev Error:", abs(g2(xd) - c_g2(xd)).max( ) ❸
    Max Chebyshev Error: 6.47574571744e-09
上一章目录下一章