目录
1 插值与拟合定义
(1)插值:求过已知有限个数据点的近似函数
(2)拟合:已知有限个数据点,求近似函数,不要求过已知数据点,只要求在某种意义
下它在这些点上的总偏差最小。
(3)数据插值与曲线拟合比较
数据插值要求逼近函数经过样本点,没有统一的逼近函数;曲线拟合只要求总体误差最小,有确定的逼近函数形式。
数据插值一般用于样本区间内的插值计算,适合有精确数据的样本;
曲线拟合不仅估算区间内其他点的函数值,还要预测;适合统计数据的样本。
2 插值方法
2.1 拉格朗日多项式插值(Lagrange插值)
2.1.1 插值多项式(代数插值)
n次多项式有n+1个待定系数,由插值条件式恰好给出n+1个方程。代入点处的函数值,得到方程组Ax=b,det(A)为范德蒙行列式,在节点互不相同的情况下,由克莱姆法则,此时的插值多项式系数是存在唯一的。
插值多项式与被插函数之间的差:
2.1.2 拉格朗日插值多项式
随着等分值的值增大,拟合精度越来越小,但是在函数的两边值会出现巨变,也就是所谓的龙格现象。
2.1.3 Matlab实现Lagrange插值
%{
n个节点数据以数组 x0, y0输入
m个插值点以数组 x 输入
输出数组 y 为m 个插值。
参数输入:
x0,y0是原有的数据
x为插值点(范围)
参数输出
y为x的插值输出
%}
function y=lagrange(x0,y0,x);
n=length(x0);
m=length(x);
for i=1:m
z=x(i);
sum=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)); %n-1次多项式
end
end
sum=p*y0(k)+sum;
end
y(i)=sum;
end
2.1.4 实操
% 拉格朗日插值法
% 原始数据
x0 = 0:0.1:1;
y0 = [-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2];
hold on;
grid on;
plot(x0,y0,'b*')
% 插值点
x = 0:0.01:1;
% 拉格朗日多项式插值
y = lagrange(x0,y0,x);
plot(x,y,'r--');
legend('原始数据','拉格朗日插值','Location','North');
运行结果
2.2 牛顿插值(Newton)
2.2.1 差商定义
2.2.2 牛顿插值公式
2.2.3 matlab代码
% newton.m
% 求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序
% 输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量,Y是纵坐标向量,
% x是以向量形式输入的m个插值点,M在[a,b]上满足|f~(n+1)(x)|≤M
% 注:f~(n+1)(x)表示f(x)的n+1阶导数
% 输出的量:向量y是向量x处的插值,误差限R,n次牛顿插值多项式L及其系数向量C,
% 差商的矩阵A
function[y,R,A,C,L] = newton(X,Y,x,M)
n = length(X);
m = length(x);
for t = 1 : m
z = x(t);
A = zeros(n,n);
A(:,1) = Y';
s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;
for j = 2 : n
for i = j : n
A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));
end
q1 = abs(q1*(z-X(j-1)));
c1 = c1 * j;
end
C = A(n, n); q1 = abs(q1*(z-X(n)));
for k = (n-1):-1:1
C = conv(C, poly(X(k)));
d = length(C);
C(d) = C(d) + A(k,k);%在最后一维,也就是常数项加上新的差商
end
y(t) = polyval(C,z);
R(t) = M * q1 / c1;
end
L = poly2sym(C);
2.2.4 实操
X = [0 pi/6 pi/4 pi/3 pi/2];
Y = [0 0.5 0.7071 0.8660 1];
x = linspace(0,pi,50);
hold on;
grid on;
M = 1;
[y,R,A,C,L] = newton(X, Y, x, M);
plot(X,Y,'r*');
plot(x,y,'b--');
legend('原始数据','牛顿插值');
运行结果
2.3 分段线性插值
将每两个相邻的节点用直线连起来,如此形成的折线就是分段线性插值函数。
高次插值多项式会发生震荡,促使人们转而寻求简单的低次多项式插值。
2.3.1 matlab代码
Matlab里有现成的一维插值函数interp1.
2.4 埃尔米特(Hermite)插值
2.4.1 定义
如果要求插值函数不仅在节点处与函数同值,而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值,这就是Hermite插值问题。
2.4.2 matlab代码
%{
设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
2.5 样条插值
许多工程技术中提出的计算问题对插值函数的光滑性有较高要求,如飞机的机翼外形,内燃机的进、排气门的凸轮曲线,都要求曲线具有较高的光滑程度,不仅要连续,而且要有连续的曲率,这就导致了样条插值的产生。
2.5.1 matlab代码
% Lagrange插值、分段线性插值、三次样条插值
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;
% 调用前面编写的Lagrange插值函数
y1 = lagrange(x0,y0,x);
% 分段线性插值
y2 = interp1(x0,y0,x);
% 三次样条插值
y3 = interp1(x0,y0,x,'spline');
% 边界为Largrange的三次样条插值
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,3,1), plot(x0,y0,'+',x,y1), title('Lagrange')
subplot(2,3,2), plot(x0,y0,'+',x,y2), title('Piecewise linear')
subplot(2,3,3), plot(x0,y0,'+',x,y3), title('Spline1')
subplot(2,3,4), plot(x0,y0,'+',x,y4), title('Spline2')
subplot(2,3,5), plot(x0,y0,'+',x,y5), title('Spline3')
dyx0=ppval(fnder(pp1),x0(1)) %求x=0处的导数
ytemp=y3(131:151);
index=find(ytemp==min(ytemp));
xymin=[x(130+index),ytemp(index)]
2.6 二维插值
2.6.1 前言
前面讲述的都是一维插值,即节点为一维变量,插值函数是一元函数(曲线)。若节点是二维的,插值函数就是二元函数,即曲面。如在某区域测量了若干点(节点)的高程(节点值),为了画出较精确的等高线图,就要先插入更多的点(插值点),计算这些点的高程(插值)。
调用格式:Z1=interp2(X,Y,Z,X1,Y1,method)
其中,X、Y是两个向量,表示两个参数的采样点,Z是采样点对应的函数值。X1、Y1是两个标量或向量,表示要插值的点。
2.6.2 matlab代码
%{
例2 在一丘陵地带测量高程,x 和 y 方向每隔100米测一个点,得高程如2表,试插
值一曲面,确定合适的模型,并由此找出最高点和该点的高程。
%}
% 二维插值
clear,clc
x=100:100:500;
y=100:100:400;
z=[636 697 624 478 450
698 712 630 478 420
680 674 598 412 400
662 626 552 334 310];
% hold on;
xi = 100:10:500;yi=100:10:400
% 三次样条插值注意得到的cz1需要到转一下
pp = csape({x,y},z')
cz1 = fnval(pp,{xi,yi})
% 二维插值
cz2 = interp2(x,y,z,xi,yi','spline')
surfl(x,y,z-400);
hold on;
surfl(xi,yi,cz1');
surfl(xi,yi,cz2+400);
[i,j] = find(cz1==max(max(cz1)))
x = xi(i),y=yi(j),zmax=cz1(i,j)
2.6.3 插值节点为散乱节点
已知n 个节点:(xi , yi ,zi)(i = 1,2,L,n)
,求点(x, y)处的插值 z 。对上述问题,Matlab 中提供了插值函数 griddata
,其格式为: ZI = GRIDDATA(X,Y,Z,XI,YI,method)
method的选择有:linear(基于三角形的线性插值)、cubic(基于三角形的三次插值)、nearest(最邻近插值法)
其中 X、Y、Z 均为 n 维向量,指明所给数据点的横坐标、纵坐标和竖坐标。向量 XI、
YI 是给定的网格点的横坐标和纵坐标,返回值 ZI 为网格(XI,YI)处的函数值。XI
与 YI 应是方向不同的向量,即一个是行向量,另一个是列向量。
2.6.4 代码示例
% 插值节点为散乱节点
x=[129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5];
y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5];
z=-[4 8 6 8 6 8 8 9 9 8 8 9 4 9];
xi=75:1:200;
yi=-50:1:150;
zi=griddata(x,y,z,xi,yi','cubic')
subplot(1,2,1), plot(x,y,'*')
subplot(1,2,2), mesh(xi,yi,zi)
2.7 B样条函数插值方法
实际中的许多问题,往往是既要求近似函数(曲线或曲面)有足够的光滑性,又要求与实际函数相同的凹凸性,一般插值函数和样条函数都不具有这种性质。如果对于一个特殊函数进行磨光处理生成磨光函数(多项式),则用磨光函数构造出样条函数作为插值函数,既有足够的光滑性,而且也具有较好的保凹凸性,因此磨光函数在一维插值(曲线)和二维插值(曲面)问题中有着广泛的应用。
由积分理论可知,对于可积函数通过积分会提高函数的光滑度,因此,我们可以利用积分方法对函数进行磨光处理。
3 拟合
已知一组(二维)数据,即平面上的n 个点(xi , yi) , i = 1,2,L,n , xi 互不相同,寻求一个函数(曲线)y = f (x) ,使 f (x)在某种准则下与所有数据点最为接近,即曲线拟合得最好。
3.1 系数a的确定
3.2 函数r(x)的选取
面对一组数据用线性最小二乘法作曲线拟合时,首要的选取合适的r(x)。如果通过机理分析,了解到y与x之间的函数关系,则r(x)容易确定。若无法知道y与x之间的关系,则需要将数据作图,直观地判断用什么样的曲线进行拟合。
对于指数曲线,拟合前需作变量代换,化为对a1,a2 的线性函数。已知一组数据,用什么样的曲线拟合最好,可以在直观判断的基础上,选几种曲线分别拟合,然后比较,看哪条曲线的最小二乘指标 J 最加粗样式小。
3.3 解方程法
3.3.1 方法解析
3.3.2 例题
3.3.3 matlab代码
% 解方程拟合
x=[19 25 31 38 44]';
y=[19.0 32.3 49.0 73.3 97.8]';
r=[ones(5,1),x.^2];
ab=r\y
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')
3.4 多项式拟合方法
函数polyfit()功能是求得最小二乘拟合多项式系数。其调用格式为:
其中输入参数 x0
,y0
为要拟合的数据,m
为拟合多项式的次数,输出参数 p
为拟合多项
式。根据样本数据X和Y,产生一个m次多项式P及其在采样点误差数据S,m为几次拟合,mu
是一个二元向量,mu(1)
是mean(X)
,而mu(2)
是std(X)。
3.4.1 例题
3.4.2 matlab代码
% 多项式拟合
% 原始数据
x0=[1990 1991 1992 1993 1994 1995 1996];
y0=[70 122 144 152 174 196 202];
% 画散点图
plot(x0,y0,'*')
hold on;
a=polyfit(x0,y0,1)
x = 1990:0.1:2000;
y = polyval(a,x);
plot(x,y,'r--');
y97=polyval(a,1997)
y98=polyval(a,1998)
a =
1.0e+04 *
0.0021 -4.0705
y97 =
233.4286
y98 =
253.9286
3.5 最小二乘优化
在无约束最优化问题中,有些重要的特殊情形,比如目标函数由若干个函数的平
方和构成。
(1)lsqin函数
% 用 lsqlin 命令求解例 4。
x=[19 25 31 38 44]';
y=[19.0 32.3 49.0 73.3 97.8]';
r=[ones(5,1),x.^2];
ab=lsqlin(r,y)
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')
(2)lsqcurvefit函数
function f=fun1(x,tdata);
f=x(1)+x(2)*exp(-0.02*x(3)*tdata); %其中 x(1)=a,x(2)=b,x(3)=k
td=100:100:1000;
cd=[4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59];
x0=[0.2 0.05 0.05];
x=lsqcurvefit(@fun1,x0,td,cd)
(3)lsqnonlin函数
function f=fun2(x);
td=100:100:1000;
cd=[4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59];
f=x(1)+x(2)*exp(-0.02*x(3)*td)-cd;
x0=[0.2 0.05 0.05]; %初始值是任意取的
x=lsqnonlin(@fun2,x0)
(4)lsqnonneg函数
c=[0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170];
d=[0.8587;0.1781;0.0747;0.8405];
x=lsqnonneg(c,d)