关键字：Matalb；曲线拟合；高次病态特性；分段低次插值
系列文章目录
数值分析——拉格朗日插值
数值分析——牛顿插值多项式
数值分析——埃尔米特（Hermit）插值
数值分析——三次样条插值
前言
一般认为对于区间 [ a , b ] \left[a,b\right] [a,b]上给出的节点做插值多项式 L n ( x ) L_n\left( x \right) Ln(x)近似 f ( x ) f\left( x \right) f(x)，一般总认为 L n ( x ) L_n\left( x \right) Ln(x)的次数 n n n越高，插值函数逼近原函数的效果越好，但实际上并非如此。这是因为对于任意的插值节点，当 n → ∞ n\rightarrow \infty n→∞时， L n ( x ) L_n\left( x \right) Ln(x)不一定收敛于 f ( x ) f\left( x \right) f(x)。为了解决高次插值的问题，学者们提出了分段低次插值的方法，即将被拟合函数 f ( x ) f\left( x \right) f(x)分为多个小区间，在小区间上利用低次多项式进行插值拟合。
本文对一些复杂的流程与概念进行了简化与省略，如果想详细了解的小伙伴可以参考李庆扬老师的《数值分析 -第5版》的第2章（需要电子版请私戳）
一、理论推导
1.高次插值的病态特性
高次插值的病态特性指的是，在对某些函数
f
(
x
)
f\left( x \right)
f(x)进行拟合时，高次多项式的插值函数拟合效果反而不好的特性。例如针对函数：
f
(
x
)
=
1
x
2
+
1
(1)
f\left( x \right)=\frac{1}{x^2+1} \tag{1}
f(x)=x2+11(1)使用拉格朗日插值法，在插值节点
x
=
−
5
,
−
4
,
−
3
,
−
2
,
−
1
,
0
,
1
,
2
,
3
,
4
,
5
x=-5,-4,-3,-2,-1,0,1,2,3,4,5
x=−5,−4,−3,−2,−1,0,1,2,3,4,5拟合十次插值函数：
% 高次插值的动态特性
x=-5:0.05:5;
y=zeros(1,length(x)); % 原函数y值
xInArr=-5:1:5; % 插值节点x值
yInArr=zeros(1,length(xInArr)); % 插值节点y值
yOutArr=zeros(1,length(x)); % 插值函数y值
% 计算插值函数y值
for i=1:1:length(xInArr)
yInArr(i)=1/(xInArr(i)^2+1);
end
% 计算十次插值多项式的poly数组
L10=func_LagrangeInterpolation(xInArr,yInArr,11);
% 计算原函数与插值函数y值
for i=1:1:length(x)
y(i)=1/(x(i)^2+1);
yOutArr(i)=polyval(L10,x(i));
end
% 画图
plot(x,y,'k-',x,yOutArr,'r:','LineWidth',2);
xlabel('x');
ylabel('y');
legend('f(x)','L10(x)');
其中func_LagrangeInterpolation是拉格朗日插值法函数，具体内容可参考数值分析——拉格朗日插值中的代码。得到的对比结果：
其中红色的点线为十次拉格朗日插值函数 L 10 ( x ) L_{10}\left( x \right) L10(x)，黑色的曲线为原函数 f ( x ) f\left( x \right) f(x)，可以看到插值函数的拟合效果并不是很好。
2.分段线性插值
分段线性插值就是通过插值点用折现段连接起来逼近
f
(
x
)
f\left( x \right)
f(x)。若在已知节点
a
=
x
0
<
x
1
<
⋯
<
x
n
=
b
a=x_0<x_1<\cdots<x_n=b
a=x0<x1<⋯<xn=b，已知各点的函数值f_0,f_1,f_2,\cdots,f_{n-1},f_n，则可以设计一个插值函数：
I
h
(
x
)
=
x
−
x
k
+
1
x
k
−
x
k
+
1
f
k
+
x
−
x
k
x
k
+
1
−
x
k
f
k
+
1
,
x
k
≤
x
≤
x
k
+
1
,
k
=
0
,
1
,
⋯
,
n
−
1
(2)
I_h\left( x \right)= \frac{x-x_{k+1}}{x_{k}-x_{k+1}}f_k+ \frac{x-x_{k}}{x_{k+1}-x_{k}}f_{k+1}, \ \ x_k \le x \le x_{k+1}, \ \ k=0,1,\cdots,n-1 \tag{2}
Ih(x)=xk−xk+1x−xk+1fk+xk+1−xkx−xkfk+1, xk≤x≤xk+1, k=0,1,⋯,n−1(2)该插值函数满足：
- 函数属于在区间 [ a , b ] [a,b] [a,b]连续的函数空间： I h ( x ) ∈ C [ a , b ] I_h\left( x \right)\in C\left[a,b\right] Ih(x)∈C[a,b]；
- 插值函数在所有插值点的数值都等于原函数： I h ( x ) = f k ( k = 0 , 1 , ⋯ , n − 1 ) I_h\left( x \right)=f_k\left(k=0,1,\cdots,n-1\right) Ih(x)=fk(k=0,1,⋯,n−1)；
- 插值函数在每个小区间上是线性函数： I h ( x ) = k x + b I_h\left( x \right)=kx+b Ih(x)=kx+b。
3.分段三次Hermit插值
分段线性插值函数的导数是间断的，如果在插值节点
a
=
x
0
<
x
1
<
⋯
<
x
n
=
b
a=x_0<x_1<\cdots<x_n=b
a=x0<x1<⋯<xn=b上，不仅已知函数值
f
0
,
f
1
,
f
2
,
⋯
,
f
n
−
1
,
f
n
f_0,f_1,f_2,\cdots,f_{n-1},f_n
f0,f1,f2,⋯,fn−1,fn，而且已知函数导数值
f
0
′
,
f
1
′
,
f
2
′
,
⋯
,
f
n
−
1
′
,
f
n
′
f^{'}_0,f^{'}_1,f^{'}_2,\cdots,f^{'}_{n-1},f^{'}_n
f0′,f1′,f2′,⋯,fn−1′,fn′，则可以设计一个插值函数：
I
h
(
x
)
=
(
x
−
x
k
+
1
x
k
−
x
k
+
1
)
2
(
1
+
2
x
−
x
k
x
k
+
1
−
x
k
)
f
k
+
(
x
−
x
k
x
k
+
1
−
x
k
)
2
(
1
+
2
x
−
x
k
+
1
x
k
−
x
k
+
1
)
f
k
+
1
+
(
x
−
x
k
+
1
x
k
−
x
k
+
1
)
2
(
x
−
x
k
)
f
k
′
+
(
x
−
x
k
x
k
+
1
−
x
k
)
2
(
x
−
x
k
+
1
)
f
k
+
1
′
,
x
k
≤
x
≤
x
k
+
1
,
k
=
0
,
1
,
⋯
,
n
−
1
(3)
\begin{align} I_h\left( x \right)= &\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^2 \left( 1+2\frac{x-x_k}{x_{k+1}-x_{k}} \right)f_k+ \left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^2\left( 1+2\frac{x-x_{k+1}}{x_{k}-x_{k+1}} \right)f_{k+1} \\ &+\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^2\left( x-x_k \right)f^{'}_k+\left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^2\left( x-x_{k+1} \right)f^{'}_{k+1}, \ \ x_k \le x \le x_{k+1}, \ \ k=0,1,\cdots,n-1 \end{align} \tag{3}
Ih(x)=(xk−xk+1x−xk+1)2(1+2xk+1−xkx−xk)fk+(xk+1−xkx−xk)2(1+2xk−xk+1x−xk+1)fk+1+(xk−xk+1x−xk+1)2(x−xk)fk′+(xk+1−xkx−xk)2(x−xk+1)fk+1′, xk≤x≤xk+1, k=0,1,⋯,n−1(3)该插值函数满足：
- 插值函数属于在区间 [ a , b ] [a,b] [a,b]上一阶导数连续的函数空间： I h ( x ) ∈ C 1 [ a , b ] I_h\left( x \right)\in C^1\left[a,b\right] Ih(x)∈C1[a,b]
- 插值函数在所有插值点的数值都等于原函数： I h ( x ) = f k ( k = 0 , 1 , ⋯ , n − 1 ) I_h\left( x \right)=f_k\left(k=0,1,\cdots,n-1\right) Ih(x)=fk(k=0,1,⋯,n−1)
- 插值函数一阶导数在所有插值点的数值都等于原函数一阶导数： I h ′ ( x ) = f k ′ ( k = 0 , 1 , ⋯ , n − 1 ) I_h^{'}\left( x \right)=f_k^{'}\left(k=0,1,\cdots,n-1\right) Ih′(x)=fk′(k=0,1,⋯,n−1)
- 插值函数在每个小区间上是三次多项式： I h ( x ) = a x 3 + b x 2 + c x + d I_h\left( x \right)=ax^3+bx^2+cx+d Ih(x)=ax3+bx2+cx+d
二、MATLAB实现
1.分段线性插值
% 分段线性插值
function [yOutArr]=func_PiecewiseLinearInterpolation(xInArr,yInArr,xQuery)
% ************************************************************
% 识别误差
% ************************************************************
xInLength=length(xInArr);
yInLength=length(yInArr);
% 插值数组不匹配
if xInLength~=yInLength
error('Error:The xInArr %d and the yInArr %d are not equal in length.',xInLength,yInLength);
end
% ************************************************************
% 计算查询节点在插值函数上的函数值
% ************************************************************
% 只一个插值节点
if xInLength==1
yOutArr=yInArr;
return;
end
% 有多个插值节点
xQueryLength=length(xQuery);
yOutArr=zeros(1,xQueryLength);
pos=1;
for i=1:1:xQueryLength
for j=pos:1:xInLength-1
% 查询节点如果超出插值节点范围则报错
if xQuery(i)>xInArr(xInLength)||xQuery(i)<xInArr(1)
error('Error:The num of the xQuery should gainer than %d and lesser than %d.',xInArr(1),xInArr(xInLength));
end
% 定位查询数组在插值节点的位置
if (xQuery(i)>xInArr(j)&&xQuery(i)<xInArr(j+1))||(xQuery(i)==xInArr(j))
pos=j;
break;
end
end
% 计算查询节点在插值函数上的函数值
yOutArr(i)=polyval(func_LinearInterpolation(xInArr(pos),yInArr(pos),xInArr(pos+1),yInArr(pos+1)),xQuery(i));
end
end
2.分段Hermit插值
% 分段三次Hermit插值
function [yOutArr]=func_PiecewiseCubeHermitInterpolation(xInArr,yInArr,dyInArr,xQuery)
% ************************************************************
% 识别误差
% ************************************************************
xInLength=length(xInArr);
yInLength=length(yInArr);
% 插值数组不匹配
if xInLength~=yInLength
error('Error:The xInArr %d and the yInArr %d are not equal in length.',xInLength,yInLength);
end
% ************************************************************
% 计算查询节点在插值函数上的函数值
% ************************************************************
% 只一个插值节点
if xInLength==1
yOutArr=yInArr;
return;
end
% 有多个插值节点
xQueryLength=length(xQuery);
yOutArr=zeros(1,xQueryLength);
pos=1;
for i=1:1:xQueryLength
for j=pos:1:xInLength-1
% 查询节点如果超出插值节点范围则报错
if xQuery(i)>xInArr(xInLength)||xQuery(i)<xInArr(1)
error('Error:The num of the xQuery should gainer than %d and lesser than %d.',xInArr(1),xInArr(xInLength));
end
% 定位查询数组在插值节点的位置
if (xQuery(i)>xInArr(j)&&xQuery(i)<xInArr(j+1))||(xQuery(i)==xInArr(j))
pos=j;
break;
end
end
% 计算查询节点在插值函数上的函数值
yOutArr(i)=polyval(func_2dotsCubicHermitInterpolation(xInArr(pos:1:pos+1),yInArr(pos:1:pos+1),dyInArr(pos),dyInArr(pos+1)),xQuery(i));
end
end
3.测试程序
% 测试分段低次插值函数——func_PiecewiseLinearInterpolatio、func_PiecewiseCubeHermitInterpolation
clc;
clear;
close all;
x=-0.2*pi:0.2*pi:10*pi;
y=sin(x);
dy=diff(y,1);
xq=0:0.01*pi:10*pi;
yq1=func_PiecewiseLinearInterpolation(x(2:1:length(x)),y(2:1:length(x)),xq);
yq2=func_PiecewiseCubeHermitInterpolation(x(2:1:length(x)),y(2:1:length(x)),dy,xq);
plot(x,y,'k*',xq,yq1,'r-',xq,yq2,'b-');
总结
以上在本文中对分段低次插值的两种常见方法的推导进行了简单介绍，并提供了MATLAB的实现代码与测试用例。
参考文献
[1]李庆扬,王能超,易大义. 数值分析.第5版[M]. 清华大学出版社, 2008.
440
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
