Python数值计算（10）—

在scipy中，scipy.interpolate下还有一个PPoly的类，用于表示插值多项式，很多插值算法的结果，都以该类的实例返回，因此有必要了解该类的使用方法。要使用该类，首先要引入相应的模块：

from scipy.interpolate import PPoly

1. 对象创建

PPoly类的构造函数为：

class PPoly(c, x, extrapolate=None, axis=0)

其中c为系数列表，需要注意的是，c是一个k×m的二维数组，c中的每一列是每个分段的系数，x表示间断点横坐标，它必须按升序或者降序方式排列，其行数必须是m+1，x[i]和x[i + 1]之间的多项式以如下形式表示：

S = sum(c[m, i] * (xp - x[i])**(k-m) for m in range(k+1))

其中k 是多项式的度。说了这么多，还是比较抽象，我们直接以一个多项式的创建为例，解释其原理。先创建一个PPoly对象，然后绘制其图形：

import numpy as np
from numpy import poly1d
from numpy.polynomial import Polynomial
from scipy.interpolate import PPoly
import matplotlib.pyplot as plt

c=np.array([[1,-2,3],[0,0,1]])
x=[0,1,2,3]
p1=PPoly(c,x)
X=np.linspace(0,3,1000)
Y=p1(X)
plt.plot(X,Y)
plt.grid()
plt.show()

所生成图形如下：

提供的系数c是一个2×3（k=2，m=3）的数组，说明多项式的度为k=2（因此最高次幂为2-1=1），提供的间断点x，第一个维度大小必须是(m+1)=4。再看看生成的分段多项式，通过简单的计算，可知三段（不含垂直线）分别是：

