- 主观认识
- 样本方程
- 首末分段
- 泰勒级数
- μ mu μ数列与z数列
- Python实现( μ mu μz数列法)
- 代码测试
分段立方是一种多项式插值，是将离散的原始数据还原为一个分段的立方函数（三次函数），英文为Cubic Spline Interpolation。其算法叫立方曲线插值算法Cubic spline interpolation algorithm。其基本思想是生成分段的三次多项式，在连接点处可导，也就是多个分段函数组成的整个函数是连续、可导的。
我举个例子，以下是9个点组成的数据集，组成了9个小向量：
v
0
=
(
1
,
1
)
v
1
=
(
2
,
2
)
v
2
=
(
3
,
1
)
v
3
=
(
4
,
−
2
)
v
4
=
(
5
,
1
)
v
5
=
(
6
,
2
)
v
6
=
(
7
,
1
)
v
7
=
(
8
,
−
2
)
v
8
=
(
9
,
1
)
v_0=(1,1)\ v_1=(2,2)\ v_2=(3,1)\ v_3=(4,-2)\ v_4=(5,1)\ v_5=(6,2)\ v_6=(7,1)\ v_7=(8,-2)\ v_8=(9,1)\
v0=(1,1)v1=(2,2)v2=(3,1)v3=(4,−2)v4=(5,1)v5=(6,2)v6=(7,1)v7=(8,−2)v8=(9,1)
在坐标系上是这个样子：
这个数据，凭直觉看，是一个周期函数，但是计算机不知道啊。如果我们用牛顿插值会是什么样子呢？会得到这样一个函数（近似值，用下面这个方程在绘图软件误差特别大）：
f
(
x
)
=
0.0032
x
8
−
0.130
x
7
+
2.2444
x
6
−
21.11
x
5
+
117.378
x
4
−
390.378
x
3
+
747.3746
x
2
−
740.381
x
+
286.0
f(x)=0.0032x^8-0.130x^7+2.2444x^6\ -21.11x^5+117.378x^4 -390.378x^3\ +747.3746x^2-740.381x+286.0
f(x)=0.0032x8−0.130x7+2.2444x6−21.11x5+117.378x4−390.378x3+747.3746x2−740.381x+286.0
算了，我给出精确方程吧：
f
(
x
)
=
0.0031746031746031746
x
8
−
0.13015873015873017
x
7
+
2.244444444444444
x
6
−
21.11111111111111
x
5
+
117.3777777777778
x
4
−
390.37777777777785
x
3
+
747.3746031746031
x
2
−
740.3809523809525
x
+
286.0
f(x)=0.0031746031746031746x^8\ -0.13015873015873017x^7\ +2.244444444444444x^6\ -21.11111111111111x^5\ +117.3777777777778x^4\ -390.37777777777785x^3\ +747.3746031746031x^2\ -740.3809523809525x+286.0
f(x)=0.0031746031746031746x8−0.13015873015873017x7+2.244444444444444x6−21.11111111111111x5+117.3777777777778x4−390.37777777777785x3+747.3746031746031x2−740.3809523809525x+286.0
在绘图软件(Gnuplot)中得到的图形是这样的：
可见偏离特别大，比如1.5，得到的是
−
3.1679687500005684
-3.1679687500005684
−3.1679687500005684，而我们预期的是一个在
(
1
,
2
)
(1,2)
(1,2)区间内的数字。
我们再试试分段立方插值。分段立方插值得到的是
1.6005154639175259
1.6005154639175259
1.6005154639175259，与真实值比较接近。分段立方插值的图像是什么样子的呢？就是下图这个样子：
与牛顿/拉格朗日插值对比可以更明显看出差异：
主观上看，分段立方插值更接近原始数据，更精确。
分段立方插值获取的是分段立方（三次）函数，假设样本有
n
n
n对数据，那么函数要分成
n
−
1
n-1
n−1，这些样本点就是边界，首先要保证分段函数要经过样本点。
为了简便，我们定义
h
i
=
x
i
+
1
−
x
i
h_i=x_{i+1}-x_{i}
hi=xi+1−xi。
分段函数的形式是
y
=
a
(
x
−
x
i
)
3
+
b
(
x
−
x
i
)
2
+
c
(
x
−
x
i
)
+
d
y=a(x-x_i)^3+b(x-x_i)^2+c(x-x_i)+d
y=a(x−xi)3+b(x−xi)2+c(x−xi)+d，定义域为
[
x
i
,
x
i
+
1
]
[x_i,x_{i+1}]
[xi,xi+1]。要穿过样本点，那么意味着当
x
=
x
i
x=x_i
x=xi时，
y
=
y
i
y=y_i
y=yi，那么我们代进去算一下：
y
i
=
f
(
x
i
)
y
i
=
a
i
(
x
i
−
x
i
)
3
+
b
i
(
x
i
−
x
i
)
2
+
c
i
(
x
i
−
x
i
)
+
d
i
y
i
=
d
i
y_i=f(x_i)\ y_i=a_i(x_i-x_i)^3+b_i(x_i-x_i)^2+c_i(x_i-x_i)+d_i\ y_i=d_i
yi=f(xi)yi=ai(xi−xi)3+bi(xi−xi)2+ci(xi−xi)+diyi=di
所以第一个未知参数
d
d
d确定了，就等于样本的值
y
i
y_i
yi。那么还剩下三个未知数
a
,
b
,
c
a,b,c
a,b,c。我们主观上看，插值函数形成的曲线是丝滑流畅的，那么在数学上，就是连续、可导的。因为每段函数经过两个点，左边界点是
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi)，右边界点是
(
x
i
+
1
,
y
i
+
1
)
(x_{i+1},y_{i+1})
(xi+1,yi+1)。
知道这个条件后，我们得到了第二个方程：
y
i
+
1
=
a
i
(
x
i
+
1
−
x
i
)
3
+
b
i
(
x
i
+
1
−
x
i
)
2
+
c
i
(
x
i
+
1
−
x
i
)
+
d
i
y
i
+
1
=
a
i
(
x
i
+
1
−
x
i
)
3
+
b
i
(
x
i
+
1
−
x
i
)
2
+
c
i
(
x
i
+
1
−
x
i
)
+
y
i
y
i
+
1
=
a
i
h
i
3
+
b
i
h
i
2
+
c
i
h
i
+
y
i
y_{i+1}=a_i(x_{i+1}-x_i)^3+b_i(x_{i+1}-x_i)^2+c_i(x_{i+1}-x_i)+d_i\ y_{i+1}=a_i(x_{i+1}-x_i)^3+b_i(x_{i+1}-x_i)^2+c_i(x_{i+1}-x_i)+y_i\ y_{i+1}=a_ih_i^3+b_ih_i^2+c_ih_i+y_i
yi+1=ai(xi+1−xi)3+bi(xi+1−xi)2+ci(xi+1−xi)+diyi+1=ai(xi+1−xi)3+bi(xi+1−xi)2+ci(xi+1−xi)+yiyi+1=aihi3+bihi2+cihi+yi
这只是连续，在连接处还需要可导。分段函数的形式是一个复合函数所以我们要用链式法则求导：
u
(
x
)
=
x
−
x
i
f
(
u
)
=
a
u
3
+
b
u
2
+
c
u
+
d
u
′
(
x
)
=
1
f
′
(
x
)
=
f
′
(
u
)
⋅
u
′
(
x
)
=
3
a
(
x
−
x
i
)
2
+
2
b
(
x
−
x
i
)
+
c
u(x)=x-x_i\ f(u)=au^3+bu^2+cu+d\ u'(x)=1\ f'(x)=f'(u)cdot u'(x)=3a(x-x_{i})^2+2b(x-x_{i})+c
u(x)=x−xif(u)=au3+bu2+cu+du′(x)=1f′(x)=f′(u)⋅u′(x)=3a(x−xi)2+2b(x−xi)+c
对于中间的分段函数，左连接点要与前一段在连接处的导数相等，也要与后一段在连接处导数相等。假设前一段的分段函数已经解出来了，成为已知数，定义为g(x)，那么我们就得到了第三个方程，起点处的导数等于前一段函数在终点的导数：
f
i
′
(
x
i
)
=
f
i
−
1
′
(
x
i
)
3
a
i
(
x
i
−
x
i
)
2
+
2
b
i
(
x
i
−
x
i
)
+
c
i
=
3
a
i
−
1
(
x
i
−
x
i
−
1
)
2
+
2
b
i
−
1
(
x
i
−
x
i
−
1
)
+
c
i
−
1
c
i
=
3
a
i
−
1
(
x
i
−
x
i
−
1
)
2
+
2
b
i
−
1
(
x
i
−
x
i
−
1
)
+
c
i
−
1
c
i
=
3
a
i
−
1
h
i
−
1
2
+
2
b
i
−
1
h
i
−
1
+
c
i
−
1
f_i'(x_i)=f_{i-1}'(x_i)\ 3a_i(x_i-x_{i})^2+2b_i(x_i-x_{i})+c_i=3a_{i-1}(x_{i}-x_{i-1})^2+2b_{i-1}(x_i-x_{i-1})+c_{i-1}\ c_i=3a_{i-1}(x_{i}-x_{i-1})^2+2b_{i-1}(x_i-x_{i-1})+c_{i-1}\ c_i=3a_{i-1}h_{i-1}^2+2b_{i-1}h_{i-1}+c_{i-1}
fi′(xi)=fi−1′(xi)3ai(xi−xi)2+2bi(xi−xi)+ci=3ai−1(xi−xi−1)2+2bi−1(xi−xi−1)+ci−1ci=3ai−1(xi−xi−1)2+2bi−1(xi−xi−1)+ci−1ci=3ai−1hi−12+2bi−1hi−1+ci−1
此外，在每个分段右边界，也是连续可导的，也就是后一段函数的一阶导相等，又可以得到一个方程：
f
i
′
(
x
i
+
1
)
=
f
i
+
1
′
(
x
i
+
1
)
3
a
i
(
x
i
+
1
−
x
i
)
2
+
2
b
i
(
x
i
+
1
−
x
i
)
+
c
i
=
3
a
i
+
1
(
x
i
+
1
−
x
i
+
1
)
2
+
2
b
i
−
1
(
x
i
+
1
−
x
i
+
1
)
+
c
i
+
1
3
a
i
h
i
2
+
2
b
i
h
i
+
c
i
=
c
i
+
1
f_i'(x_{i+1})=f_{i+1}'(x_{i+1})\ 3a_i(x_{i+1}-x_{i})^2+2b_i(x_{i+1}-x_{i})+c_i=3a_{i+1}(x_{i+1}-x_{i+1})^2+2b_{i-1}(x_{i+1}-x_{i+1})+c_{i+1}\ 3a_ih_{i}^2+2b_ih_{i}+c_i=c_{i+1}
fi′(xi+1)=fi+1′(xi+1)3ai(xi+1−xi)2+2bi(xi+1−xi)+ci=3ai+1(xi+1−xi+1)2+2bi−1(xi+1−xi+1)+ci+13aihi2+2bihi+ci=ci+1
这样，
c
c
c这个未知数也出来了。分段立方插值还有个设定，分段函数的一阶导必须连续可导，丝滑流畅。其实就是二阶导要相等。为了直观展示，我把上面例子的一阶导用绿色线条画出来：
设前一段函数的二阶导为
f
i
−
1
′
′
(
x
)
f_{i-1}''(x)
fi−1′′(x)，因为二阶导的一般方程是
f
′
′
(
x
)
=
6
a
(
x
−
x
i
)
+
2
b
f''(x)=6a(x-x_{i})+2b
f′′(x)=6a(x−xi)+2b，连接点处二阶导相等，我们又得到一个方程：
f
i
′
′
(
x
)
=
6
a
i
(
x
−
x
i
)
+
2
b
i
f
i
′
′
(
x
i
)
=
f
i
−
1
′
′
(
x
i
)
2
b
i
=
6
a
i
−
1
(
x
i
−
x
i
−
1
)
+
2
b
i
−
1
b
i
=
3
a
i
−
1
h
i
−
1
+
b
i
−
1
f_i''(x)=6a_i(x-x_{i})+2b_i\ f_i''(x_i)=f_{i-1}''(x_i)\ 2b_i=6a_{i-1}(x_i-x_{i-1})+2b_{i-1}\ b_i=3a_{i-1}h_{i-1}+b_{i-1}\
fi′′(x)=6ai(x−xi)+2bifi′′(xi)=fi−1′′(xi)2bi=6ai−1(xi−xi−1)+2bi−1bi=3ai−1hi−1+bi−1
也等同于：
a
i
−
1
=
b
i
−
b
i
−
1
3
h
i
−
1
a_{i-1}=frac{b_i-b_{i-1}}{3h_{i-1}}
ai−1=3hi−1bi−bi−1
与后一段在连接处二阶导相等，我们有：
f
i
′
′
(
x
i
+
1
)
=
f
i
+
1
′
′
(
x
i
+
1
)
6
a
i
(
x
i
+
1
−
x
i
)
+
2
b
i
=
6
a
i
+
1
(
x
i
+
1
−
x
i
+
1
)
+
2
b
i
+
1
6
a
i
h
i
+
2
b
i
=
2
b
i
+
1
3
a
i
h
i
+
b
i
=
b
i
+
1
f_{i}''(x_{i+1})=f_{i+1}''(x_{i+1})\ 6a_{i}(x_{i+1}-x_{i})+2b_{i}=6a_{i+1}(x_{i+1}-x_{i+1})+2b_{i+1}\ 6a_{i}h_{i}+2b_{i}=2b_{i+1}\ 3a_{i}h_{i}+b_{i}=b_{i+1}
fi′′(xi+1)=fi+1′′(xi+1)6ai(xi+1−xi)+2bi=6ai+1(xi+1−xi+1)+2bi+16aihi+2bi=2bi+13aihi+bi=bi+1
也等同于：
a
i
=
b
i
+
1
−
b
i
3
h
i
a_{i}=frac{b_{i+1}-b_{i}}{3h_{i}}
ai=3hibi+1−bi
总结下就是四个方程：
d
i
=
y
i
h
i
3
a
i
+
h
i
2
b
i
+
h
i
c
i
=
y
i
+
1
−
y
i
3
h
i
2
a
i
+
2
h
i
b
i
+
c
i
=
c
i
+
1
3
h
i
a
i
+
b
i
=
b
i
+
1
d_i=y_i\ h_i^3a_i+h_i^2b_i+h_ic_i=y_{i+1}-y_i\ 3h_{i}^2a_i+2h_{i}b_i+c_i=c_{i+1}\ 3h_{i}a_{i}+b_{i}=b_{i+1}
di=yihi3ai+hi2bi+hici=yi+1−yi3hi2ai+2hibi+ci=ci+13hiai+bi=bi+1
除去已经有解的
d
i
d_i
di，还剩三个方程，却有5个未知数，所以还差一些方程。
这么依赖前一段函数，也就是要求出前段函数的
a
,
b
,
c
a,b,c
a,b,c三个未知数，但是对于第一段函数，是没有前一段函数的，怎么办呢？
这就需要讲讲分段立方插值的两个重要设定了，起点二阶导为0与终点二阶导为0。为什么要这两个设定呢？我们可以用物理学的概念来理解，二阶导相当于加速度，而加速度源于力，这其实假设起点和终点是受力平衡的状态。在有些文章里是写错了的，我举个例子，国外著名的网站，就把首末端二阶导为0写成了首末端一阶导为0，这就不够严谨了，哈哈：
但是对于其他的英文网站，我推荐这篇加州伯克利写的文章。
对于第一段函数，二阶导方程如下：
f
′
′
(
x
0
)
=
0
6
a
0
(
x
0
−
x
0
)
+
2
b
0
=
0
2
b
0
=
0
f''(x_0)=0\ 6a_0(x_0-x_0)+2b_0=0\ 2b_0=0
f′′(x0)=06a0(x0−x0)+2b0=02b0=0
所以第一段函数
b
0
=
0
b_0=0
b0=0，这就减少一个了一个未知数。对于终点的二阶导为0，假设终点为
(
x
n
,
y
n
)
(x_n,y_n)
(xn,yn)我们计算一下：
f
′
′
(
x
n
)
=
6
a
n
−
1
(
x
n
−
x
n
−
1
)
+
2
b
n
−
1
=
0
3
a
n
−
1
(
x
n
−
x
n
−
1
)
+
b
n
−
1
=
0
f''(x_{n})=6a_{n-1}(x_{n}-{x_{n-1}})+2b_{n-1}=0\ 3a_{n-1}(x_{n}-{x_{n-1}})+b_{n-1}=0\ \
f′′(xn)=6an−1(xn−xn−1)+2bn−1=03an−1(xn−xn−1)+bn−1=0
哈哈，但是知道了这个方程，怎么与第一分段二阶导联系起来呢?
因为n个数据点，那么有
n
−
1
n-1
n−1个分段函数，除去
b
0
=
0
b_0=0
b0=0总共有
3
n
−
4
3n-4
3n−4个未知数。每一个段都有这个方程：
h
i
3
a
i
+
h
i
2
b
i
+
h
i
c
i
=
y
i
+
1
−
y
i
h_i^3a_i+h_i^2b_i+h_ic_i=y_{i+1}-y_i
hi3ai+hi2bi+hici=yi+1−yi
那么就现在就有了
n
−
1
n-1
n−1个方程。每一段和后一段的关系，有以下两个方程：
3
h
i
2
a
i
+
2
h
i
b
i
+
c
i
=
c
i
+
1
3
h
i
a
i
+
b
i
=
b
i
+
1
3h_{i}^2a_i+2h_{i}b_i+c_i=c_{i+1}\ 3h_{i}a_{i}+b_{i}=b_{i+1}
3hi2ai+2hibi+ci=ci+13hiai+bi=bi+1
那么就有
2
n
−
4
2n-4
2n−4个方程，加起来有
3
n
−
5
3n-5
3n−5个方程。再加上
3
a
n
−
1
(
x
n
−
x
n
−
1
)
+
b
n
−
1
=
0
3a_{n-1}(x_{n}-{x_{n-1}})+b_{n-1}=0
3an−1(xn−xn−1)+bn−1=0，这个方程，总共有
3
n
−
4
3n-4
3n−4个方程。未知数与方程的数量恰好相等，然后就可以解方程求出所有系数了。
大学时微积分最后一课就是泰勒级数和麦克劳林级数。大部分时候我们觉得这除了用来计算sin(x)和cos(x)的近似值外没什么用，而且真正到了计算机数值计算领域，已经是用切比雪夫近似值来计算了，不用泰勒级数了。原因很简单，泰勒级数一不精确，二性能低下。
我们对任意一段函数进行泰勒级数展开：
f
i
(
x
)
=
f
i
(
x
i
)
+
f
i
′
(
x
i
)
(
x
−
x
i
)
+
f
i
′
′
(
x
i
)
2
!
(
x
−
x
i
)
2
+
f
i
′
′
′
(
x
i
)
3
!
(
x
−
x
i
)
3
f_i(x)=f_i(x_i)+f'_i(x_i)(x-x_i)+frac{f''_i(x_i)}{2!}(x-x_i)^2+frac{f'''_i(x_i)}{3!}(x-x_i)^3
fi(x)=fi(xi)+fi′(xi)(x−xi)+2!fi′′(xi)(x−xi)2+3!fi′′′(xi)(x−xi)3
再看看我们的公式：
f
i
(
x
)
=
a
i
(
x
−
x
i
)
3
+
b
i
(
x
−
x
i
)
2
+
c
i
(
x
−
x
i
)
+
d
i
f_i(x)=a_i(x-x_i)^3+b_i(x-x_i)^2+c_i(x-x_i)+d_i
fi(x)=ai(x−xi)3+bi(x−xi)2+ci(x−xi)+di
对比一下，换个顺序，就会发现：
a
i
=
f
i
′
′
′
(
x
i
)
3
!
b
i
=
f
i
′
′
(
x
i
)
2
!
c
i
=
f
i
′
(
x
i
)
d
i
=
f
i
(
x
i
)
a_i=frac{f'''_i(x_i)}{3!}\ b_i=frac{f''_i(x_i)}{2!}\ c_i=f'_i(x_i)\ d_i=f_i(x_i)
ai=3!fi′′′(xi)bi=2!fi′′(xi)ci=fi′(xi)di=fi(xi)
哈哈，数学是不是特别有意思？
直接这样解方程，程序性能比较低，在实际生产中，我们使用
μ
mu
μ数列与z数列法。我们换个思路，找到一种更简便的方法，而不是强解
3
n
−
4
3n-4
3n−4元一次方程组。
μ
mu
μ数列的定义是：
μ
i
=
{
x
i
+
1
−
x
i
2
(
x
i
+
1
−
x
i
−
1
)
−
μ
i
−
1
(
x
i
−
x
i
−
1
)
,
i
∈
Z
+
0
,
i
=
0
mu_i=left{ begin{aligned} frac{x_{i+1}-x_{i}}{2(x_{i+1}-x_{i-1})-mu_{i-1}(x_i-x_{i-1})},i in mathbb{Z^+}\ 0,i=0 end{aligned} right.
μi=⎩⎪⎨⎪⎧2(xi+1−xi−1)−μi−1(xi−xi−1)xi+1−xi,i∈Z+0,i=0
z
z
z数列的定义是：
z
i
=
{
3
[
y
i
+
1
(
x
i
−
x
i
−
1
)
−
y
i
(
x
i
+
1
−
x
i
−
1
)
+
y
i
−
1
(
x
i
+
1
−
x
i
)
]
(
x
i
−
x
i
−
1
)
(
x
i
+
1
−
x
i
)
−
(
x
i
−
x
i
−
1
)
z
i
−
1
2
(
x
i
+
1
−
x
i
−
1
)
−
μ
i
−
1
(
x
i
−
x
i
−
1
)
,
i
∈
Z
+
0
,
i
=
0
z_i=left{ begin{aligned} frac { frac{3[y_{i+1}(x_{i}-x_{i-1})-y_i(x_{i+1}-x_{i-1})+y_{i-1}(x_{i+1}-x_i)]} {(x_i-x_{i-1})(x_{i+1}-x_i)} -(x_i-x_{i-1})z_{i-1}}{2(x_{i+1}-x_{i-1})-mu_{i-1}(x_i-x_{i-1})},i in mathbb{Z^+}\ 0,i=0 end{aligned} right.
zi=⎩⎪⎨⎪⎧2(xi+1−xi−1)−μi−1(xi−xi−1)(xi−xi−1)(xi+1−xi)3[yi+1(xi−xi−1)−yi(xi+1−xi−1)+yi−1(xi+1−xi)]−(xi−xi−1)zi−1,i∈Z+0,i=0
因为公式过于复杂，所以使用我们之前定义的
h
i
=
x
i
+
1
−
x
i
h_i=x_{i+1}-x_{i}
hi=xi+1−xi，将公式改写一下：
μ
i
=
{
h
i
2
(
x
i
+
1
−
x
i
−
1
)
−
μ
i
−
1
h
i
−
1
,
i
∈
Z
+
0
,
i
=
0
z
i
=
{
3
[
y
i
+
1
h
i
−
1
−
y
i
(
x
i
+
1
−
x
i
−
1
)
+
y
i
−
1
h
i
]
h
i
−
1
h
i
−
h
i
−
1
z
i
−
1
2
(
x
i
+
1
−
x
i
−
1
)
−
μ
i
−
1
h
i
−
1
,
i
∈
Z
+
0
,
i
=
0
mu_i=left{ begin{aligned} frac{h_{i}}{2(x_{i+1}-x_{i-1})-mu_{i-1}h_{i-1}},i in mathbb{Z^+}\ 0,i=0 end{aligned} right.\ z_i=left{ begin{aligned} frac { frac{3[y_{i+1}h_{i-1}-y_i(x_{i+1}-x_{i-1})+y_{i-1}h_i]} {h_{i-1}h_i} -h_{i-1}z_{i-1}}{2(x_{i+1}-x_{i-1})-mu_{i-1}h_{i-1}},i in mathbb{Z^+}\ 0,i=0 end{aligned} right.
μi=⎩⎪⎨⎪⎧2(xi+1−xi−1)−μi−1hi−1hi,i∈Z+0,i=0zi=⎩⎪⎨⎪⎧2(xi+1−xi−1)−μi−1hi−1hi−1hi3[yi+1hi−1−yi(xi+1−xi−1)+yi−1hi]−hi−1zi−1,i∈Z+0,i=0
搞这么复杂有什么用？哈哈，因为：
b
n
−
1
=
z
n
−
1
b_{n-1}=z_{n-1}
bn−1=zn−1
并且有：
b
i
=
z
i
−
μ
i
b
i
+
1
,
i
<
n
−
1
b_i=z_i-mu_ib_{i+1},i
首先简单实现一个霍纳多项式
# _*_ coding:utf-8 _*_
class HornerPolynomial:
def __init__(self, items: list):
self.__polynomial = items
@property
def polynomial(self):
return self.__polynomial
# 本身是个函数
def __call__(self, *args, **kwargs):
result = 0
for coefficient in self.polynomial:
result = result * args[0] + coefficient
return result
def __str__(self):
result = ''
first = True
n= len(self.polynomial)
for i, coefficient in enumerate(self.polynomial):
if coefficient == 0:
continue
if first:
first = False
else:
result += '+'
if coefficient != 1:
result += str(coefficient) + '*'
order = n - i - 1
if order > 0:
result += 'x'
else:
result += '1'
if order > 1:
result += '**' + str(order)
return result
再实现立方曲线插值：
# _*_ coding:utf-8 _*_
from com.youngthing.mathalgorithm.interpolation.horner_polynominal import HornerPolynomial
class PeacewisePolynomials:
def __init__(self, konts, polynomials):
self.__konts = konts
self.__polynomials = polynomials
def __call__(self, *args, **kwargs):
x = args[0]
for i, e in enumerate(self.__polynomials):
if x >= self.__konts[i] and x <= self.__konts[i + 1]:
return e(x - self.__konts[i])
return None
def __str__(self):
x = ''
for i, e in enumerate(self.__polynomials):
s = '' + str(e) + ',[' + str(self.__konts[i]) + ',' + str(self.__konts[i + 1]) + ']n'
x += s
return x
def getx(value):
return value[0]
def interpolate(values):
values = sorted(values, key=getx)
n = len(values) - 1
h = [0] * n
for i in range(0, n):
h[i] = getx(values[i + 1]) - getx(values[i])
mu = [0] * n
z = [0] * (n + 1)
x = [getx(e) for e in values]
y = [e[1] for e in values]
for i in range(1, n):
g = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1]
mu[i] = h[i] / g
temp = 3 * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1]) + y[i - 1] * h[i])
z[i] = (temp / (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g
c = [0] * n
b = [0] * (n + 1)
a = [0] * n
for j in range(n - 1, -1, -1):
b[j] = z[j] - mu[j] * b[j + 1];
c[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (b[j + 1] + 2 * b[j]) / 3
a[j] = (b[j + 1] - b[j]) / (3 * h[j])
polys = [HornerPolynomial([a[i], b[i], c[i], y[i]]) for i in range(0, n)]
return PeacewisePolynomials(x, polys)
代码测试
from com.youngthing.mathalgorithm.interpolation.cubic_spline import *
if __name__ == '__main__':
values=[
[1,1],
[2,2],
[3,1],
[4, -2],
[5, 1],
[6, 2],
[7, 1],
[8, -2],
[9, 1],
]
#测试数据
cubic = interpolate(values)
print(cubic(1.5))
print(cubic)
测试结果完全正确：
1.6005154639175259 -0.26804123711340205*x**3+1.268041237113402*x+1,[1,2] -0.6597938144329897*x**3+-0.8041237113402062*x**2+0.463917525773196*x+2*1,[2,3] 2.9072164948453607*x**3+-2.783505154639175*x**2+-3.1237113402061856*x+1,[3,4] -2.9690721649484533*x**3+5.9381443298969065*x**2+0.030927835051546726*x+-2*1,[4,5] 0.9690721649484534*x**3+-2.9690721649484533*x**2+3.0*x+1,[5,6] -0.9072164948453607*x**3+-0.06185567010309301*x**2+-0.030927835051546282*x+2*1,[6,7] 2.6597938144329896*x**3+-2.783505154639175*x**2+-2.8762886597938144*x+1,[7,8] -1.731958762886598*x**3+5.195876288659794*x**2+-0.463917525773196*x+-2*1,[8,9]
