插值与拟合的Matlab实现

王正盛编写

在科技工程中，除了要进行一定的理论分析外，通过实验、观测数据，做分析、处理也是必不可少的一种途径。由于实验测定实际系统的数据具有一定的代表性，因此在处理时必须充分利用这些信息；又由于测定过程中不可避免会产生误差，故在分析经验公式时又必须考虑这些误差的影响。两者相互制约。据此合理建立实际系统数学模型的方法成为数值逼近法。

一、插值法

1、数学原理

工程实践和科学实验中，常常需要从一组实验观测数据中，求自变量与因变量的一个近似的函数关系式。

例如：观测行星的运动，只能得到某时刻所对应的行星位置(用经纬度表示)，想知道任何时刻的行星位置。

例如：大气压测定问题；导弹发射问题；程序控制铣床加工精密工件问题；飞机船舶制造问题等等。都属于此类问题。因为考虑到代数多项式既简单又便于计算，所以人们就用代数多项式近似地表示满足个点的函数关系式——插值法建模。

(1)计算方法课程中学习了两种多项式插值：Lagrange插值和Newton均差插值：

已知n+1个数据点：

n次Lagrange插值公式：

特别地，当n=1时，————线性插值

当n=2时，

———————抛物线插值或二次插值

Newton均差插值公式：，其中是k阶均差，可由均差表方便计算得到。

Lagrange插值和Newton均差插值本质上是一样的，只是形式不同而已，因为插值多项式是唯一的。

(2)Runge现象和分段低次插值：

如在[-5，5]上各阶导数存在，但在此区间取n个节点构造的Lagrange插值多项式在区间并非都收敛，而且分散得很厉害。(matlab\bin\ Lagrange.m是自己编写的M文件)

[例]取n=-10

y=1./(1+x.^2);

x0=[-5:0.1:5];

y0=lagrange(x,y,x0);

y1=1./(1+x0.^2);

plot(x0,y0)

hold on

plot(x0,y1,'b:')

legend('插值曲线','原数据曲线')

因此插值多项式一般不要超过四次为宜。为避免Runge现象提出分段插值，所谓分段插值就是首先把插值点分开，在每一段上用低次插值，再连接起来。

如分段线性插值、分段二次插值、分段三次插值等等。

(3)三次样条插值：

三次样条插值算法的插值精度较高，所构造的曲线比较光滑(即在节点处连续可导)。因此，在许多工程设计或制造业中，例如飞机、导弹、高速机车、万吨级轮船等外形设计常常利用本算法进行插值计算。

三次样条插值函数的定义：

设已知n+1个数据，如果函数满足如下性质：

(1)在每个子区间上，是次数不超过三次的多项式；

(2)

(3)在插值区间上二次连续可微，记为

则称为三次样条插值函数，简称三次样条。

三次样条函数在每个子区间上可由4个系数唯一确定，因此，在上有4n个待定系数。由于，

故有

这里给出了3n-3个条件，加上插值条件(2)，共有4n-2个条件。为了确定，通常还需要补充两个边界条件。

常用的边界条件有三类：

第一种边界条件：给定两边界节点处的一阶导数，并要求满足：

第二种边界条件：给定两边界节点处的二阶导数，并要求满足：

特别地，若，则所得的样条称为自然样条(MATLAB中即是自然样条)。

第三种边界条件：被插函数是以为周期的周期函数，要求满足：

定理：三次样条插值问题的解存在且唯一。(证明略)

求三次样条函数有多种方法，这里介绍一种著名的三弯矩法。

1、表达式：设三次样条函数为，且，在上，由定义中的(2)及处的二阶导数为，则有

(1)

(2)

由于是三次多项式，所以是线性函数，于是线性插值，得

令，得

(3)

对(3)式积分两次，得

(4)

其中为积分常数，利用(1)，有

解得

代入(4)并化简，得

——————————————————(5)这就是系数用二阶导数表示的三次样条的表达式。因此，只要确定，即得三次样条函数。

2、结点的关系式

为了求，利用的一阶导数的连续性：；

将(5)式对求导数，得

类似地，可得

由，得

上式两边同乘以，得

令

则得到方程组

(6)

称为结点的M关系式(力学上称为三弯矩方程组)

3、边界条件