$\left\{\begin{matrix} y=x,0 \leqslant x<1\\ y=-2x+2,1\leqslant x < 2\\ y=3x-5, 2\leqslant x \leqslant 3 \end{matrix}\right.$

可以看到各段的最高次系数与c的第0列相同，但是第二列就不一样了，原因在哪里呢？PPoly的coefs中的多项式系数是每个区间的本地系数，它的起点是按 $x_i$ 为基础的，上面的函数关系修改为：

$\left\{\begin{matrix} y=1(x-0)+0,0 \leqslant x<1\\ y=-2(x-1)+0,1\leqslant x < 2\\ y=3(x-2)+1, 2\leqslant x \leqslant 3 \end{matrix}\right.$

就可以发现，各段的系数与之保持一致了，且采用了降幂方式排列。

通过上面的过程，验证一下该分段多项式：

c=np.array([[1,-2,3],[0,0,1]])
x=[0,1,2]
p1=PPoly(c.T,x)

注意这里c使用了转置，x也去掉了中断点3，因此k=3,m=2，函数最高次为k-1=2，函数将会分2段，各段的函数应该是：

$y=1(x-0)^2-2\ast(x-0)^1+3=x^2-2x+3,\ 0\le x\le1$

$y=0(x-1)^2+0\ast(x-1)^1+1=1,\ 1\le x\le2$

绘制的图形如下：

当然，如果只有一个函数（不分段），例如常见的多项式， $px=2x^3-3x^2-4x+1$ ，可以表示为：

x=np.array([[2,-3,-4,1]]).T
y=np.array([0,1])
p=PPoly(x,y)

和比较常用的Numpy下的Polynomial对比一下：

from numpy.polynomial import Polynomial
from scipy.interpolate import PPoly
import matplotlib.pyplot as plt

x=np.array([[2,-3,-4,1]]).T
y=np.array([0,1])
p=PPoly(x,y)

X=np.linspace(-1,1,50)
Y=p(X)
p1=Polynomial([1,-4,-3,2])
Y1=p1(X)
plt.plot(X,Y,'r')
plt.plot(X,Y1,'g*')
plt.grid()
plt.show()

运行效果为：

上面的例子先通过现有数据构造PPoly对象p，然后根据p生成点对(X,Y)并上绘制其图形（红色实线），然后构造一个Polynomial对象p2绘制图形（绿色点），可见两者是完全相符的。此外，注意到我们的分段曲线p是在[0,1]上生成的，但是在区间[-1,1]上都是可以正常使用的。

2. 对象属性

通过c和x可以获取其系数以及间断点。例如，通过CubicSpline对数据进行拟合：

x=np.array([0,1,2,3,4])
y=np.array([0,2,-1,1,0])

cs=CubicSpline(x,y)
print(cs.c)
'''
[[ 2.41666667  2.41666667 -2.08333333 -2.08333333]
 [-9.75       -2.5         4.75       -1.5       ]
 [ 9.33333333 -2.91666667 -0.66666667  2.58333333]
 [ 0.          2.         -1.          1.        ]]
'''
print(cs.x) # [0. 1. 2. 3. 4.]

其中x属性很容易理解，就是间断点，c返回的是一个列表，每一列就是每一段的系数，对于分段三次样条曲线而言，分别表示：

$y=2.41666667(x-0)^3-9.75(x-0)^2+9.33333333(x-0)+0,0\le x<1\\ y=2.41666667(x-1)^3-2.5(x-1)^2-2.91666667(x-1)+2,1\le x<2\\ y=-2.08333333(x-2)^3+4.75(x-2)^2-0.66666667(x-2)-1,2\le x<3\\ y=-2.08333333(x-3)^3-1.5(x-3)^2+2.58333333(x-3)+1,3\le x\le4$

3. 数值计算

3.1 单点值

向PPoly对象提供一个标量或者数组，就和函数调用一样，可以计算对应的函数值，以上面的分段三次样条曲线为例：

X=np.linspace(0,4,50)
Y=cs(X)

绘制该曲线如下：

3.2 区间定积分

通过integrate计算PPoly函数在给定区间的积分值，例如，上面的分段多项式在[0.5,1.5]之间的积分：

fint=cs.integrate(0.5,1.5)
print(fint) # 1.791666666666667

3.3 求根

通过solve方法，可以求解方程 $f(x)=y$ 的根，其中y作为solve的参数，依旧以上述分段多项式为例，在区间[0,4]上找到 $y=3$ 时对应的 $x$ 的值，从图形中可知，应该有4个根：

cs=CubicSpline(x,y)
rts=cs.solve(1)
print(rts) # [0.12229237 1.29076033 3. 3.81029911]

测试确实将4个根都找到了。

对于一般的 $f(x)=0$ 的求根，也可以使用roots方法，等效于solve(0)：

cs=CubicSpline(x,y)
rts=cs.roots()
print(rts) # [0. 1.5620559  2.64950957 4.]

3.4 多项式的运算

和Poly1d和Polynomial不同，分段多项式由于其断点的存在，因此不支持一般意义上的多项式四则运算，但是仍旧可以对他进行积分和微分运算。

3.4.1 微分

对PPoly对象进行微分，可以使用PPoly. derivative方法，默认为一阶导数，也可以通过参数指定其他阶导数，返回值仍旧是一个PPoly对象，例如对前面的分段三次样条曲线：

x=np.array([0,1,2,3,4])
y=np.array([0,2,-1,1,0])

cs=CubicSpline(x,y)
X=np.linspace(0,4,50)
Y=cs(X)
csd1=cs.derivative()#一阶导
csd2=cs.derivative(2)#二阶导
plt.plot(X,Y,'r')
plt.plot(X,csd1(X),'g')
plt.plot(X,csd2(X),'b')
plt.grid()
plt.show()

效果如下：

可以看到一阶导数（绿色）是连续可导的，二阶导数（蓝色）是连续的。

3.4.2 积分

对分段多项式进行多项式积分，并不是integrate函数，该函数在前面介绍过，它实际上的功能是数值积分，实际完成多项式积分的是antiderive函数，该函数是原分段函数的不定积分，测试如下：

x=np.array([0,1,2,3,4])
y=np.array([0,2,-1,1,0])

cs=CubicSpline(x,y)
X=np.linspace(0,4,50)
Y=cs(X)
csd1=cs.derivative()#一阶导
csd1i=csd1.antiderivative()
plt.plot(X,Y,'r')
plt.plot(X,csd1i(X),'g*')
plt.grid()
plt.show()

效果如下：

测试代码先计算其导函数csd1，然后又对其进行积分得到csd1i，并将原函数cs和csd1i绘图对比，如果给derive方法提供负整数，其起到的效果和antiderive是一样的。

4. 额外的问题

最后一个问题是，如何通过PPoly对象，构造一个Polynomial对象？

例如，有个这样的分段函数：

c=np.array([[3,1],[0,1],[-2,0],[1,-1]])
x=np.array([0,1,2])
p=PPoly(c,x)

按PPoly的定义，其表达式为：

$y=3(x-0)^3+0(x-0)^2-2(x-0)^1+1=3x^3-2x+1,0\le x<1$

$y=1(x-1)^3+1(x-1)^2+0(x-1)^1-1=x^3-2x^2+x+1,1\le x\le2$

编写转换的函数如下：

def ppolyCoff_to_normal(p:PPoly):
    '''
    由于 PPoly的coefs中的多项式系数是每个区间的本地系数，
    因此必须减去对应节点区间的较低端点，以使用传统多项式方程中的系数。
    即对于区间 [x1,x2] 上的系数 [a,b,c,d]，对应的多项式为
    f(x)=a*(x−x1)**3+b*(x−x1)**2+c*(x−x1)+d .
    此外，PPoly采用的降幂排列，而返回的系数中，
    按Numpy中的Polynomial组织形式升幂排列
    '''
    ppoly_coef=p.c
    x=p.x
    k,m=ppoly_coef.shape
    p_coef=np.zeros((k,m))
    for c in range(m):
        S=P([0])
        for r in range(k):
            S+=ppoly_coef[r,c]*P([-x[c],1])**(k-r-1)
            p_coef[:,c]=S.coef
return p_coef,p.x

测试代码：

c=np.array([[3,1],[0,1],[-2,0],[1,-1]])
x=np.array([0,1,2])
p=PPoly(c,x)
pc,px=ppolyCoff_to_normal(p)
print(pc)
‘’’
[[ 1. -1.]
 [-2.  1.]
 [ 0. -2.]
 [ 3.  1.]]
‘’’

print(px) # [0. 1. 2.]

函数返回系数和间断点，每一列表示每个区间段的多项式，并按升幂排列，绘制上述函数返回的多项式图形：

c=np.array([[3,1],[0,1],[-2,0],[1,-1]])
x=np.array([0,1,2])
p=PPoly(c,x)
pc,px=ppolyCoff_to_normal(p)

X=np.linspace(0,2,100)
Y=p(X)
plt.plot(X,Y,'r')
r,c=pc.shape
for i in range(c):
    p=Polynomial(pc[:,i])
    X=np.linspace(px[i],px[i+1],20)
    Y=p(X)
    plt.plot(X,Y,'b*')
plt.grid()
plt.show()

结果如下图所示：