【数学建模笔记 09-1】数学建模的插值

09-1. 插值

定义

插值：求过已知有限个数据点的近似函数。

插值方法

拉格朗日插值

用多项式作为插值工具，称代数插值。

已知函数 $f(x)$ 在区间 $[a,b]$ 上 $n+1$ 个不同点 $x_0,x_1,\dots ,x_n$ 处的函数值 $y_i=f(x_i),i=0,1,\dots ,n$，求一个至多 $n$ 次多项式
$${\phi }_n(x)=a_0+a_{1}x+\cdots +a_n{x}^n$$
使其在给定点处与 $f(x)$ 同值，即
$${\phi }_n(x_i)=f(x_i)=y_i,i=0,1,\dots ,n.$$
根据插值条件，构造方程组
$$\{\begin{array}{rl}& a_0+a_1{x}_0+a_2{x}_0^2+\cdots +a_n{x}_0^n=y_0,\\ & a_0+a_1{x}_2+a_2{x}_2^2+\cdots +a_n{x}_2^n=y_2,\\ & ⋮\\ & a_0+a_1{x}_n+a_2{x}_n^2+\cdots +a_n{x}_n^n=y_n.\end{array}$$
记变量矩阵为 $A$，则
$$\text{det}(A)=\left|\begin{array}{ccccc}1& x_0& x_0^2& \dots & x_0^n\\ 1& x_1& x_1^2& \dots & x_1^n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_n& x_n^2& \dots & x_n^n\end{array}\right|$$
是范德蒙特行列式。当 $x_0,x_1,\dots ,x_n$ 互不相同时，值不为 0，因此方程组有唯一解。

构造基函数
$$l_i(x)=\frac{(x-x_0)\dots (x-x_{i-1})(x-x_{i+1})\dots (x-x_n)}{(x_i-x_0)\dots (x_i-x_{i-1})(x_i-x_{i+1})\dots (x_i-x_n)}$$

$$=\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j},i=0,1,\dots ,n.$$

满足
$$l_i(x_j)=\{\begin{array}{rl}& 1,j=i,\\ & 0,j\ne i.\end{array}$$
这意味着，$l_i(x_j)$ 在 $x_j$ 上值为 1，在其他点上值都为 0。

于是有 $y_i{l}_i(x_j)$ 在 $x_j$​ 上值为 $y_j$，在其他点上值都为 0。

因此有多项式
$$L_n(x)=\sum _{i=0}^n{y}_i{l}_i(x)=\sum _{i=0}^n{y}_i(\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j})$$
满足
$$L_n(x_j)=0+0+\cdots +y_j+\cdots +0=y_j.$$
这意味着多项式 $L_n(x)$ 在 $x_0,x_1,\dots ,x_n$ 各点上，都有 $L_n(x_j)=y_j$，满足唯一的多项式插值条件。

例子

已知有某个二次多项式函数 $f$，在三个点上的取值为：

• $f(4)=10,$
• $f(5)=5.25,$
• $f(6)=1,$

求 $f(18)$ 的值。

首先写出基函数
$$l_0(x)=\frac{(x-5)(x-6)}{(4-5)(4-6)},$$

$$l_1(x)=\frac{(x-4)(x-6)}{(5-4)(4-6)},$$

$$l_2(x)=\frac{(x-4)(x-5)}{(6-4)(6-5)}.$$

于是有
$$L(x)=f(4)l_0(x)+f(5)l_1(x)+f_6{l}_2(x)$$

$$=\frac{1}{4}(x^2-28x+136).$$

则 $f(18)=L(18)=-11$。

分段线性插值

即用折线代替曲线，多项式为
$$P_1(x)=\frac{x-x_i}{x_{i+1}-x_i}{y}_{i+1}+\frac{x-x_{i+1}}{x_i-x_{i+1}},$$

$$x\in [x_i,x_{i+1}],i=0,1,\dots ,n-1.$$

分段二次插值

即用抛物线代替曲线，多项式为
$$P_2(x)=\frac{(x-x_{2i+1})(x-x_{2i+2})}{(x_{2i}-x_{2i+1})(x_{2i}-x_{2i+2})}{y}_{2i}$$