(6)式中是n+1个未知数的n-1个方程的方程组，不能唯一的确定，要唯一确定，必须附加两个条件，可由实际情况在两端点处给出的条件，称为边界条件(端点条件)。常见两种边界条件：

(1)给定两端点的二阶导数，(7)

这时由结点的M关系式(6)可求出，特别地当时为自然样条。

(2)给出两端点的一阶导数(斜率)：，此时由一阶导数的连续性，可得到两个方程：

(8)

边界条件(7)和(8)可以写成统一形式

(9)

其中

当时，即为(7)式；当时，即为(8)式。

将(6)式与(9)式合在一起得到以下三对角方程组：

此方程组对角占优，用追赶法一定可解出。

4、解题步骤

第一步：确定边界条件；

第二步：建立结点的M关系式，求出各点的；

第三步：将所得的代回样条函数表达式(5)，即得在个子区间上的表达式。

例题：某飞机头部外形曲线数据如下：x 0 70 130 210 337 578 776 1012 1142 1462 1841

y 0 57 78 103 135 182 214 244 256 272 275

设给定端点条件：时，；时，

求插值点处的函数值及一阶、二阶导数值。

答案：x 100 250 400 500 800

y 69.440456 114.369173 148.323244 167.884110 217.485051

y' 0.319709 0.265181 0.204950 0.186896 0.143324

y" -0.004312 -0.000855 -0.000199 -0.000162 -0.000159

2、建模实例

例：丙烷导热系数的估计

通过实验得到如下数据T P K

68

87

106

140 9798.1

13324

9007.8

13355

9791.8

14277

9656.3

12463 0.0848

0.0897

0.0762

0.0807

0.0696

0.0753

0.0611

0.0651

表中T、P和K分别表示温度、压力和导热系数，并假设在这个范围内导热系数近似地随压力线性变化，求P=10130和T=99下的K。

对于此问题进行分析，不难发现这是一个两个自变量的插值问题，即要求同时对压力和温度进行内插。但是考虑到该问题假设，内插将分为二部分完成。首先对不同温度求出在P=10130下的一组导热系数，利用线性插值求出T=68和P=10130下的导热系数

类似可得到在T=87、106和140时的导热系数：T 68 87 106 140

K 0.0853 0.0774 0.0699 0.0618

其次对温度的插值，采用三次Lagrange插值多项式计算得到：

P=10130和T=99下的K=0.0725

3、在MATLAB中的实现

MATLAB提供一维、二维和N维插值函数interp1、interp2和interpn，以及三次样条插值函数spline等。

(1)一维数据插值(一元函数)

MATLAB提供的函数interp1可以根据已知数据表[ X,Y]，用各种不同的算法计算XI各点上的函数近似值。该函数有三种调用格式。

(1.1)YI=interp1(X,Y,XI)

根据数据表[X,Y]，用分段线性插值算法计算XI各点(或一个点)上的函数近似值YI。当Y是向量时，则对Y向量插值，得到结果YI是与XI同样大小的向量；当Y是矩阵时，则对Y的每一列向量插值，得到结果YI是一矩阵：它的列数与Y的列数相同，它的行数与向量XI的大小相同。

(1.2)YI=interp1(Y,XI)

格式调用方法与(1.1)相同，只是插值节点用其节点序号：X=1:N，这里的N=size(Y),即X的大小与Y一样。

(1.3) YI=interp1(X,Y,Xi,method)

函数调用与上述相同，只是需要指定输入参数method的具体的算法：

'linear'——分段线性插值，可缺省

'cubic'——分段三次多项式插值

'spline'——三次样条插值：即在每个分段(子区间)内构造一个三次多项式，使其插值函数满足插值条件外，还要求在各个节点处具有光滑的条件(导数存在)。

'nearest'——最邻近区域插值：即在就近插值节点的区域上的函数取值该节点之函数值，该插值函数为一个阶梯函数：

例1：已知数据表如下，试利用interp1的不同插值算法求xi=1,1.5,2,2.5,…,5各点的函数的近似值。X 1.0 2.0 3.0 4.0 5.0

Y 3.5 4.6 5.5 3.2 2.0

解：

hold off

xx=1:1:5;

yx=[3.5,4.6,5.5,3.2,2];

xxi=1:0.5:5;

f0=interp1(xx,yx,xxi);

