分段低次插值
分段低次插值
为什么需要分段低次插值
之前拉格朗日插值时，可以发现，当插值节点很多时，会发生龙格现象，即在数据点的两端边缘出，多项式值与被插值的原函数差值很大。
解决这种现象可以通过所谓的分段低次插值来逼近原函数，所谓分段低次插值，就是分段低次多项式，对每两个数据点之间的小区间进行低次多项式插值。
常用的分段低次插值有：分段线性插值，分段三次Hermite插值，三次样条插值。
分段线性插值
分段线性插值很简单，实际上就是将所有数据点用折线连起来，得到的分段函数就是分段线性插值。但很明显，这样插值得到的函数在数据点处一阶导数不存在。
分段三次Hermite插值
同样很简单，只要在任意一个小区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk,xk+1]进行两点三次Hermite插值即可。
分段三次Hermite插值虽然更复杂，但保证了插值函数一阶连续可导，更光滑。
三次样条插值
三次样条插值对插值函数光滑性要求更高：
记插值函数 S ( x ) S(x) S(x)，要求：
- S ( x ) S(x) S(x)整体二阶连续可导。
- 每个小区间上都是三次多项式。
现在来看看如何求三次样条插值，首先，假定数据点 a = x 0 < x 1 < . . . < x n = b a=x_0\lt x_1\lt ...\lt x_n=b a=x0<x1<...<xn=b，插值数据点的原函数为 f ( x ) f(x) f(x)，则必有以下条件：
S ( x k ) = f ( x k ) , k = 0 , 1 , 2 , . . . , n S(x_k)=f(x_k),k=0,1,2,...,n S(xk)=f(xk),k=0,1,2,...,n
S ′ ( x k − ) = S ′ ( x k + ) , k = 1 , 2 , . . . , n − 1 S'(x_k^-)=S'(x_k^+),k=1,2,...,n-1 S′(xk−)=S′(xk+),k=1,2,...,n−1
S ′ ′ ( x k − ) = S ′ ′ ( x k + ) , k = 1 , 2 , . . . , n − 1 S''(x_k^-)=S''(x_k^+),k=1,2,...,n-1 S′′(xk−)=S′′(xk+),k=1,2,...,n−1
实际上，如果先将区间分段，共分为 n n n段，每段有 4 4 4个待定系数，需要 4 n 4n 4n个方程。而上面三个方程实际上是 4 n − 2 4n-2 4n−2个方程。费解方程个数的话只需要注意到：第一个式子实际上是 S ( x k + ) = f ( x k ) , k = 0 , 1 , . . . , n − 1 S(x_k^+)=f(x_k),k=0,1,...,n-1 S(xk+)=f(xk),k=0,1,...,n−1和 S ( x k − ) = f ( x k ) , k = 1 , . . . , n S(x_k^-)=f(x_k),k=1,...,n S(xk−)=f(xk),k=1,...,n。所以总方程数目是 2 n + ( n − 1 ) + ( n − 1 ) = 4 n − 2 2n+(n-1)+(n-1)=4n-2 2n+(n−1)+(n−1)=4n−2
现在还缺少两个方程，因为上面后两个方程里， k k k不能取 0 0 0和 n − 1 n-1 n−1。
想要确定插值函数，需要额外两个方程才行，这引出了边界条件。
第一类边界条件
给定两端一阶导数值：
S ′ ( x 0 + ) = f 0 ′ S'(x_0^+)=f_0' S′(x0+)=f0′
S ′ ( x n − ) = f n ′ S'(x_n^-)=f_n' S′(xn−)=fn′
第二类边界条件
给定两端二阶导数值
S ′ ′ ( x 0 + ) = f 0 ′ ′ S''(x_0^+)=f_0'' S′′(x0+)=f0′′
S ′ ′ ( x n − ) = f n ′ ′ S''(x_n^-)=f_n'' S′′(xn−)=fn′′
其中， f 0 ′ ′ = f n ′ ′ = 0 f_0''=f_n''=0 f0′′=fn′′=0称为自然边界条件。
第三类边界条件
周期边界条件，即 x n − x 0 x_n-x_0 xn−x0是一个周期，那么就有：
S ( x 0 ) = S ( x n ) , S ′ ( x 0 + ) = S ′ ( x n − ) , S ′ ′ ( x 0 + ) = S ′ ′ ( x n − ) S(x_0)=S(x_n),S'(x_0^+)=S'(x_n^-),S''(x_0^+)=S''(x_n^-) S(x0)=S(xn),S′(x0+)=S′(xn−),S′′(x0+)=S′′(xn−)$
推导样条插值的方法
首先，由于 S ( x ) S(x) S(x)二阶连续可导，且是三次多项式，那么 S ′ ′ ( x ) S''(x) S′′(x)实际上是分段线性函数。可以假设每一段的斜率是 M k , k = 0 , 1 , . . . , n M_k,k=0,1,...,n Mk,k=0,1,...,n
在区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk,xk+1]里，必有： S k ′ ′ ( x ) = x k + 1 − x x k + 1 − x k M k + x − x k x k + 1 − x k M k + 1 S_k''(x)=\dfrac{x_{k+1}-x}{x_{k+1}-x_k}M_k+\dfrac{x-x_k}{x_{k+1}-x_k}M_{k+1} Sk′′(x)=xk+1−xkxk+1−xMk+xk+1−xkx−xkMk+1
积分两次得： S k ( x ) = ( x k + 1 − x ) 3 6 ( x k + 1 − x k ) M k + ( x − x k ) 3 6 ( x k + 1 − x k ) M k + 1 + c 1 x + c 2 S_k(x)=\dfrac{{(x_{k+1}-x)}^3}{6(x_{k+1-x_k})}M_k+\dfrac{{(x-x_k)}^3}{6(x_{k+1-x_k})}M_{k+1}+c_1x+c_2 Sk(x)=6(xk+1−xk)(xk+1−x)3Mk+6(xk+1−xk)(x−xk)3Mk+1+c1x+c2
积分两次，得到两个待定常数 c 1 , c 2 c_1,c_2 c1,c2
这两个未知数只需要两个方程就可以确定，那么直接带入 S k ( x k ) = y k , S k ( x k + 1 ) = y k + 1 S_k(x_k)=y_k,S_{k}(x_{k+1})=y_{k+1} Sk(xk)=yk,Sk(xk+1)=yk+1就可以得出结果。
接着只需要确定所有的 M k M_k Mk就可以。
这里用Mathematica计算 x 0 , x 1 x_0,x_1 x0,x1上的结果，表达式还是相当复杂的。
{ { c1 → − − 1 6 M0 ( x1 − x0 ) 2 + 1 6 M1 ( x1 − x0 ) 2 + y0 − y1 x1 − x0 , c2 → − − M0 x0 2 x1 + 2 M0x0 x1 2 − M0 x1 3 + M1 x0 3 − 2 M1 x0 2 x1 + M1x0 x1 2 − 6 x0y1 + 6 x1y0 6 ( x0 − x1 ) } } \left\{\left\{\text{c1}\to -\dfrac{-\frac{1}{6} \text{M0} (\text{x1}-\text{x0})^2+\frac{1}{6} \text{M1} (\text{x1}-\text{x0})^2+\text{y0}-\text{y1}}{\text{x1}-\text{x0}},\\\text{c2}\to -\dfrac{-\text{M0} \text{x0}^2 \text{x1}+2 \text{M0} \text{x0} \text{x1}^2-\text{M0} \text{x1}^3+\text{M1} \text{x0}^3-2 \text{M1} \text{x0}^2 \text{x1}+\text{M1} \text{x0} \text{x1}^2-6 \text{x0} \text{y1}+6 \text{x1} \text{y0}}{6 (\text{x0}-\text{x1})}\right\}\right\} {{c1→−x1−x0−61M0(x1−x0)2+61M1(x1−x0)2+y0−y1,c2→−6(x0−x1)−M0x02x1+2M0x0x12−M0x13+M1x03−2M1x02x1+M1x0x12−6x0y1+6x1y0}}
不过课件里提供了一种化简，直接抄过来：
S k ( x ) = S k ( x ) = ( x k + 1 − x ) 3 6 ( x k + 1 − x k ) M k + ( x − x k ) 3 6 ( x k + 1 − x k ) M k + 1 + ( y k − M k ( x k + 1 − x k ) 2 6 x k + 1 − x x k + 1 − x k ) + ( y k + 1 − M k + 1 ( x k + 1 − x k ) 2 6 x − x k x k + 1 − x k ) S_k(x)=S_k(x)=\dfrac{{(x_{k+1}-x)}^3}{6(x_{k+1-x_k})}M_k+\dfrac{{(x-x_k)}^3}{6(x_{k+1-x_k})}M_{k+1}+(y_k-\dfrac{M_k{(x_{k+1}-x_k)}^2}{6}\dfrac{x_{k+1}-x}{x_{k+1}-x_k})+\\(y_{k+1}-\dfrac{M_{k+1}{(x_{k+1}-x_k)}^2}{6}\dfrac{x-x_k}{x_{k+1}-x_k}) Sk(x)=Sk(x)=6(xk+1−xk)(xk+1−x)3Mk+6(xk+1−xk)(x−xk)3Mk+1+(yk−6Mk(xk+1−xk)2xk+1−xkxk+1−x)+(yk+1−6Mk+1(xk+1−xk)2xk+1−xkx−xk)
对这个表达式应用方程： S k − 1 ′ ( x k − ) = S k ′ ( x k + ) S_{k-1}'(x_k^-)=S_k'(x_k^+) Sk−1′(xk−)=Sk′(xk+)可得到方程：
μ k M k − 1 + 2 M k + λ k M k + 1 = d k \mu_kM_{k-1}+2M_k+\lambda_kM_{k+1}=d_k μkMk−1+2Mk+λkMk+1=dk
其中： μ k = h k − 1 h k − 1 + h k , h k = x k + 1 − x k , λ k = h k h k − 1 + h k , d k = 6 f [ x k − 1 , x k , x k + 1 ] \mu_k=\dfrac{h_{k-1}}{h_{k-1}+h_k},h_k=x_{k+1}-x_k,\lambda_k=\dfrac{h_k}{h_{k-1}+h_k},\\d_k=6f[x_{k-1},x_k,x_{k+1}] μk=hk−1+hkhk−1,hk=xk+1−xk,λk=hk−1+hkhk,dk=6f[xk−1,xk,xk+1]
这里有 n − 1 n-1 n−1个方程，加上两个边界条件即可以确定所有的 M k M_k Mk。
第一类边界条件解法
带入边界条件 S ′ ( x 0 + ) = f 0 ′ , S ′ ( x n − ) = f n ′ S'(x_0^+)=f_0',S'(x_n^-)=f_n' S′(x0+)=f0′,S′(xn−)=fn′到原始的表达式里，可以得到， 2 M 0 + M 1 = d 0 , M n − 1 + 2 M n = d n 2M_0+M_1=d_0,M_{n-1}+2M_n=d_n 2M0+M1=d0,Mn−1+2Mn=dn
注意这里的 d n , d 0 d_n,d_0 dn,d0需要从原始表达式里直接求得，不能用 d k = 6 f [ x k − 1 , x k , x k + 1 ] d_k=6f[x_{k-1},x_k,x_{k+1}] dk=6f[xk−1,xk,xk+1]
由此，可以得到方程：
[ 2 1 μ 1 2 λ 1 μ 2 2 λ 2 ⋱ ⋱ ⋱ μ n − 1 2 λ n − 1 1 2 ] [ M 0 M 1 M 2 ⋮ M n − 1 M n ] = [ d 0 d 1 d 2 ⋮ d n − 1 d n ] \begin{bmatrix}2& 1& & &\\\mu_1&2&\lambda_1&&\\&\mu_2&2&\lambda_2&\\&&\ddots&\ddots&\ddots&\\&&&\mu_{n-1}&2&\lambda_{n-1}\\&&&&1&2\end{bmatrix}\begin{bmatrix}M_0\\M_1\\M_2\\\vdots\\M_{n-1}\\M_n\end{bmatrix}=\begin{bmatrix}d_0\\d_1\\d_2\\\vdots\\d_{n-1}\\d_n\end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎡2μ112μ2λ12⋱λ2⋱μn−1⋱21λn−12⎦⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎡M0M1M2⋮Mn−1Mn⎦⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎡d0d1d2⋮dn−1dn⎦⎥⎥⎥⎥⎥⎥⎥⎤
第二类边界条件解法
第二类边界条件相当于直接给出了 M 0 , M n M_0,M_n M0,Mn，那便直接用边界条件带入 μ k M k − 1 + 2 M k + λ k M k + 1 = d k \mu_kM_{k-1}+2M_k+\lambda_kM_{k+1}=d_k μkMk−1+2Mk+λkMk+1=dk确定 d 1 , d n − 1 d_1,d_{n-1} d1,dn−1这两个方程，可以解得，第二类边界条件的方程为：
[ 2 λ 1 μ 2 2 λ 2 ⋱ ⋱ ⋱ μ n − 2 2 λ n − 2 μ n − 1 2 ] [ M 1 M 2 ⋮ M n − 2 M n − 1 ] = [ d 1 − μ 1 M 0 d 2 ⋮ d n − 2 d n − 1 − λ n − 1 M n ] \begin{bmatrix}2&\lambda_1&&\\\mu_2&2&\lambda_2&\\&\ddots&\ddots&\ddots&\\&&\mu_{n-2}&2&\lambda_{n-2}\\&&&\mu_{n-1}&2\end{bmatrix}\begin{bmatrix}M_1\\M_2\\\vdots\\M_{n-2}\\M_{n-1}\end{bmatrix}=\begin{bmatrix}d_1-\mu_1M_0\\d_2\\\vdots\\d_{n-2}\\d_{n-1}-\lambda_{n-1}M_n\end{bmatrix} ⎣⎢⎢⎢⎢⎡2μ2λ12⋱λ2⋱μn−2⋱2μn−1λn−22⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡M1M2⋮Mn−2Mn−1⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡d1−μ1M0d2⋮dn−2dn−1−λn−1Mn⎦⎥⎥⎥⎥⎥⎤
第三类边界条件解法
第三类边界条件，即周期性边界条件，可以理解为 0 , 1 , . . . , n − 1 , 0 / n , 1 0,1,...,n-1,0/n,1 0,1,...,n−1,0/n,1构成了一个闭环，相当于指定了 M 0 = M n M_0=M_n M0=Mn。那么同样只需要看看 0 / n 0/n 0/n附近的方程形式就好。
选择 k = n k=n k=n，那么 μ n M n − 1 + 2 M n + λ n M n + 1 = d n \mu_nM_{n-1}+2M_{n}+\lambda_nM_{n+1}=d_n μnMn−1+2Mn+λnMn+1=dn，可以将 M n + 1 M_{n+1} Mn+1视为 M 1 M_1 M1
选择 k = 1 k=1 k=1,那么 μ 1 M 0 + 2 M 1 + λ 1 M 2 = d 1 \mu_1M_{0}+2M_{1}+\lambda_1M_{2}=d_1 μ1M0+2M1+λ1M2=d1，可以将 M 0 M_0 M0视为 M n M_n Mn
那么方程为：
[ 2 λ 1 μ 1 μ 2 2 λ 2 ⋱ ⋱ ⋱ μ n − 1 2 λ n − 1 λ n μ n 2 ] [ M 1 M 2 ⋮ M n − 1 M n ] = [ d 1 d 2 ⋮ d n − 1 d n ] \begin{bmatrix}2&\lambda_1&&&\mu_1\\\mu_2&2&\lambda_2&\\&\ddots&\ddots&\ddots&\\&&\mu_{n-1}&2&\lambda_{n-1}\\\lambda_n&&&\mu_{n}&2\end{bmatrix}\begin{bmatrix}M_1\\M_2\\\vdots\\M_{n-1}\\M_{n}\end{bmatrix}=\begin{bmatrix}d_1\\d_2\\\vdots\\d_{n-1}\\d_{n}\end{bmatrix} ⎣⎢⎢⎢⎢⎡2μ2λnλ12⋱λ2⋱μn−1⋱2μnμ1λn−12⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡M1M2⋮Mn−1Mn⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡d1d2⋮dn−1dn⎦⎥⎥⎥⎥⎥⎤
解这三类边界条件的方程叫做三弯矩方程，它们都有唯一解。
三次样条插值示例
对函数 f ( x ) = 1 1 + x 2 f(x)=\dfrac{1}{1+x^2} f(x)=1+x21进行插值，区间 [ − 5 , 5 ] [-5,5] [−5,5]，取11个等距节点，比较三次样条插值和10次拉格朗日插值的函数。
为了方便，直接导入了前面博客里的拉格朗日插值的程序。https://blog.csdn.net/lan_777/article/details/111146113
import numpy as np
import matplotlib.pyplot as plt
# 对函数1/(1+x^2)插值，区域[-5,5]，取11个数据点
# 使用第二类边界条件来解
n = 11 - 1
xdata = np.linspace(-5, 5, n + 1)
ydata = 1 / (1 + xdata ** 2)
M0 = 37 / 4396
Mn = 37 / 4396
# 计算三弯矩方程所需的系数矩阵
# hk的k从0到n-1
hk = np.diff(xdata)
# uk=h_{k-1}/(h_{k-1}+h_k),k从1到n-1
uk = hk[1:] / (hk[1:] + hk[:-1])
lambdak = hk[:-1] / (hk[1:] + hk[:-1])
# 计算dk，k同样从1到n-1
dk = 6 * np.diff(np.diff(ydata) / np.diff(xdata)) / np.diff(xdata)[:-1]
# 计算M
# 创建系数矩阵
factor_mat = 2 * np.eye(n - 1) + np.diag(lambdak[:-1], k=1) + np.diag(uk[1:], k=-1)
dk[0] = dk[0] - uk[0] * M0
dk[-1] = dk[-1] - lambdak[-1] * Mn
# 计算二阶导数
Mk = np.linalg.inv(factor_mat) @ dk
# print(type(Mk))
Mk = np.insert(Mk, 0, M0)
# print((Mk))
Mk = np.insert(Mk, len(Mk), Mn)
# print((Mk))
# 分段求函数
def sx(mk, hk, xdata, ydata, plot_x):
if np.max(plot_x) > np.max(xdata) or np.min(plot_x) < np.min(xdata):
print("输出范围需要在节点内。")
return False
plot_y = np.array([])
for i in range(n):
# print((plot_x < xdata[i+1]) & (plot_x >= xdata[i]))
segment_x = plot_x[((plot_x < xdata[i + 1]) & (plot_x >= xdata[i]))]
if i == n-1:
segment_x = plot_x[((plot_x <= xdata[i + 1]) & (plot_x >= xdata[i]))]
segment_y = (mk[i + 1] - mk[i]) / (6 * hk[i]) * (segment_x - xdata[i]) ** 3 + \
mk[i] / 2 * (segment_x - xdata[i]) ** 2 + \
((ydata[i + 1] - ydata[i]) / hk[i] - (hk[i] * (mk[i + 1] + 2 * mk[i])) / 6) * \
(segment_x - xdata[i]) + ydata[i]
plot_y = np.insert(plot_y, len(plot_y), segment_y)
return plot_y
# 导入之前的拉格朗日插值程序，在下面直接绘制
import lagrange
plot_x = np.linspace(-5, 5, 400)
plot_y = sx(Mk, hk, xdata, ydata, plot_x)
# print(plot_y)
p1, = plt.plot(plot_x, plot_y)
p2, = plt.plot(plot_x, 1 / (1 + plot_x ** 2))
p4, = plt.plot(lagrange.x_plot,lagrange.l_matrix @ lagrange.ydata)
p3 = plt.scatter(xdata, ydata)
plt.legend([p1, p2, p4], ["Spline", "f(x)", "Lagrange"])
plt.show()
运行结果：
这里可以看出，三次样条插值的结果要明显优秀于高次拉格朗日插值。
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
