第3章 曲线拟合的最小二乘法
3.1 实验目的
结合最小二乘原理用MATLAB编写相应的曲线拟合程序。
3.2 MATLAB命令
函数 | 含义 |
---|---|
polyfit | 多项式拟合 |
lsqcurvefit | 曲线拟合 |
lsqnonlin | 最小二乘法 |
csaps | 样条拟合 |
MATLAB解决曲线拟合问题的三种方法:
调用函数,以表格形式输出拟合函数;调用函数,以公式形式输出拟合函数;由用户图形界面进行曲线拟合。
3.3 实验3例题:曲线拟合
例1:对以下数据进行1次、2次、3次多项式拟合,并求出其总体误差。
x | 0 | 0.3 | 0.6 | 0.9 | 1.2 | 1.5 | 1.8 | 2.1 | 2.4 | 2.7 | 3 |
---|---|---|---|---|---|---|---|---|---|---|---|
y | 0 | 0.29 | 0.56 | 0.78 | 0.93 | 0.99 | 0.97 | 0.86 | 0.67 | 0.42 | 0.14 |
用以下多项式进行拟合:
f
(
x
)
=
a
1
x
m
+
a
2
x
m
−
1
+
.
.
.
+
a
m
x
+
a
m
+
1
f(x)=a_1x^m+a_2x^{m-1}+...+a_mx+a_{m+1}
f(x)=a1xm+a2xm−1+...+amx+am+1
clc;
clear all;
xi=0:0.3:3;
yi=[0 0.29 0.56 0.78 0.93 0.99 0.97 0.86 0.67 0.42 0.14];
plot(xi,yi,'.K','MarkerSize',15); %画出数据点的散点图
hold on;
x=0:0.01:3;
p1=polyfit(xi,yi,1);%一次多项式拟合
y1=polyval(p1,x);
yy1=polyval(p1,xi);
err1=norm(yy1-yi,2);
p2=polyfit(xi,yi,2);%二次多项式拟合
y2=polyval(p2,x);
yy2=polyval(p2,xi);
err2=norm(yy2-yi,2);
p3=polyfit(xi,yi,3);%二次多项式拟合
y3=polyval(p3,x);
yy3=polyval(p3,xi);
err3=norm(yy3-yi,2);
plot(x,y1,'Y');
hold on;
plot(x,y2,'R');
hold on;
plot(x,y3,'B');
hold on;
err1
err2
err3
结果:
err1 =
1.0847
err2 =
0.0686
err3 =
0.0652
polyfit()和polyval()函数
p=polyfit(x,y,n)
[p,s]= polyfit(x,y,n)
说明:x,y为数据点,n为多项式阶数,返回p为幂次从高到低的多项式系数向量p。x必须是单调的。矩阵s用于生成预测值的误差估计。
多项式曲线求值函数:polyval( )
调用格式:
y=polyval(p,x)
[y,DELTA]=polyval(p,x,s)
说明:y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。
[y,DELTA]=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计Y DELTA。它假设polyfit函数数据输入的误差是独立正态的,并且方差为常数。则Y DELTA将至少包含50%的预测值。
norm()函数
格式:n=norm(A,p) 功能:norm函数可计算几种不同类型的矩阵范数,根据p的不同可得到不同的范数。
例2:对以下数据做二次多项式拟合。
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 |
法一:解超定方程
x=0:0.1:1;
y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];
R=[(x.^2)',x',ones(11,1)];
A=R\y'
结果:
A =
-9.8108
20.1293
-0.0317
因此得到的拟合二次多项式是:
f
(
x
)
=
−
9.8108
x
2
+
20.1293
x
−
0.0317
f(x)=-9.8108x^2+20.1293x-0.0317
f(x)=−9.8108x2+20.1293x−0.0317
法二:使用多项式拟合命令
x=0:0.1:1;
y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];
A=polyfit(x,y,2);
z=polyval(A,x);
plot(x,y,'.k','MarkerSize',15);
hold on;
plot(x,z,'r')
结果:
例3:用下面给出的数据拟合给出的y(t)中的参数a,b,k。
y
(
t
)
=
a
+
b
e
0.02
k
t
y(t)=a+be^{0.02kt}
y(t)=a+be0.02kt
数据:
t | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 |
---|---|---|---|---|---|---|---|---|---|---|
y | 4.54 | 4.99 | 5.35 | 5.65 | 5.90 | 6.10 | 6.26 | 6.39 | 6.50 | 6.59 |
法一:使用lsqcurvefit
curvefun1.m:
tdata=100:100:1000;
ydata=[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('curvefun1',x0,tdata,ydata)
命令及结果:
>> tdata=100:100:1000;
ydata=1e-3*[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('curvefun1',x0,tdata,ydata);
f=curvefun1(x,tdata);
x =
0.0069 -0.0029 0.0809
f =
Columns 1 through 7
0.0044 0.0048 0.0051 0.0054 0.0056 0.0058 0.0060
Columns 8 through 10
0.0061 0.0063 0.0064
拟合函数为:
y
(
t
)
=
0.0069
−
0.0029
e
0.0016
t
y(t)=0.0069-0.0029e^{0.0016t}
y(t)=0.0069−0.0029e0.0016t