f1=interp1(xx,yx,xxi,'linear');

f2=interp1(xx,yx,xxi,'cubic');

f3=interp1(xx,yx,xxi,'spline');

f4=interp1(xx,yx,xxi,'nearest');

f0,f1,f2,f3,f4

Columns 1 through 7

3.50004.05004.60005.05005.50004.35003.2000

Columns 8 through 9

2.60002.0000

f1 =

Columns 1 through 7

3.50004.05004.60005.05005.50004.35003.2000

Columns 8 through 9

2.60002.0000

f2 =

Columns 1 through 7

3.50004.07504.60005.26255.50004.48133.2000

Columns 8 through 9

2.46252.0000

f3 =

Columns 1 through 7

3.50003.77344.60005.37665.50004.59533.2000

Columns 8 through 9

2.07972.0000

f4 =

Columns 1 through 7

3.50004.60004.60005.50005.50003.20003.2000

Columns 8 through 9

2.00002.0000

legend('线性插值','三次插值','样条插值','最近区域插值')

例2：某观测站测得某日6：00~18：00之间每隔2小时的室内外温度列表如下，要求用三次样条插值分别求得该日室内外6：30~17：30之间每隔2小时各点的近似温度。时间h 6 8 10 12 14 16 18

室内温度t1 18.0 20.0 22.0 25.0 30.0 28.0 24.0

室外温度t2 15.0 19.0 24.0 28.0 34.0 32.0 30.0

解：设时间变量h为一行向量，温度t为一个2列矩阵，第一列存储室内温度，第二列存储室外温度。

t=[18 20 22 25 30 28 24;15 19 24 28 34 32 30]';

xi=6.5:2:17.5;

yi=interp1(h,t,xi,'spline')

18.502015.6553

20.498620.3355

22.519324.9089

26.377529.6383

30.205134.2568

26.817830.9594

例3：编制分段二次插值程序，即在每一个插值子区间上用抛物线插值。

MATLAB\BIN\lag2.m

以例1数据为例：

yx=[3.5 4.6 5.5 3.2 2];

xxi=1:0.5:5;

f0=interp1(xx,yx,xxi);

f1=interp1(xx,yx,xxi,'linear');

f2=interp1(xx,yx,xxi,'cubic');

f3=interp1(xx,yx,xxi,'spline');

f4=interp1(xx,yx,xxi,'nearest');

f5=lag2(xx,yx,xxi)

plot(xxi,f1,xxi,f2,xxi,f3,xxi,f4,xxi,f5,'b:')

Columns 1 through 7

3.50004.07504.60005.45005.50004.21253.2000

Columns 8 through 9

2.46252.0000

legend('线性插值','三次插值','样条插值','最近区域插值','二次插值')

(2)三次样条插值

三次样条插值算法的插值精度较高，所构造的曲线比较光滑。因此，在许多工程设计或制造业中，例如飞机、导弹、高速机车、万吨级轮船等外形设计常常利用本算法进行插值计算。它有两种调用格式：

(2.1)yi=spline(x,y,xi)

本格式根据数据表[x,y]用spline插值计算各节点处的函数近似值。这种格式与使用interp1函数做三次样条插值作用一样。

例：

yx=[3.5,4.6,5.5,3.2,2];

xxi=1:0.5:5;

f0=interp1(xx,yx,xxi,)

f1=interp1(xx,yx,xxi,'spline')

yi=spline(xx,yx,xxi)

Columns 1 through 7

3.50004.05004.60005.05005.50004.35003.2000

Columns 8 through 9

2.60002.0000

f1 =

Columns 1 through 7

3.50003.77344.60005.37665.50004.59533.2000

Columns 8 through 9

2.07972.0000

yi =

Columns 1 through 7

3.50003.77344.60005.37665.50004.59533.2000

Columns 8 through 9

2.07972.0000

后两者计算结果完全一样。

(2.2) pp=spline(x,y)

本格式根据数据表[x,y]用spline插值算法获取三次样条插值pp形式，即在pp中存储了各分段的三次样条插值多项式的系数和其他有关必要的信息。欲求得xi的插值结果yi，应利用yi=ppval(pp,xi)格式调用之。

例：

yx=[3.5,4.6,5.5,3.2,2];

xxi=1:0.5:5;

f0=interp1(xx,yx,xxi)

