Python数值计算（8）

接下来打算讨论一下多项式（Polynomial），以及随之相关的几种插值（interpolation）方法。

首先介绍一下NumPy中与多项式相关的模块及其使用方式。

1. Python中多项式概述

在数学中，多项式（polynomial）是指由变量、系数以及它们之间的加、减、乘、幂运算（非负整数次方）得到的表达式，例如，等等。，也是多项式，当然后面这两个是包含了多个变量，我们主要讨论前面一种，仅含有一个变量的多项式。

在Numpy中，子模组块polynomial提供了多项式相关的操作，包括了创建及其运算等，将数学上的多项式抽象为Polynomial类，然后每个具体的多项式就是这个类的实例。这个对象的位置很深，完整的类是numpy.polynomial.Polynomial，为何有两个polynomial呢？因为这个包下面还有其他特殊的多项式，例如切比雪夫多项式，勒让德多项式等，这些后面在用到的时候再说。在使用的时候为了方便，在不引起语法歧义的情况下，可以这样：

from numpy.polynomial import Polynomial as P

2. 多项式的创建

一种最简单直观的方式是使用其系数，在该方式下，多项式以其系数的升幂排列，例如，而不是，注意创建的时候缺少的项要用0补齐，该多项式的创建如下：

p=P([1,-2,0,3])
print(p) # 1.0 - 2.0 x + 0.0 x**2 + 3.0 x**3

另外一种方式是使用其根，通常多项式也可以写成的形式，模块也支持这种方式创建，需要使用到fromroots方法：

f=P.fromroots([-1, 1])
print(f) # -1.0 + 0.0 x + 1.0 x**2

当然，如果最高幂的系数不是1，可以整体乘以一个标量：

f=2*P.fromroots([-1, 1])
print(f) # -2.0 + 0.0·x + 2.0·x²

通过上面的两个例子还可以看到，该模块比较“实在”，它严格按照升幂方式排列，而且在某一项的系数为0时也将其打印出来。

另外，在打印多项式的时候，可以看到使用的是Python中的**运算符（Windows环境下是这样）可以修改其显示方式，例如：

a=P([1,-2,0,3])
print(f'{a:unicode}') # 1.0 - 2.0·x + 0.0·x² + 3.0·x³

将会使用unicode方式，将幂用数学上的形式表达。

如果对于所有的多项式，显示都采用Unicode方式，则可以这样：

numpy.polynomial.set_default_printstyle("unicode")

则在此之后默认的显示方式都是使用Unicode形式，本文章后续都采用Unicode方式。

此外，如果也可以看到，默认的多项式变量是x，其实在创建的时候通过指定参数symbol，可以设置改值，并通过实例对象的属性symbol获得该字符串：

a=P([1,-2,0,3],symbol='a')
print(a) # 1.0 - 2.0·a + 0.0·a² + 3.0·a³
print(a.symbol) # a

注意，对于实例对象，symbol在创建的时候如果指定了，就不能再修改，实例对象的symbol属性是只读的。

3. 多项式的属性

除了前面介绍到symbol之外，最直观和最重要的自然是其系数，使用coef属性：

a=P([1,-2,0,3])
print(a.coef) # [ 1. -2.  0.  3.]

另外一个比较意思的是degree函数，获取该多项式的度，这个概念可以理解为是多项式的系数列表的长度减去1，也通常和多项式的最高幂系数相等，但不是绝对的，测试如下：

a=P([1,-2,0,3])
print(a)
print(len(a.coef)-1) # 3
print(a.degree()) # 3

b=P([1,-2,0,3,0])
print(b) # 1.0 - 2.0·x + 0.0·x² + 3.0·x³ + 0.0·x⁴
print(len(b.coef)-1) # 4
print(b.degree()) # 4

从数学角度来说，a和b这两个多项式应该是一样的，但是从代码运行的结果可见，b在输出时，最高次幂虽然系数为0，但是仍旧显示出来了，而且b的度和a的度也不一样。可以将多项式进行裁剪，去掉为0的最高次幂：

print(b.trim().degree()) # 3

另外，还可以计算多项式的根，这需要使用roots函数，还是以上面的多项式为例：

a=P([1,-2,0,3])
print(a.roots()) # [-1. +0.j          0.5-0.28867513j  0.5+0.28867513j]

可以看到，求解是在复数域内完成的，这个也是代数基本定理的：n次复系数多项式方程在复数域内有且只有n个根（m重根算m个），这个定理最先由德国数学家罗特，达朗贝尔、欧拉和拉格朗日都先后尝试证明但不完善，但是在高斯的博士论文中，给出了第一个严格证明（不愧是数学王子），据说目前已经有两百多种证明方法。

4. 多项式的值替换

我们在有了一个多项式之后，对于特定的值，可以通过该多项式计算对应的表达式值，例如函数：

这种计算在Numpy中被称为Evaluation，使用起来也非常简单，和调用函数一样：

a=P([1,-2,0,3])
f0=a(0)
f1=a(1)
f2=a(2)
print([f0,f1,f2]) # [1.0, 2.0, 21.0]

除了对单个点的计算，其实也可以对一个序列进行计算，例如，整体计算并进行绘图：

from numpy.polynomial import Polynomial as P
import numpy as np
import matplotlib.pyplot as plt

np.polynomial.set_default_printstyle("unicode")
f=P([1,-2,0,3])
x=np.linspace(-2,2,50)
y=f(x)
plt.plot(x,y)
plt.plot(1,f(1),'r*') # 显示单个点
plt.grid()
plt.show()

效果如下：