《计算方法学习指导》实验(3)

第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}$$

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.8108x^2+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.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.0029e^{0.0016t}$$