ps=spline(xx,yx);

yyi=ppval(ps,xxi)

Columns 1 through 7

3.50004.05004.60005.05005.50004.35003.2000

Columns 8 through 9

2.60002.0000

yyi =

Columns 1 through 7

3.50003.77344.60005.37665.50004.59533.2000

Columns 8 through 9

2.07972.0000

注意：两者计算结果不同的。

(3)二维数据插值(二元函数)

a.插值基点为网格节点

二维数据插值是构造一个二元插值函数去近似，即曲面插值。如测山高水深，画出更为精确的等高线图，就需要先插入更多的插值点。与interp1类似，指令interp2可以根据数据表[x,y,z]，用各种不同的算法计算[xi,yi]各点上的函数近似值zi。该函数有两种调用格式：

(3.1) zi=interp2(x,y,z,xi,yi)

本指令格式根据数据表[x,y,z]，用双线性插值算法计算坐标平面x-y上[xi,yi]各点(或一个点)的二元函数近似值zi。这里x可以是一行向量，它与z(矩阵)的各列向量相对应；y可以为一个列向量，它与z的各行向量相对应。对于[xi,yi]与zi间的对应关系，则和[x,y]与z的对应关系相同。

这里的双线性插值，是指按x,y的双向线性插值。

(3.2) i=interp2(z,xi,yi)

本指令格式与(3.1)的功能相同，只是插值节点用其节点序号：x=1:n,y=1:m，[m,n]=size(z).

(3.3) zi=interp2(x,y,z,xi,yi,method)

本指令格式的调用方法与上述指令格式相同，只是规定在格式中指定具体的算法。这里输入参数method有：

'linear'——双线性插值。默认值，即格式(3.1)

'cubic'——双三次插值

'nearest'——最邻近区域插值

例：某实验对一根长10 m的钢轨进行热源的温度传播测试。用x表示测试点0:2.5:10(m)，用h表示测试时间0:30:60(s)，用T表示测试所得各点的温度。

试用双线性插值在一分钟内每隔20s、钢轨每隔1m处的温度Ti。

解：

h=[0:30:60]'

T=[95 14 0 0 0 ;88 48 32 12 6;67 64 54 48 41]

mesh(x,h,T)

xi=[0:10];%在x,y轴取小的步长以平滑数据。

hi=[0:2:60]';

Ti=interp2(x,h,T,xi,hi)

02.50005.00007.500010.0000

h =

0

30

60

T =

9514000

884832126

6764544841

Ti =

Columns 1 through 7

95.000062.600030.200011.20005.600000

94.533363.226731.920013.44007.78672.1333 1.6000

94.066763.853333.640015.68009.97334.26673.2000

93.600064.480035.360017.920012.16006.40004.8000

93.133365.106737.080020.160014.34678.53336.4000

92.666765.733338.800022.400016.533310.66678.0000

92.200066.360040.520024.640018.720012.80009.6000

91.733366.986742.240026.880020.906714.933311.2000

91.266767.613343.960029.120023.093317.066712.8000

90.800068.240045.680031.360025.280019.200014.4000

90.333368.866747.400033.600027.466721.333316.0000

89.866769.493349.120035.840029.653323.466717.6000

89.400070.120050.840038.080031.840025.600019.2000

88.933370.746752.560040.320034.026727.733320.8000

88.466771.373354.280042.560036.213329.866722.4000

88.000072.000056.000044.800038.400032.000024.0000

86.600071.586756.573345.946739.706733.466725.8400

85.200071.173357.146747.093341.013334.933327.6800

83.800070.760057.720048.240042.320036.400029.5200

82.400070.346758.293349.386743.626737.866731.3600

81.000069.933358.866750.533344.933339.333333.2000

79.600069.520059.440051.680046.240040.800035.0400

78.200069.106760.013352.826747.546742.266736.8800

76.800068.693360.586753.973348.853343.733338.7200

75.400068.280061.160055.120050.160045.200040.5600

74.000067.866761.733356.266751.466746.666742.4000

72.600067.453362.306757.413352.773348.133344.2400

71.200067.040062.880058.560054.080049.600046.0800

69.800066.626763.453359.706755.386751.066747.9200

68.400066.213364.026760.853356.693352.533349.7600

67.000065.800064.600062.000058.000054.000051.6000