$$+\frac{(x-x_{2i})(x-x_{2i+2})}{(x_{2i+1}-x_{2i})(x_{2i+1}-x_{2i+2})}{y}_{2i+1}$$

$$+\frac{(x-x_{2i})(x-x_{2i+1})}{(x_{2i+2}-x_{2i+1})(x_{2i+2}-x_{2i+1})}{y}_{2i+2},$$

$$x\in [x_{2i},x_{2i+2}],i=0,1,2,\dots ,n-1.$$

牛顿插值

设函数 $f(x)$ 的一系列相异的结点 $x_0\le x_1\le \cdots \le x_n$，则称
$$\frac{f(x_i)-f(x_j)}{x_i-x_j},i\ne j$$
为 $f(x)$ 关于 $x_i,x_j$ 的一阶差商，记 $f[x_i,x_j]$。

同样的，称一阶差商的差商
$$\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}$$
为 $f(x)$​ 关于 $x_i,x_j,x_k$ 的二阶差商，记 $f[x_i,x_j,x_k]$。

同样有
$$f[x_0,x_1,\dots ,x_k]=\frac{f[x_0,\dots ,x_{k-1}]-f[x_1,\dots ,x_k]}{x_0-x_k}$$
为 $f(x)$ 关于 $x_0,x_1,\dots ,x_k$ 的 $k$ 阶差商。

根据差商定义，有
$$f(x)=f(x_0)+(x-x_0)f[x,x_0],$$

$$f[x,x_0]=f[x_0,x_1]+(x-x_1)f[x,x_0,x_1],$$

$$f[x,x_0,x_1]=f[x_0,x_1,x_2]+(x-x_2)f[x,x_0,x_1,x_2],$$

$$\dots$$

$$f[x,x_0,\dots ,x_{n-1}]=f[x_0,\dots ,x_n]+(x-x_n)f[x,\dots ,x_n]$$

即得
$$f(x)=f(x_0)+(x-x_0)f[x_0,x_1]+\dots$$

$$+(x-x_0)(x-x_1)\dots (x-x_{n-1})f[x_0,\dots ,x_n]$$

$$+(x-x_0)(x-x_1)\dots (x-x_n)f[x,x_0,\dots ,x_n].$$

记
$$N_n(x)=f(x_0)+(x-x_0)f[x_0,x_1]+\dots$$

$$+(x-x_0)(x-x_1)\dots (x-x_{n-1})f[x_0,\dots ,x_n],$$

$$R_n(x)=(x-x_0)(x-x_1)\dots (x-x_n)f[x,x_0,\dots ,x_n].$$

$N_n(x)$ 为 $f(x)$ 的 $n$ 次插值多项式，$R_n(x)$ 称牛顿插值余项。

三次样条插值

三次样条方程 $S(x)$ 满足以下条件：

• 在每个分段小区间 $[x_i,x_{i+1}]$​ 上，$S(x)=S_i(x)$ 都是一个三次方程；
• 满足插值条件，即 $S(x_i)=y_i,i=0,1,\dots ,n$；
• 曲线光滑，即 $S(x),S^{\prime }(x),S^{\prime \prime }(x)$ 连续。

则对于第 $i$ 个小区间 $[x_i,x_{i+1}]$ 有
$$S_i(x)=a_i+b_{i}x+c_i{x}^2+d_i{x}^3,$$
包含 $a_i,b_i,c_i,d_i$ 4 个未知数。

有对于点 $x_0,x_1,\dots ,x_n$，有 $n$​ 个小区间，因此要找出 $4n$ 个方程求解 $4n$ 个未知数。

为满足插值条件，除两个端点外，所有内部点都满足
$$S_i(x_{i+1})=y_{i+1},S_{i+1}(x_{i+1})=y_{i+1},$$

$$i=0,1,2,\dots ,n-2,$$

共 $2(n-1)$​ 个方程，

加上两个端点
$$S_0(x_0)=y_0,S_{n-1}(x_n)=y_n,$$
共 $2n$ 个方程。

