Python数值计算（8）——Polynomial

C囧囧

已于 2024-09-26 23:09:33 修改

阅读量1.3k

点赞数 10

分类专栏： Python 文章标签： python 开发语言

于 2024-07-27 22:14:29 首次发布

本文链接：https://blog.youkuaiyun.com/cdinten/article/details/140740471

版权

Python 专栏收录该内容

45 篇文章

订阅专栏

最近工作中会用到插值和拟合，所以接下来打算先讨论一下多项式（Polynomial），以及随之相关的几种插值（interpolation）方法。

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

1. Python中多项式概述

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

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

from numpy.polynomial import Polynomial as P

2. 多项式的创建

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

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

另外一种方式是使用其根，通常多项式也可以写成 $f(x)=a_n(x-x_1)...(x-x_{n-1})(x-x_n)$ 的形式，模块也支持这种方式创建，需要使用到fromroots方法：

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

当然，如果最高幂的系数 $a_n$ 不是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. 多项式的值替换

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

$f(x)=3x^3-2x+1\\ f(0)=1\\ f(1)=2\\ f(2)=21\\$

这种计算在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()

效果如下：

4. 多项式的计算

除了进行值替换，多项式之间可以进行代数的基本运算，如加减乘除等，甚至还可以进行微分和积分等常见运算。

分别使用传统的运算法就可以了，唯一需要注意的是除法，应该使用//，而不是/，测试如下：

a=P([1,-2,0,3])
b=P([-2,1])
print(a-b) # 3.0 - 3.0·x + 0.0·x² + 3.0·x³
print(a+b) # -1.0 - 1.0·x + 0.0·x² + 3.0·x³
print(a*b) # -2.0 + 5.0·x - 2.0·x² - 6.0·x³ + 3.0·x⁴
print(a//b) # 10.0 + 6.0·x + 3.0·x²

通常两个多项式的除法不大可能刚好整除，因此，会有余项，计算余项很简单，使用%运算符，例如上例中：

print(a%b) # 21.0

由于b是一个一次式，因此余数为一个常数。

在需要同时获取商和余式的时，最好的做法是使用divmod函数：

quo, rem = divmod(a,b)

甚至还可以使用** 计算多项式的幂：

a=P([1,-2,0,3])
print(a*a)
print(a**2)
# both display: 1.0 - 4.0·x + 4.0·x² + 6.0·x³ - 12.0·x⁴ + 0.0·x⁵ + 9.0·x⁶

积分和微分，使用函数integ和deriv，积分和微分之后仍旧是一个多项式：

a=P([1,-2,0,3])
print(a.integ()) # 0.0 + 1.0·x - 1.0·x² + 0.0·x³ + 0.75·x⁴
print(a.deriv()) # -2.0 + 0.0·x + 9.0·x²