Columns 8 through 11

0000

1.06670.72000.56000.4000

2.13331.44001.12000.8000

3.20002.16001.68001.2000

4.26672.88002.24001.6000

5.33333.60002.80002.0000

6.40004.32003.36002.4000

7.46675.04003.92002.8000

8.53335.76004.48003.2000

9.60006.48005.04003.6000

10.66677.20005.60004.0000

11.73337.92006.16004.4000

12.80008.64006.72004.8000

13.86679.36007.28005.2000

14.933310.08007.84005.6000

16.000010.80008.40006.0000

18.213313.186710.76008.3333

20.426715.573313.120010.6667

22.640017.960015.480013.0000

24.853320.346717.840015.3333

27.066722.733320.200017.6667

29.280025.120022.560020.0000

31.493327.506724.920022.3333

33.706729.893327.280024.6667

35.920032.280029.640027.0000

38.133334.666732.000029.3333

40.346737.053334.360031.6667

42.560039.440036.720034.0000

44.773341.826739.080036.3333

46.986744.213341.440038.6667

49.200046.600043.800041.0000

双线性插值得到的钢轨温度立体图

b.插值基点为散乱节点

问题：已知n个节点，求点处的插值。

指令：cz=griddata(x,y,z,cx,cy,'method')

这里x,y,z都是n维向量，表明数据点的横坐标、纵坐标和竖坐标。向量cx,cy是给定的网格点的横坐标和纵坐标，指令cz=griddata(x,y,z,cx,cy,'method')返回在网格(cx,cy)处的函数值。cx与cy应是方向不同的向量，即一个是行向量，一个是列向量。

这里输入参数method有：

'linear'——双线性插值。默认值，即格式(3.1)

'cubic'——双三次插值

'nearest'——最邻近区域插值

例：在某海域测得一些点(x，y)处的水深z，在矩形区域(75，200)*(-50，150)内画出海底曲面图形(即设计航线问题——绘制等高线)

x 129 140 103.5 88 185.5 195 105 158 108 77 81 162 162 118

y 7.5 142 23 147 23 138 86 -6.5 -81 3 57 -66 84 -34

z 4 8 6 8 6 8 8 9 9 8 8 9 4 9

解：

x=[129 140 103.5 88 185.5 195 105 158 108 77 81 162 162 118];

y=[7.5 142 23 147 23 138 86 –6.5 –81 3 57 –66 84 –34];

z=[-4 –8 –6 –8 –6 –8 –8 –9 –9 –8 –8 –9 –4 –9];

cx=75:5:200;

cy=-70:5:150;

cz=griddata(x,y,z,cx,cy','cubic')

meshz(cx,cy,cz)

contour3(cx,cy,cz,16)

contour3(cx,cy,cz,16)

4、习题与上机实验

1)

t=[5 8 9 15 25 29 31 30 22 25 27 24];

x1=1:0.2:12;

y=interp1(h,t,x1,'spline');

plot(h,t,h,t,'+',x1,y,'b:')

legend('线性插值','原数据点','样条插值')

线性插值与样条插值的比较

2)对在[-5，5]上用年n(=11)个节点(等分)分别作Lagrage插值、分段线性插值和三次样条插值，用m(=21)个插值点(等分)作图，比较结果。

解：

n=11;m=21;

x=-5:10/(m-1):5;

y=1./(1+x.^2);

z=0*x;

x0=-5:10/(n-1):5;

y0=1./(1+x0.^2);

y1=lagrange(x0,y0,x);

y2=interp1(x0,y0,x);

y3=interp1(x0,y0,x,'spline');

[x',y',y1',y2',y3']

plot(x,z,'k',x,y,'b',x,y1,'c',x,y2,'k',x,y3,'b:')

legend('','y=1/(1+x^2)','Lagrange','Piece-linear','Spline')

-5.00000.03850.03850.03850.0385

-4.50000.04711.57870.04860.0484

-4.00000.05880.05880.05880.0588

-3.50000.0755-0.22620.07940.0745

-3.00000.10000.10000.10000.1000

-2.50000.13790.25380.15000.1401

-2.00000.20000.20000.20000.2000

-1.50000.30770.2353 0.35000.2973

-1.00000.50000.50000.50000.5000

-0.50000.80000.84340.75000.8205

01.00001.00001.00001.0000

0.50000.80000.84340.75000.8205

1.00000.50000.50000.50000.5000

1.50000.30770.23530.35000.2973

2.00000.20000.20000.20000.2000

2.50000.13790.25380.15000.1401

3.00000.10000.10000.10000.1000

3.50000.0755-0.22620.07940.0745

4.00000.05880.05880.05880.0588

4.50000.04711.57870.04860.0484

5.00000.03850.03850.03850.0385

3)机床加工问题

