定义
在实际问题中,一般通过实际观测得出一个函数
y=f(x)
,并知道有限个点。
yi=f(xi),i=0,1,...,n
当需要知道 x0,x1,...,xn 之间的点x的函数值,那么就需要进行插值,常用一些较简单的,满足条件的函数 g(x) 来代替 f(x) ,这就是插值法。要经过已知的数据点。
拟合也是已知有限个数据点,求近似函数,不要求过已知数据点,只要求在某种意义下它在这些点上的总偏差最小,即最佳地拟合数据。
分段线性插值
将每个相邻点用直线连起来。如此形成的折线就是分段线性插值。
样条插值
使得插值函数更加光滑。将具有一定光滑性的分段多项式称为样条函数。
MATLAB实现插值
拉格朗日插值法:
%% 设n个节点数据以数组x0, y0输入(注意Matlat 的数组下标从1 开始),m 个插值点以数组x 输入,输出数组y 为m 个插值
function y=lagrange(x0,y0,x);
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
埃尔米特(Hermite)插值
%% 设n个节点的数据以数组x0(已知点的横坐标),y0(函数值),y1(导数值)输入(注意Matlat 的数组下标从1 开始),m 个插值点以数组x 输入,输出数组y 为m个插值。
function y=hermite(x0,y0,y1,x);
n=length(x0);m=length(x);
for k=1:m
yy=0.0;
for i=1:n
h=1.0;
a=0.0;
for j=1:n
if j~=i
h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;
a=1/(x0(i)-x0(j))+a;
end
end
yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));
end
y(k)=yy;
end
一维插值函数
语法:y = interp1(x0, y0, x, ‘method’)
method指定插值的方式,默认为线性插值。有以下几种插值方法:
‘nearest’ 最近项插值
‘linear’ 线性插值
‘spline’ 立方样条插值
‘cubic’ 立方插值
所有插值方法要求
x0
是单调的。
当
x0
位等距时可以用快速插值法,使用快速插值法的格式为 ‘*nearest’、 ‘*linear’、 ‘*spline’、 ‘*cubic’。
三次样条插值
一维插值函数:
利用样条函数进行插值,即取插值函数为样条函数,称为样条插值。
y = interp1(x0, y0, x, ‘spline’);
y = spline(x0, y0, x);
pp = csape(x0, y0, conds);
pp = csape(x0, y0, conds, valconds); y = ppval(pp, x);
其中
x0,y0
为已知数据点,x是插值点,y是插值点的函数值。
举例:12小时内,每小时测量一次室外温度,并进行样条插值
几种一维插值方式比较:
clc,clear
x0=[0 3 5 7 9 11 12 13 14 15];
y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x=0:0.1:15;
y1=lagrange(x0,y0,x); % 调用前面编写的Lagrange插值函数
y2=interp1(x0,y0,x);
y3=interp1(x0,y0,x,'spline');
pp1=csape(x0,y0); y4=ppval(pp1,x);
pp2=csape(x0,y0,'second'); y5=ppval(pp2,x);
fprintf('比较一下不同插值方法和边界条件的结果:\n')
fprintf('x y1 y2 y3 y4 y5\n')
xianshi=[x',y1',y2',y3',y4',y5'];
fprintf('%f\t%f\t%f\t%f\t%f\t%f\n',xianshi')
subplot(2,2,1), plot(x0,y0,'+',x,y1), title('Lagrange')
subplot(2,2,2), plot(x0,y0,'+',x,y2), title('Piecewise linear')
subplot(2,2,3), plot(x0,y0,'+',x,y3), title('Spline1')
subplot(2,2,4), plot(x0,y0,'+',x,y4), title('Spline2')
dyx0=ppval(fnder(pp1),x0(1)) %求x=0处的导数
ytemp=y3(131:151);
index=find(ytemp==min(ytemp));
xymin=[x(130+index),ytemp(index)]
效果:
二维插值函数:
z = interp2(x0, y0, z0, x, y, ‘method’)
其中x0,y0分别为想m维和n维向量,表示节点,z0为n*m维矩阵,表示节点值,x,y为一维数组,表示插值点,x与y应是方向不同的向量,即一个是行向量,一个是列向量,z为矩阵。
如果使用三次样条插值,可以使用命令
pp = csape({x0, y0}, z0, conds, valconds); z = fnval(pp, {x, y});
举例:
MATLAB实现拟合
线性最小二乘法
问题:已知一组数据,即平面上n个点,i=1,2,…,n,xi互不相同,寻求一个函数y = f(x),使f(x)在某种准则下与所有数据点最为接近,即曲线拟合得最好。
令
f(x)=a1r1(x)+a2r2(x)+...+amrm(x)
,其中
rk(x)
是事先选定的一组线性无关的函数,
ak
是待定系数。
拟合准则是使
yi
,
i=1,2,...,n
,与
f(xi)
的距离
δi
的平方和最小,称为最小二乘准则。
多项式拟合方法
a = ployfit(x0, y0, m),其中输入参数x0, y0为要拟合的数据,m为拟合多项式的次数,输出参数a为拟合多项式
y=a(1)xm+...+a(m)x+a(m+1)
的系数向量a = [a(1),…,a(m),a(m+1)]。
多项式在x处的值y可用下面的函数计算y = polyval(a,x)
举例:
最小二乘优化
用于求解最小二乘优化问题的函数有lsqlin、lsqcurvefi、lsqnonlin、lsqnonneg
x = lsqcurvefit(‘fun’, x0, xdata, ydata, options);
x = lsqnonlin(‘fun’, x0, options);