卷帘工程验收单怎么写:插值与拟合matlab实现

插值与拟合的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

hold off  

x=[-5:1:5];

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

f0 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

    2.6000    2.0000

f1 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

    2.6000    2.0000

f2 =

  Columns 1 through 7

    3.5000    4.0750    4.6000    5.2625    5.5000    4.4813    3.2000

  Columns 8 through 9

    2.4625    2.0000

f3 =

  Columns 1 through 7

    3.5000    3.7734    4.6000    5.3766    5.5000    4.5953    3.2000

  Columns 8 through 9

    2.0797    2.0000

f4 =

  Columns 1 through 7

    3.5000    4.6000    4.6000    5.5000    5.5000    3.2000    3.2000

  Columns 8 through 9

    2.0000    2.0000  

plot(xxi,f1,xxi,f2,xxi,f3,xxi,f4) 

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列矩阵，第一列存储室内温度，第二列存储室外温度。

h=6:2:18;

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

yi =

   18.5020   15.6553

   20.4986   20.3355

   22.5193   24.9089

   26.3775   29.6383

   30.2051   34.2568

   26.8178   30.9594  

plot(h,t,xi,yi)  

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

MATLAB\BIN\lag2.m

以例1数据为例：

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

f5=lag2(xx,yx,xxi)

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

f5 =

  Columns 1 through 7

    3.5000    4.0750    4.6000    5.4500    5.5000    4.2125    3.2000

  Columns 8 through 9

    2.4625    2.0000

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

（2）           三次样条插值

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

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

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

例：

 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,'spline')

  yi=spline(xx,yx,xxi)  

f0 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

2.6000                2.0000

f1 =

  Columns 1 through 7

    3.5000    3.7734    4.6000    5.3766    5.5000    4.5953    3.2000

  Columns 8 through 9

    2.0797    2.0000  

yi =

  Columns 1 through 7

    3.5000    3.7734    4.6000    5.3766    5.5000    4.5953    3.2000

  Columns 8 through 9

    2.0797    2.0000  

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

      (2.2) pp=spline(x,y)

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

例：

 xx=1:1:5;

 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)  

f0 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

    2.6000    2.0000

yyi =

  Columns 1 through 7

    3.5000    3.7734    4.6000    5.3766    5.5000    4.5953    3.2000

  Columns 8 through 9

    2.0797    2.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。

解：

x=0:2.5:10

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)  

x =

         0    2.5000    5.0000    7.5000   10.0000

h =

     0

    30

    60

T =

    95    14     0     0     0

    88    48    32    12     6

    67    64    54    48    41

Ti =

  Columns 1 through 7

   95.0000   62.6000   30.2000   11.2000    5.6000         0         0

   94.5333   63.2267   31.9200   13.4400    7.7867    2.1333    1.6000

   94.0667   63.8533   33.6400   15.6800    9.9733    4.2667    3.2000

   93.6000   64.4800   35.3600   17.9200   12.1600    6.4000    4.8000

   93.1333   65.1067   37.0800   20.1600   14.3467    8.5333    6.4000

   92.6667   65.7333   38.8000   22.4000   16.5333   10.6667    8.0000

   92.2000   66.3600   40.5200   24.6400   18.7200   12.8000    9.6000

   91.7333   66.9867   42.2400   26.8800   20.9067   14.9333   11.2000

   91.2667   67.6133   43.9600   29.1200   23.0933   17.0667   12.8000

   90.8000   68.2400   45.6800   31.3600   25.2800   19.2000   14.4000

   90.3333   68.8667   47.4000   33.6000   27.4667   21.3333   16.0000

   89.8667   69.4933   49.1200   35.8400   29.6533   23.4667   17.6000

   89.4000   70.1200   50.8400   38.0800   31.8400   25.6000   19.2000

   88.9333   70.7467   52.5600   40.3200   34.0267   27.7333   20.8000

   88.4667   71.3733   54.2800   42.5600   36.2133   29.8667   22.4000

   88.0000   72.0000   56.0000   44.8000   38.4000   32.0000   24.0000

   86.6000   71.5867   56.5733   45.9467   39.7067   33.4667   25.8400

   85.2000   71.1733   57.1467   47.0933   41.0133   34.9333   27.6800

   83.8000   70.7600   57.7200   48.2400   42.3200   36.4000   29.5200

   82.4000   70.3467   58.2933   49.3867   43.6267   37.8667   31.3600

   81.0000   69.9333   58.8667   50.5333   44.9333   39.3333   33.2000

   79.6000   69.5200   59.4400   51.6800   46.2400   40.8000   35.0400

   78.2000   69.1067   60.0133   52.8267   47.5467   42.2667   36.8800

   76.8000   68.6933   60.5867   53.9733   48.8533   43.7333   38.7200

   75.4000   68.2800   61.1600   55.1200   50.1600   45.2000   40.5600

   74.0000   67.8667   61.7333   56.2667   51.4667   46.6667   42.4000

   72.6000   67.4533   62.3067   57.4133   52.7733   48.1333   44.2400

   71.2000   67.0400   62.8800   58.5600   54.0800   49.6000   46.0800

   69.8000   66.6267   63.4533   59.7067   55.3867   51.0667   47.9200

   68.4000   66.2133   64.0267   60.8533   56.6933   52.5333   49.7600

   67.0000   65.8000   64.6000   62.0000   58.0000   54.0000   51.6000

  Columns 8 through 11

         0         0         0         0

    1.0667    0.7200    0.5600    0.4000

    2.1333    1.4400    1.1200    0.8000

    3.2000    2.1600    1.6800    1.2000

    4.2667    2.8800    2.2400    1.6000

    5.3333    3.6000    2.8000    2.0000

    6.4000    4.3200    3.3600    2.4000

    7.4667    5.0400    3.9200    2.8000

    8.5333    5.7600    4.4800    3.2000

    9.6000    6.4800    5.0400    3.6000

   10.6667    7.2000    5.6000    4.0000

   11.7333    7.9200    6.1600    4.4000

   12.8000    8.6400    6.7200    4.8000

   13.8667    9.3600    7.2800    5.2000

   14.9333   10.0800    7.8400    5.6000

   16.0000   10.8000    8.4000    6.0000

   18.2133   13.1867   10.7600    8.3333

   20.4267   15.5733   13.1200   10.6667

   22.6400   17.9600   15.4800   13.0000

   24.8533   20.3467   17.8400   15.3333

   27.0667   22.7333   20.2000   17.6667

   29.2800   25.1200   22.5600   20.0000

   31.4933   27.5067   24.9200   22.3333

   33.7067   29.8933   27.2800   24.6667

   35.9200   32.2800   29.6400   27.0000

   38.1333   34.6667   32.0000   29.3333

   40.3467   37.0533   34.3600   31.6667

   42.5600   39.4400   36.7200   34.0000

   44.7733   41.8267   39.0800   36.3333

   46.9867   44.2133   41.4400   38.6667

   49.2000   46.6000   43.8000   41.0000