待加工零件的外形根据工艺要求由一组数据(x,y)给出(在平面情况下)，用程控铣床加工时每一刀只能沿x方向和y方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的(x,y)的坐标。下表给出机翼截面下轮廓线上的部分数据，假设需要得到x坐标每改变0.1时的y的坐标。试完成加工所需数据，画出曲线，并求出x=0处的曲线的斜率和13

机翼截面下轮廓线上的部分数据x 0 3 5 7 9 11 12 13 14 15

y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6

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);

y2=interp1(x0,y0,x);

y3=interp1(x0,y0,x,'spline');

[x',y1',y2',y3'];

subplot(3,1,1)

plot(x0,y0,'k+',x,y1,'r')

grid

title('Lagrange')

subplot(3,1,2)

plot(x0,y0,'k+',x,y2,'r')

grid

title('Piece-linear')

subplot(3,1,3)

plot(x0,y0,'k+',x,y3,'r')

grid

title('Spline')

可见Lagrange插值的结果根本不能用，分段线性插值光滑性较差(特别当x=14时)，建立采用三次样条插值。可以直接从数据结果中得到：x=0时的曲线斜率为0.0579/0.1=0.579；13

二、曲线(面)拟合法建模

1、数学原理

已知一组(二维)数据，即平面上的n个点互不相同。寻求一个函数(曲线)，使在某种准则下与所有数据点最为接近，即曲线拟合得最好。首先，确定所求曲线的形式(即经验公式)，线性最小二乘法是解决曲线拟合最常用的方法。其基本思路：令

