09-1. 插值
定义
插值：求过已知有限个数据点的近似函数。
插值方法
拉格朗日插值
用多项式作为插值工具，称代数插值。
已知函数 在区间 上 个不同点 处的函数值 ，求一个至多
次多项式 使其在给定点处与
同值，即 根据插值条件，构造方程组 记变量矩阵为 ，则 是范德蒙特行列式。当 互不相同时，值不为
0，因此方程组有唯一解。
构造基函数
满足 这意味着， 在
上值为 1，在其他点上值都为
0。
于是有 在 上值为 ，在其他点上值都为 0。
因此有多项式 满足 这意味着多项式 在
各点上，都有
，满足唯一的多项式插值条件。
例子
已知有某个二次多项式函数 ，在三个点上的取值为：
求 的值。
首先写出基函数
于是有
则 。
分段线性插值
即用折线代替曲线，多项式为
分段二次插值
即用抛物线代替曲线，多项式为
牛顿插值
设函数 的一系列相异的结点
，则称
为 关于 的一阶差商，记 。
同样的，称一阶差商的差商 为 关于 的二阶差商，记 。
同样有 为 关于 的 阶差商。
根据差商定义，有
即得
记
为 的 次插值多项式， 称牛顿插值余项。
三次样条插值
三次样条方程
满足以下条件：
- 在每个分段小区间
上，
都是一个三次方程；
- 满足插值条件，即 ；
- 曲线光滑，即 连续。
则对于第 个小区间 有 包含 4
个未知数。
有对于点 ，有
个小区间，因此要找出 个方程求解 个未知数。
为满足插值条件，除两个端点外，所有内部点都满足
共 个方程，
加上两个端点 共 个方程。
为满足连续条件，有每个内部点对应两个方程的一阶导数和二阶导数相等，即
共 个方程，总共 个方程，还差两个。
最后两个方程通过边界条件得到：
- 自然边界：；
- 固定边界：；
- 非扭结边界：。
选一种边界条件即可。
若选用自然条件，则构成方程组 对方程组求解即可。
求解方程组的过程较繁琐，思路是将方程组转化为线性形式，此处不过多赘述，大致思路为，令
则 。
设 ，
则 于是有
得到线性方程组，求解即可。
Python 代码
拉格朗日插值
使用 scipy 对上述例题进行拉格朗日插值，代码如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
|
import numpy as np
from scipy.interpolate import lagrange
x = np.array([4, 5, 6])
y = np.array([10, 5.25, 1])
poly = lagrange(x, y)
res = poly(18)
print('poly =\n', poly)
print('res =', res)
import matplotlib.pyplot as plt
x1 = [i for i in range(-5,20)]
y1 = []
for each in x1:
y1.append(poly(each))
plt.plot(x1, y1)
plt.scatter(x,y)
plt.scatter(18, res)
|
输出如下：
1
2
3
4
| poly =
2
0.25 x - 7 x + 34
res = -11.0
|
可以看出，结果为 并在图上标注出来。
样条插值
使用 scipy 对 (1,4), (2,3), (5,7), (8,11), (9,5), (12,3), (15,13),
(17,10) 进行一次、二次、三次样条插值，代码如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
|
import numpy as np
from scipy.interpolate import interp1d
x = np.array([1, 2, 5, 8, 9, 12, 15, 17])
y = np.array([4, 3, 7, 11, 5, 3, 13, 10])
p1 = interp1d(x, y, kind='linear')
p2 = interp1d(x, y, kind='quadratic')
p3 = interp1d(x, y, kind='cubic')
x1 = np.linspace(1, 17, 100)
y1 = p1(x1)
y2 = p2(x1)
y3 = p3(x1)
import matplotlib.pyplot as plt
plt.scatter(x,y)
plt.plot(x1,y1, label='linear')
plt.plot(x1,y2, label='quadratic')
plt.plot(x1,y3, label='cubic')
plt.legend()
|
输出如下：
其中蓝线、黄线、绿线分别对应一次、二次、三次样条插值。