mesh(xi,hi,Ti)

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

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

解：

clear

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）

h=1:12;

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)个插值点（等分）作图，比较结果。

解：

hold off

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

ans =

   -5.0000    0.0385    0.0385    0.0385    0.0385

   -4.5000    0.0471    1.5787    0.0486    0.0484

   -4.0000    0.0588    0.0588    0.0588    0.0588

   -3.5000    0.0755   -0.2262    0.0794    0.0745

   -3.0000    0.1000    0.1000    0.1000    0.1000

   -2.5000    0.1379    0.2538    0.1500    0.1401

   -2.0000    0.2000    0.2000    0.2000    0.2000

   -1.5000    0.3077    0.2353    0.3500    0.2973

   -1.0000    0.5000    0.5000    0.5000    0.5000

   -0.5000    0.8000    0.8434    0.7500    0.8205

         0    1.0000    1.0000    1.0000    1.0000

    0.5000    0.8000    0.8434    0.7500    0.8205

    1.0000    0.5000    0.5000    0.5000    0.5000

    1.5000    0.3077    0.2353    0.3500    0.2973

    2.0000    0.2000    0.2000    0.2000    0.2000

    2.5000    0.1379    0.2538    0.1500    0.1401

    3.0000    0.1000    0.1000    0.1000    0.1000

    3.5000    0.0755   -0.2262    0.0794    0.0745

    4.0000    0.0588    0.0588    0.0588    0.0588

    4.5000    0.0471    1.5787    0.0486    0.0484

    5.0000    0.0385    0.0385    0.0385    0.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

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;

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

解：

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

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

a1 =    2.7937   -0.1540

a2 =    0.5614    0.8287    1.1560

a3 =   -0.1156    1.1681   -0.0871    1.5200

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

p1=polyval(a1,x)  

p1 =

    1.2429    2.6397    4.0366    5.4334    6.8303    8.2271  

p2=polyval(a2,x)

p2 =

    1.7107    2.5461    3.6623    5.0591    6.7367    8.6950  

  p3=polyval(a3,x) 

p3 =

1.7540    2.4855    3.6276    5.0938    6.7974    8.6517

   v1=y-p1;

 v2=y-p2;

 v3=y-p3;

  s1=norm(v1,'fro')

  s2=norm(v2,'fro')

  s3=norm(v3,'fro')  

s1 =

    0.9558

s2 =

    0.4220

s3 =

    0.4057  

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

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

解：

x=0:0.1:2.5;

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

p=polyfit(x,y,6)

px=poly2str(p,'x')

p =

    0.0084   -0.0983    0.4217   -0.7435    0.1471    1.1064    0.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]区间数据）

x=0:0.1: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('拟合曲线','原数据线')  

x1 =

  Columns 1 through 7

         0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000

  Columns 8 through 14

    0.7000    0.8000    0.9000    1.0000    1.1000    1.2000    1.3000

  Columns 15 through 21

    1.4000    1.5000    1.6000    1.7000    1.8000    1.9000    2.0000

  Columns 22 through 26

    2.1000    2.2000    2.3000    2.4000    2.5000

p =

0.0084   -0.0983    0.4217   -0.7435    0.1471    1.1064    0.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中的实现：

x=[19 25 31 38 44];

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

x1 =

         361         625         961        1444        1936

x1 =

           1         361

           1         625

           1         961

           1        1444

           1        1936

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”

 hold off 

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

ff0 =

  Columns 1 through 7

    3.5000    4.0500    4.6000    5.0500    5.5000    4.3500    3.2000

  Columns 8 through 9

    2.6000    2.0000

ff1 =

  Columns 1 through 7

    3.5257    4.2807    4.7571    4.9550    4.8743    4.5150    3.8771

  Columns 8 through 9

    2.9607    1.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

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

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

p1 =

   10.3185    1.4400

p2 =

   -9.8108   20.1293   -0.0317

p3 =

  1.0e+006 *

  Columns 1 through 7

   -0.4644    2.2965   -4.8773    5.8233   -4.2948    2.0211   -0.6032

  Columns 8 through 11

    0.1090   -0.0106    0.0004   -0.0000

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

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

                                      此文完成于2001年暑假。

                                           王 正 盛