其中是事先选定的一组函数，是待定系数(。拟合的准则是使n个点与的距离的平方和最小，称为最小二乘准则。

记

为求使J达到最小，只需利用极值的必要条件得到关于的超定线性方程组。

当时，即为多项式插值。

比较常用的经验公式：

多项式：(一般m=2,3，不宜过高)

特别地，一次多项式拟合又称为“线性回归”。

指数曲线：

双曲线(一支)：

注意：指数曲线拟合首先需作变量代换，化成多项式拟合进行。

2、在MATLAB中的实现

在MATLAB中，实现最小二乘拟合通常采用两种途径：

(1)利用polyfit函数进行多项式拟合。

P=polyfit(x,y,n)——用n阶多项式拟合x,y向量给定的数据

PA=polyval(p,xi)——求xi点上的拟合函数的近似值。

(2)利用常用的矩阵除法解决复杂型函数的拟合。

例1：以一次、二次和三次多项式拟合下数据x 0.5 1.0 1.5 2.0 2.5 3.0

y 1.75 2.45 3.81 4.80 7.00 8.60

解：

y=[1.75 2.45 3.81 4.80 7.00 8.60];

a1=polyfit(x,y,1)

a2=polyfit(x,y,2)

a3=polyfit(x,y,3)

x1=[0.5:0.05:3.0];

y1=a1(2)+a1(1)*x1;

y2=a2(3)+a2(2)*x1+a2(1)*x1.^2;

y3=a3(4)+a3(3)*x1+a3(2)*x1.^2+a3(1)*x1.^3;

plot(x,y,'*')

hold on

plot(x1,y1,'b:',x1,y2,'k',x1,y3,'g')

a2 =0.56140.82871.1560

a3 =-0.11561.1681-0.08711.5200

legend('原数据图','一次拟合','二次拟合','三次拟合')

1.24292.63974.03665.43346.83038.2271

1.71072.54613.66235.05916.73678.6950

1.75402.48553.62765.09386.79748.6517

v1=y-p1;

v2=y-p2;

v3=y-p3;

s1=norm(v1,'fro')

s2=norm(v2,'fro')

s3=norm(v3,'fro')

0.9558

s2 =

0.4220

s3 =

0.4057

可见，三次拟合方差最小。

例2：用6阶多项式对区间[0，2.5]上的误差函数进行最小二乘拟合。

解：

y=erf(x);%计算“误差函数”在[0,2.5]内的数据点

p=polyfit(x,y,6)

px=poly2str(p,'x')

0.0084-0.09830.4217-0.74350.14711.10640.0004

px =

0.0084194 x^6 - 0.0983 x^5 + 0.42174 x^4 - 0.74346 x^3 + 0.1471 x^2

+ 1.1064 x + 0.00044117

例3有效拟合的区间性图示(用[0,2.5]区间数据拟合曲线拟合[0,5]区间数据)

x1=0:0.1:2.5

y=erf(x);

y1=erf(x1);

p=polyfit(x1,y1,6)

f=polyval(p,x);

plot(x,y,'bo',x,f,'r-')

axis([0,5,0,2])

legend('拟合曲线','原数据线')

Columns 1 through 7

00.10000.20000.30000.40000.50000.6000

Columns 8 through 14

0.70000.80000.90001.00001.10001.20001.3000

Columns 15 through 21

1.40001.50001.60001.70001.80001.90002.0000

Columns 22 through 26

2.10002.20002.30002.40002.5000

p =

0.0084-0.09830.4217-0.74350.14711.10640.0004

说明：在[0,2.5]区间之外，两条曲线的表现完全不同。

例4：用最小二乘法求一个形如的经验公式，拟合下数据。x 19 25 31 38 44

y 19.0 32.3 49.0 73.3 97.8

下面介绍另一种方法解此拟合问题。用矩阵的方法求解，把a,b看成是为知量，已知，求解一超定方程：

在MATLAB中的实现：

y=[19.0 32.3 49.0 73.3 97.8];

x1=x.^2

x1=[ones(5,1),x1']

ab=x1\y'

x0=[19:0.2:44];

y0=ab(1)+ab(2)*x0.^2;

clf

plot(x,y,'bo')

hold on%图形合在同一坐标下

plot(x0,y0,'r-')

36162596114441936

x1 =

1361

1625

1961

11444

11936

ab =

0.9726

0.0500

3、习题与上机实验

1)已知数据表如下，试利用二次多项式拟合下数据，求出xi=1,1.5,2,2.5,…,5各点的函数的近似值，并与用interp1的线性插值算法求xi=1,1.5,2,2.5,…,5各点的函数的近似值作比较。X 1.0 2.0 3.0 4.0 5.0

Y 3.5 4.6 5.5 3.2 2.0

解：%断开与上例题“hold on”

xx=1:1:5;

yx=[3.5,4.6,5.5,3.2,2];

xxi=1:0.5:5;

ff0=interp1(xx,yx,xxi)

yy=polyfit(xx,yx,2);

ff1=polyval(yy,xxi)

plot(xx,yx,'*',xxi,ff0,'b:',xxi,ff1,'k')

Columns 1 through 7

3.50004.05004.60005.05005.50004.35003.2000

Columns 8 through 9

2.60002.0000

ff1 =

Columns 1 through 7

3.52574.28074.75714.95504.87434.51503.8771

Columns 8 through 9

2.96071.7657

legend('线性插值','二次拟合')

2)分别用线性回归、二次多项式和十次多项式拟合下数据：x 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

y -0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2

y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];

p1=polyfit(x,y,1)

p2=polyfit(x,y,2)

p3=polyfit(x,y,10)

xi=linspace(0,1,100);

z1=polyval(p1,xi);

z2=polyval(p2,xi);

z3=polyval(p3,xi);

plot(x,y,'*',xi,z1,'b:',xi,z2,'g',xi,z3,'k')

10.31851.4400

p2 =

-9.810820.1293-0.0317

p3 =

1.0e+006 *

Columns 1 through 7

-0.46442.2965-4.87735.8233-4.29482.0211-0.6032

Columns 8 through 11

0.1090-0.01060.0004-0.0000

legend('原数据图','线性拟合','二次拟合','高阶拟合')

可见，高次拟合曲线在数据点附近更加接近数据点的测量值，但从整体上来说，曲线波动比较大，并不一定适合实际应用的需要。所以，在进行高次曲线拟合时，“越高越好”的观点是不一定对的。

此文完成于2001年暑假。

王 正 盛