为满足连续条件，有每个内部点对应两个方程的一阶导数和二阶导数相等，即
$$S_i^{\prime }(x_{i+1})=S_{i+1}^{\prime }(x_{i+1}),S_i^{\prime \prime }(x_{i+1})=S_{i+1}^{\prime \prime }(x_{i+1})$$

$$i=0,1,2,\dots ,n-2,$$

共 $2n-2$ 个方程，总共 $4n-2$ 个方程，还差两个。

最后两个方程通过边界条件得到：

• 自然边界：$S_0^{\prime \prime }(x_0)=S_{n-1}^{\prime \prime }(x_n)=0$​​；
• 固定边界：$S_0^{\prime \prime }(x_0)=A,S_{n-1}^{\prime \prime }(x_n)=B$​​；
• 非扭结边界：$S_0^{\prime \prime \prime }(x_0)=S_1^{\prime \prime \prime }(x_1),S_{n-1}^{\prime \prime \prime }(x_n)=S_{n-2}^{\prime \prime \prime }(x_{n-1})$。

选一种边界条件即可。

若选用自然条件，则构成方程组
{ S i ( x i ) = y i , i = 0 , 1 , … , n , S i ( x i + 1 ) = S i + 1 ( x i + 1 ) , i = 0 , 1 , … , n − 2 , S i ′ ( x i + 1 ) = S i + 1 ′ ( x i + 1 ) , i = 0 , 1 , … , n − 2 , S i ′ ′ ( x i + 1 ) = S i + 1 ′ ′ ( x i + 1 ) , i = 0 , 1 , … , n − 2 , S 0 ′ ′ ( x 0 ) = S n − 1 ′ ′ ( x n ) = 0. \left\{\begin{aligned} &S_i(x_i)=y_i,i=0,1,\dots,n,\\ &S_i(x_i+1)=S_{i+1}(x_{i+1}),i=0,1,\dots,n-2,\\ &S_i'(x_i+1)=S_{i+1}'(x_{i+1}),i=0,1,\dots,n-2,\\ &S_i''(x_i+1)=S_{i+1}''(x_{i+1}),i=0,1,\dots,n-2,\\ &S_0''(x_0)=S_{n-1}''(x_n)=0. \end{aligned}\right. ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​​Si​(xi​)=yi​,i=0,1,…,n,Si​(xi​+1)=Si+1​(xi+1​),i=0,1,…,n−2,Si′​(xi​+1)=Si+1′​(xi+1​),i=0,1,…,n−2,Si′′​(xi​+1)=Si+1′′​(xi+1​),i=0,1,…,n−2,S0′′​(x0​)=Sn−1′′​(xn​)=0.​
对方程组求解即可。

求解方程组的过程较繁琐，思路是将方程组转化为线性形式，此处不过多赘述，大致思路为，令
$$S_j(x)=a_j+b_j(x-x_j)+c_j(x-x_j{)}^2+d_j(x-x_3{)}^3$$
则 $a_{j+1}=S_j(x_{j+1})=S_{j+1}(x_{j_1})$​。

设 $h_j=x_{j+1}-x_j$，​​

则
$$a_{j+1}=a_j+b_j{h}_j+c_j{h}_j^2+d_j{h}_j^3,$$
于是有
$$b_{j+1}=b_j+2c_j{h}_j+3d_j{h}_j^2,$$

$$c_{j+1}=c_j+3d_j{h}_j.$$

得到线性方程组，求解即可。

Python 代码

拉格朗日插值

使用 scipy 对上述例题进行拉格朗日插值，代码如下：

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-25
# @ function: 使用 scipy 包进行拉格朗日插值

# %%

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)

输出如下：

poly =
       2
0.25 x - 7 x + 34
res = -11.0

可以看出，结果为
$$f(x)=0.25x^2-7x+34,f(18)=-11.$$
并在图上标注出来。

样条插值

使用 scipy 对 (1,4), (2,3), (5,7), (8,11), (9,5), (12,3), (15,13), (17,10) 进行一次、二次、三次样条插值，代码如下：

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-25
# @ function: 使用 scipy 包进行样条插值

# %%

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()

输出如下：

其中蓝线、黄线、绿线分别对应一次、二次、三次样条插值。