函数插值是差值函数p(x)与被插函数f(x)在节点处函数值相同,即p( )=f( )
不都严格地等于零。但是,为了使近似曲线能尽量反应所给数据点的变化趋势,要求| |按某种度量标准最小。
即
为最小。这种要求误差平方和最小的拟合称为曲线拟合的最小二乘法。
根据线性最小二乘拟合理论,我们得知关于系数矩阵A的解法为A=R\Y。
例题 假设测出了一组,由下面的表格给出,且已知函数原型为
y(x)=c1+c2*e^(-3*x)+c3*cos(-2*x)*exp(-4*x)+c4*x^2
x |
0 |
0.2 |
0.4 |
0.7 |
0.9 |
0.92 |
0.99 |
1.2 |
1.4 |
1.48 |
1.5 |
y |
2.88 |
2.2576 |
1.9683 |
1.9258 |
2.0862 |
2.109 |
2.1979 |
2.5409 |
2.9627 |
3.155 |
3.2052 |
试用已知数据求出待定系数的值。
x=[0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5]';
y=[2.88;2.2576;1.9683;1.9258;2.0862;2.109;
A=[ones(size(x)) exp(-3*x),cos(-2*x).*exp(-4*x)x.^2];
c=A\y;
c'
运行结果为
ans =
x1=[0:0.01:1.5]';
A1=[ones(size(x1))exp(-3*x1),cos(-2*x1).*exp(-4*x1) x1.^2];
y1=A1*c;
plot(x1,y1,x,y,'o')
y(x)= 0.8700-0.6797*e^(-3*x)+2.3397*cos(-2*x)*exp(-4*x)+ 1.2200*x^2
产生的,由上图可见拟合效果较好。
在Matlab的线性最小二乘拟合中,用得较多的是多项式拟合,其命令是
A=polyfit(x,y,m)
其中 表示函数中的自变量矩阵, 表示因变量矩阵,是输出的系数矩阵,即多项式的系数。
多项式在自变量x处的函数值y可用以下命令计算:
y=polyval(A,x)
例题 对下面一组数据作二次多项式拟合,即要求出二次多项式 中的 ,使最小。
x |
0 |
0.1 |
0.2 |
0.3 |
0.4 |
0.5 |
0.6 |
0.7 |
0.8 |
0.9 |
1 |
y |
-0.447 |
1.978 |
3.28 |
6.16 |
7.08 |
7.34 |
7.66 |
9.56 |
9.48 |
9.30 |
11.2 |
在Matlab中输入以下命令
x=0:.1:1;
y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.489.30 11.2];
a=polyfit(x,y,2)
运行结果为
a =
f=vpa(poly2sym(a),5)
%vpa(polyval2sym(),n)只适用于关于多项式函数的拟合。因为此函数对于自变量统一规定为“x”,将由polyfit()所得出的系数按自变量幂次升降放在相应的位置。
运行结果为
f =
下面画出由拟合得到的曲线及已知的数据散点图
y1=polyval(a,x);
plot(x,y,'o',x,y1)
(1)lsqcurvefit( )
lsqcurvefit( )是非线性最小二乘拟合函数,其本质上是求解
最优化问题。其使用格式为
x=lsqcurvefit(fun,x0,xdata,ydata)
其中,fun是要拟合的非线性函数,x0是初始参数,xdata,ydata是拟合点的数据,该函数最终返回系数矩阵。
例题 假设已知
并已知该函数满足原型为 ,其中 为待定系数。
在Matlab中输入以下命令
x=0:.1:10;
y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);
f=inline('a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x)','a','x');
%建立函数原型,则可以根据他来进行下面的求取系数的计算
[a,res]=lsqcurvefit(f,[1,1,1,1,1],x,y);
a',res
运行结果为
ans =
res =
所求得的系数与原式中的系数相近。
如果向进一步提高精度,则需修改最优化的选项,函数的调用格式也将随之改变。
在Matlab中输入以下命令
ff=optimset;ff.TolFun=1e-20;ff.TolX=1e-15;%修改精度,即是修改其限制条件
[a,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff);%两个空矩阵表示系数向量的上下限
a',res
运行结果为
ans =
res =
9.5035e-021
下面绘图
x1=0:0.01:10;
y1=f(a,x1);
plot(x1,y1,x,y,'o')
(2)lsqnonlin( )
lsqnonlin( )函数是另一种求最小二乘拟合的函数,其本质上是求解最优化问题
最优化解。它的应用格式为
x=lsqnonlin(fun,x0)
其中,fun的定义与lsqnonlin( )函数中fun的定义有差别,x0仍为初始参数向量,将输出的系数结果放在变量 中。
说明lsqnonlin()函数的使用方法。
首先编写目标函数 (curve_fun.m)
function y=curve_fun(p)%非线性最小二乘拟合函数
x=[0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0.560.56 1.10 1.10];
y=[76 47 97 107 123 139 159 152 191 201 207200];
y=p(1)*x./(p(2)+x)-y;
再用lsqnonlin()函数求解,输入
[p,resnorm,residual]=lsqnonlin(@curve_fun,[200,0.1])
运行结果为
p =
resnorm =
residual =
上面的两种方法都可以作非线性最小二乘曲线拟合
(3)非线性函数的线性化
在进行非线性拟合时,以往由于计算机和相关软件水平有限,常常先把非线性的曲线拟合线性化,然后再进行拟合。下面比较一下先线性化再进行拟合和直接进行非线性拟合的差异。
例题 已知数据
t |
0.25 |
0.5 |
1 |
1.5 |
2 |
3 |
4 |
6 |
8 |
c |
19.21 |
18.15 |
15.36 |
14.10 |
12.89 |
9.32 |
7.45 |
5.24 |
3.01 |
满足曲线 通过数据拟合求出参数和 。
方法一:先将非线性函数转化为线性函数
编写Matlab程序如下
d=300;
t=[0.25 0.5 1 1.5 2 3 4 6 8];
c=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.243.01];
运行结果为
a =
k =
0.2347
v =
由此也可以求出相关系数。
方法二:应用非线性拟合直接求解系数
建立m文件:
function f=curvefun3(x,tdata)
d=300
f=(x(1)\d)*exp(-x(2)*tdata)%x(1)=v;x(2)=k
运行程序
tdata=[0.25 0.5 1 1.5 2 3 4 6 8];
cdata=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.243.01];
x0=[10 0.5];
x=lsqcurvefit('curvefun3',x0,tdata,cdata)
运行结果为
x =
f=curvefun3(x,tdata);
plot(tdata,cdata,'o',tdata,f)
我们发现两种求法求出的系数很接近。
(三)线性拟合和非线性拟合区别与联系
在许多实际问题中,变量之间内在的关系并不想前面说的那样简单。呈线性关系,但有些非线性拟合曲线可以通过适当的变量替换转化为线性曲线,从而用线性拟合进行处理。对于一个实际的曲线拟合问题,一般先根据观测值在直角坐标平面上描出散点图,看一看散点的分布同哪类曲线图形接近,让后选用相接近的曲线拟合方程,再通过适当的变量替换转化为线性拟合问题,按线性拟合解出后再还原为原变量所表示的曲线拟合方程。
表1.1
线性拟合方程 |
变量变换 |
变换后线性拟合方程 |
Y= |
, |
|
Y=a |
|
Y=a |
Y= |
|
|
Y= |
|
|
Y= |
|
|
例题
测出一组实际数据 见下表 是对其进行函数拟合。
X |
1.1052 |
1.2214 |
1.3499 |
1.4918 |
1.6478 |
3.6693 |
Y |
0.6795 |
0.6006 |
0.5309 |
0.4693 |
0.4148 |
0.1546 |
X |
1.8221 |
2.0138 |
2.2255 |
2.4596 |
2.7183 |
|
Y |
0.3666 |
0.3241 |
0.2864 |
0.2532 |
0.2238 |
|
>>x=[1.1052,1.2214,1.3499,1.4918,1.6478,1.8221,2.0138,2.2255,2.4596,2.7183,3.6693];
>>y=[0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,0.2864,0.2532,0.2238,0.1546];
>>plot(x,y,x,y,'*')
见下图
由上图可以看出是非线性曲线y=a
Y= |
, |
|
分别对x,y进行对数变换
>>x1=log(x);y1=log(y);plot(x1,y1)
用线性函数拟合的方法可以求出现行参数,使得ln y=aln x+b,即y= ,求解系数a,b及 。
>>A=[x1',ones(size(x1'))];c=[A\y1']
c =
>> exp(c(2))
ans =
拟合函数为y(x)=0.768 703 388 199 24