数值分析上机作业详细分析

函数插值与曲线拟合

实验2.1 Runge现象及如何避免

考虑区间[-1,1]的一个等距划分，分点为

2i

x i =-1+, i =0, 1, 2, , n

n

则拉格朗日插值多项式为

L n (x ) =∑

1

l (x ) 2i

1+25x i =0j

n

其中的l i (x ), i =0, 1, 2, , n 是n 次拉格朗日插值基函数。

（1） 选择不断增大的分点数目n=2,3…. ，画出原函数f(x)及插值多项式函数

L n (x ) 在[-1，1]上的图像，比较并分析实验结果。

（2）选择其他的函数，例如定义在区间[-5，5]上的函数 x

h (x ) =, g (x ) =arctan x 4

1+x

重复上述的实验，看其结果如何。

（3）区间[a,b]上切比雪夫点的定义为

x k =

⎛(2k -1) πb +a b -a

+c o s 22⎝2(n +1) ⎫

⎪⎪, k =1, 2, , n +1 ⎭

以x 1, x 2, x n +1为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果, 并分析原因。 数值计算结果

（1）首先，对函数在[-1,1]上进行拉格朗日插值，算出各阶的插值函数，然后在区间上等距地取100个点计算插值结果和理论结果，然后做图和求二者之差绝对值的无穷范数。下图1图2给出了对龙格函数进行不同阶数等距插值的插值结果。下表1为龙格函数不同阶数插值时插值结果与理论值偏差的无穷范数。

表1 插值结果与理论值偏差的无穷范数

图1 龙格函数不同阶等距插值结果

从图1我们已经可以看到龙格现象，例如20阶插值比10阶插值得到的插值函数在插值区间出现了更严重的震荡，它们偏离理论值得无穷范数分别为0.2979和8.1906。即随着插值点的增多Ln(x)没有一致趋近f(x)，而是在部分区间出现了严重的偏离被插值函数的震荡现象，这就是龙格现象。

为了进一步观察龙格现象，下面图2、图3和图4分别显示了11~18阶，33~40阶和40~47阶的等距插值结果。试图看看会出现什么结果。

图2龙格函数不同阶等距插值结果

图3 龙格函数不同阶等距插值结果

图4 龙格不同阶等距插值结果

从图2可以看明显地看到龙格现象，从11阶、13阶一直到17阶奇数阶随着阶数的增大震荡的峰值越来越大；同时，12阶、14阶一直到18阶偶数阶随着阶数的增大震荡的峰值同样也是越来越到。我们清楚地看到了龙格现象。

然而，随着阶数的大幅升高，图3和图4中插值函数的震荡现象变得不那么理想，变得没什么规律，甚至看不出震荡现象，只是在左端发生了偏移。例如40阶插值时只有左端发生了向下的极大偏移，而且图像左右极不对称。然而，从理论上来讲，插值函数左右应该是对称的。为什么会出现这种现象？此现象由什么引起？它对我们判断龙格现象又会产生什么样的干扰？

实际上，从40阶到47阶插值图像的变化情况我们已经不能判断是否有龙格现象了。为了说明这个结论，我们继续进行插值实验。下面图5到图9分别是10阶、18阶、20阶和30阶

从图3和图4好像也观察不出什么规律。下面继续作进一步计算。下面图5、图6、图、7图8和图9分别作了10阶、18阶、20阶、30阶和40阶插值。

图5 龙格函数10阶等距插值

图6 龙格函数18阶等距插值

图7 龙格函数20阶等距插值

图8 龙格函数30阶等距插值

从上图5可以看出左右端点处（±1）插值结果和理论结果相同，同为0.03846。但是，通过Matlab 中fprintf 函数控制位数输出，精确到小数点后15发现二者在更高精度下是不等的，理论值为0.[**************]，而插值函数计算结果为0.[**************]。由此可以看出计算得出的插值函数系数项与理论上插值函数的系数值存在舍入误差。这种误差在插值阶数不高时表现不明显，可是当阶数高到30阶时这种舍入误差带来的误差积累结果已经开始显著，如图8所示，左右端点的插值函数计算结果分别为0.6047和0.03849，本应相等的两个值，现在已经出现明显的偏差。同样如图1所示，当阶数到40时，数值计算过程中的舍入误差的积累已经相当剧烈，左右表现出了极大的不对称性。

在自变量更大的范围内，插值函数是不是仍然对称呢？下面图9在自变量[-1.5,1.5]范围内画出了[-1,1]范围内的插值函数的图像。此时可以观察到插值函数的对称性。通过图上给的几个函数值可以看出，此时的插值函数，因为阶数很高，对自变量极其敏感。尤其在插值区间外插值函数变化剧烈，严重偏离龙格函数。

图9 龙格函数40阶等距插值

通过上面的数值计算和分析，我们得到如下几个结论:（a ）在[-1,1]范围内龙格函数的拉格朗日插值函数存在龙格现象；实际上龙格已经证明在大约[-0.72,0.72]或以0为中心更小的区间上对龙格函数进行拉格朗日插值才是收敛的，即没有龙格现象。（b ）插值区间是影响龙格现象的重要因素之一。后面还会看到插值点的选取也是影响龙格现象的重要因素之一。（c ）决定数值计算结果的因素除了算理以外还与计算机的舍入误差有关，当阶数很高时，积累起来的舍入误差会大到淹没算理的程度。此时再利用数值计算的结果分析算理是不合适的。具体到本实验，当阶数很高时，插值函数与理论函数产生的数值偏差由龙格现象和数值计算过程中的舍入误差积累共同产生，因此我们不能根据高阶时的数值计算结果来判断到底是不是龙格现象。数值计算中我们应该在保证舍入误差的积累影响不明显时，进行不同阶数的插值进行判断是否存在龙格现象。（d ）计算机在进行高阶或多次计算时，舍入误差可能多实验结果产生重大影响。

（2）下面给出了h(x)和g(x)的等距插值结果。

图10为h(x)函数在[-1,1]上不同阶等距插值结果，图11为h(x)在[-1,1]上30阶等距插值的结果。表2、表3、表4为h(x) 在[-1,1]上各阶插值结果与理论值偏差的无穷范数。

表2 h(x) 各阶插值结果与理论值偏差的无穷范数

表3 h(x) 各阶插值结果与理论值偏差的无穷范数

表4 h(x) 各阶插值结果与理论值偏差的无穷范数

图10 h(x)函数不同阶等距插值

图11 h(x)函数30阶等距插值

从表2到表4和图10 和图11结合上面的结论可以判断，此情况下并没有出现龙格现象。图11中没有出现震荡现象，左端的偏差是数值计算中舍入误差的累积产生的。可以不严格地设想，如果没有舍入误差，左端插值结果就会上移到理论值，同时左端右侧的偏离也会被带动到理论值，即没有龙格现象。

为了研究h(x)等距插值的收敛区间，即不发生龙格现象的区间大小。下面继续进行实验，改变插值区间的大小和插值阶数。

图11 h(x)在区间[-5,5]上的30阶等距插值

图12 h(x)在区间[-2,2]上的30阶等距插值

图13 h(x)在区间[-2,2]上的等距插值

图14 h(x)在区间[-1.5,1.5]上的30阶等距插值

图15 h(x)在区间[-1.5,1.5]上的33阶等距插值

图11到图15分别计算了区间[-5,5]，[-2,2]和[-1.5,1.5]区间上30阶和33阶的插值结果。从图11可以看出h(x)在区间[-5,5]上30阶插值具有明显的龙格现象。然和缩小区间，在[-2,2]上进行30阶等距插值，发现龙格现象比在区间[-5,5]上减小了10000被左右，但当插值阶数增加到33时，震荡加剧，说明此时仍有明显的龙格现象。继续缩小区间，在区间[-1.5,1.5]上进行30阶等距插值，发现龙格现象几乎为零；增加插值阶数到33发现震荡现象没有增加，说明此时基本没有了龙格现象。即函数h(x)的拉格朗日等距插值收敛区间大约为[-1.5,1.5]。

下图16为g(x)在区间[-1,1]上各阶等距插值图像，图17为g(x)在区间[-1,1]上30阶等距插值图像。下表5、表6和表7为g(x)在区间[-1,1]上的各阶等距插值与理论值之差的绝对值的无穷范数。

表5 函数g(x)各阶等距插值与理论值偏差的无穷范数

表6 函数g(x)各阶等距插值与理论值偏差的无穷范数

表7 函数g(x)各阶等距插值与理论值偏差的无穷范数

图16 g(x)函数不同阶等距插值

图17 g(x)函数30阶等距插值

从上表5、6、7和图16、16可以看出函数g(x)在区间[-1,1]上的等距拉格朗日插值没有龙格现象。为了寻找函数g(x)的等距插值收敛区间，继续进行数值实验。

图18 函数g(x)在区间[-5,5]上的30阶等距插值

图19 函数g(x)在区间[-3,3]上的30阶等距插值

图20 函数g(x)在区间[-2.5,2.5]上的30阶等距插值

图21 函数g(x)在区间[-2.5,2.5]上的33阶等距插值

从图18、19可以看出函数g(x)在区间[-5,5]和[-3,3]上的等距拉格朗日插值都是存在龙格现象的，即插值函数不收敛到被插值函数。继续缩小插值区间，从图20和图21可以看出在区间[-2.5,2.5]上的30阶插值稍微有一点震荡，而插值阶数增加到33阶时震荡基本被抑制。即函数g(x)的等距拉格朗日插值收敛区间大约为[-2.5,2.5]。

（3）利用切比雪夫点对上面三个函数进行拉格朗日插值

从上面图8、图11和图18中可以看到，函数f(x)在区间[-1,1]的30阶等距拉格朗日插值，函数h(x)在区间[-5,5]的30阶等距拉格朗日插值和函数g(x)在区间[-5,5]的30阶等距拉格朗日插值都存在严重的龙格现象。为了研究切比雪夫点拉格朗日插值的性质，下面进行进行等距插值和切比雪夫点插值的比较。

首先研究函数f(x)。下图22为函数f(x)在区间[-1,1]上的30阶等距和切比雪夫点插值 ，图23为函数f(x)在区间[-1,1]上的33阶等距和切比雪夫点插值，图24为函数f(x)在区间[-20,20]上的30阶等距和切比雪夫点插值。

图22 函数f(x)在区间[-1,1]上的30阶等距和切比雪夫点插值

图23 函数f(x)在区间[-1,1]上的33阶等距和切比雪夫点插值

图24 函数f(x)在区间[-20,20]上的30阶等距和切比雪夫点插值

对比图22和图23，首先分别对比30阶插值时，等距插值出现严重龙格现象而切比雪夫点插值没有出现龙格现象，33阶插值时出现相同现象。然后比较30阶和33阶的等距插值随插值节点增多，龙格现象越来越严重；然而切比雪夫点插值一直没有出现龙格现象。可见，对于函数f(x)，切比雪夫点插值可以有效避免龙格现象。

为了考察更大插值区间内，切比雪夫点插值是不是可以有效避免龙格现象，上面图24对区间[-20,20]内函数f(x)进行了30阶的等距了切比雪夫点插值。发现，在更大插值区域内，切比雪夫点插值仍然可以有效避免龙格现象，然而等距插值的龙格现象加剧。

然后研究函数h(x)。下图25为函数h(x)在区间[-5,5]上的30阶等距和切比雪夫点插值 ，图26为函数h(x)在区间[-5,5]上的33阶等距和切比雪夫点插值，图27(1)(2)分别为函数h(x)在区间[-20,20]上的30阶、33阶等距和切比雪夫点插值。

图25 函数h(x)在区间[-5,5]上的30阶等距和切比雪夫点插值

图26 函数h(x)在区间[-5,5]上的33阶等距和切比雪夫点插值

图27(1) 函数h(x)在区间[-20,20]上的30阶等距和切比雪夫点插值

图27(2) 函数h(x)在区间[-20,20]上的30阶等距和切比雪夫点插值

对比图25中等距插值和切比雪夫点插值，等距插值存在龙格现象，而切比雪夫点插值不存在，图26相同。对比图25和图26，当增大插值阶数时，等距插值龙格现象加剧，而切比雪夫点插值仍然没有龙格现象。说明对于函数h(x)可以有效避免龙格现象。

同样为了观察更大插值区间，切比雪夫点插值对龙格现象的影响，进行了区间[-20,20]上的30阶、33阶等距和切比雪夫点插值。发现两种情况下，等距插值龙格现象严重，而33阶的切比雪夫点插值相比30阶的切比雪夫点插值函数是收敛的，所以在此区间函数h(x)的切比雪夫点插值没有龙格现象。

最后研究函数g(x)。下图28为函数g(x)在区间[-5,5]上的30阶等距和切比雪夫点插值，图29为函数g(x)在区间[-5,5]上的33阶等距和切比雪夫点插值，图30为函数g(x)在区间[-20,20]上的30阶等距和切比雪夫点插值。

图28 函数g(x)在区间[-5,5]上的30阶等距和切比雪夫点插值

图29 函数g(x)在区间[-5,5]上的33阶等距和切比雪夫点插值

图30 函数g(x)在区间[-20,20]上的30阶等距和切比雪夫点插值

对比图28和图29，首先分别对比30阶插值时，等距插值出现严重龙格现象而切比雪夫点插值没有出现龙格现象，33阶插值时出现相同现象。然后比较30阶和33阶的等距插值随插值节点增多，龙格现象依然严重；然而切比雪夫点插值一直没有出现龙格现象。可见，对于函数g(x)，切比雪夫点插值可以有效避免龙格现象。

为了考察更大插值区间内，切比雪夫点插值是不是可以有效避免龙格现象，上面图30对区间[-20,20]内函数g(x)进行了30阶的等距了切比雪夫点插值。发现，在更大插值区域内，切比雪夫点插值仍然可以有效避免龙格现象，然而等距插值的龙格现象加剧。

综合上面对函数f,h,g 的切比雪夫点插值的计算结果来看，切比雪夫点插值可以有效地避免龙格现象。理论上的原因如下。

有一个定理如下，如果插值节点为切比雪夫点，被插值函数f ∈C n+1[−1,1]，Ln(x)为相应的插值多项式，则

1

max |f−Ln|≤||fn+1||∞ −1≤x≤1对于一般区间可以通过伸缩和平移变换到区间[-1,1]上来。

可以看到，插值函数与理论值的偏差被右端项控制。对于函数f,g,h 的高阶导数值得最大值变化缓慢，远比2n (n +1) 慢，随插值点增多，右端项分母快速增大，右端项很快趋近于0，即可以保证当插值点趋向无穷时，插值函数收敛到被插值函数。这也是为什么切比雪夫点插值可以避免龙格现象。实际上，它可以保证在全空间上插值函数收敛到被插值函数。

实验2.2样条函数插值

（1）考虑函数

f (x ) =

x

−5≤x ≤5 1+x 分别选择等距节点和切比雪夫点进行样条函数插值，插值节点个数分别取10，20。随节点

个数增加，比较被逼近函数和样条插值函数误差的变化情况，并与拉格朗日多项式插值比较(可以用MATLAB 的函数“spline ”作此函数的三次样条插值，取n=10、20，分别画出插值函数及原函数的图形) 。 （2）样条插值的思想是早产生于工业部门。作为工业应用的例子考虑如下问题：

ii. 主函数myspline(x,y,边界类型，边界值，x i ) 其中：x 节点 y 节点上的函数值 x i 未知节点 返回：S(xi )

iii. 三对角方程组用追赶法求解(书P160) 。

数值计算过程及结果分析

（1）图31为函数f(x)在[-5,5]上的10阶样条插值和拉格朗日插值结果，图32为函数f(x)在[-5,5]上的20阶样条插值和拉格朗日插值结果。下表8给出了不同阶插值不同插值方式下插值结果与理论值误差向量的无穷范数。

从图31、32和表8的数据可以看出，样条函数插值时，等距节点比切比雪夫节点插值结果

更好；拉格朗日插值时切比雪夫点插值比等距插值效果要好。从表8可以看出随阶数增大，等距样条插值与理论值的偏差收敛到0。

图31 函数f(x)在[-5,5]上的10阶样条插值和拉格朗日插值

图32 函数f(x)在[-5,5]上的20阶样条插值和拉格朗日插值

（2）三弯矩方程组求解样条插值问题。利用自编样条插值函数Myspline 函数进行样条插值。同时利用Matlab 库函数spline 进行插值作为对照，已验证Myspline 函数的可靠性。下图33为利用spline 函数和Myspline 函数对题目中第一类边界问题进行的三次样条插值结果，二者偏差的无穷范数为4.4409×10-16。图34为利用spline 函数和Myspline 函数对题目中第二类边界问题（自由边界条件）进行的三次样条插值结果，二者偏差的无穷范数为0.0149。通过上面的计算可以证明Myspline 函数的可靠性。

图33 用spline 函数和Myspline 函数对题目中第一类边界问题进行的三次样条插值结果

图34 利用spline 函数和Myspline 函数对题目中第二类边界问题（自由边界条件）进行的三

次样条插值结果

实验2.3一维插值的应用-作图

利用Matlab 对自己手型进行屏幕采样，然后利用样条插值画出自己的手型。 实验过程如下：

首先利用下面一段Matlab 的代码进行手型屏幕采样。采集了40个样点。

figure('position',get(0,'screensize')) axes('position',[0 0 1 1]) [x,y]=ginput;

然后利用函数spline 进行三次样条插值。由于图形比较复杂，图形走势不满足函数一一对应要求。需要定义第三个变量，建立某种一一对应的关系。令采样点的顺序为自变量，次数采样点可以是小数。这样构造满足同一个点的横坐标和纵坐标对应同一个自变量，而不必在意自变量有何性质。下面是三次样条插值的代码：

tri=1:1:40;%第三个自变量

pp1=spline(tri,x);%横坐标与自变量的样条插值函数 pp2=spline(tri,y);%纵坐标与自变量的样条插值函数 N=400;%重新绘制手型时的样点数

trig=linspace(1,40,N);%把自变量在1到40之间平分为400份 xg=ppval(pp1,trig);%计算插值横坐标 yg=ppval(pp2,trig);%计算插值纵坐标 plot(xg,yg);%绘图

最终绘制结果如图35所示。

图35 手型效果图

思考题：二维插值

（1）在一丘陵地带测量高程，x 和y 方向每隔100米测一个点，得高程数据如下。试用MATLAB 的二维插值函数“interp2”进行插值，并由此找出最高点和该

利用MATLAB 的peaks 函数生成某山区的一些地点及其高度三维数据（单位：m ）。命令格式：[x,y,z]=peaks(n)，生成的n 阶矩阵x,y,z 为测量的山区地点三维数据（n>=30）。根据peaks 函数生成的数据，利用Matlab 二维插值画出该山区的地貌图和等值线图(提示函数：interp2、meshgrid 、plot3等) 。

数值计算过程及结果 （1）通过interp2插值得到图37所示的地形图。最高点坐标为(167.3,179.6)，高程为721.1。图36为未插值时的地形图。图37为插值后的地形图。

图36 未插值时的地形图

图37 插值后的地形图

由于代码比较简短，故将Matlab 代码附于下面：

x=[100 200 300 400];%测量数据 y=[100；200；300；400]; z=[636 697 624 478 698 712 630 478 680 674 598 412 662 626 552 334];

temp=linspace(100,400,50);%插值点个数 [xx,yy]=meshgrid(temp,temp);%生成插值网格 zz=interp2(x,y,z,xx,yy,'spline' );%进行插值 surf(xx,yy,zz);%绘图 m=max(max(zz));%寻找最大值

[xxx,yyy]=find(zz==m);%确定最大值位置, 返回行数和列数 Xmax=xx(1,yyy);%确定最大值的横坐标 Ymax=yy(xxx,1);%确定最大值得纵坐标

（2）图38是利用Matlab 做出的某一地区的地貌图。图39为图38地形对应的等高线。

图38 某地区地貌图

Matlab 代码如下：

[x,y,z]=peaks(40); temp=linspace(100,400,50); [xx,yy]=meshgrid(temp,temp); zz=interp2(x,y,z,xx,yy,'spline' ); surf(xx,yy,zz);

图39 图38地形对应的等高线

Matlab 代码如下：

[x,y,z]=peaks(40); temp=linspace(100,400,50); [xx,yy]=meshgrid(temp,temp); zz=interp2(x,y,z,xx,yy,'spline' ); contour(x,y,z,5);%参数5为指定的条数

%contour(x,y,z,[1,4]); %参数1,4为指定的等高线的高度

实验3.1曲线逼近方法的比较

考虑实验2.1中的著名问题。下面的MATLAB 程序给出了该函数的二次和三次拟合多项式。

x=-1:0.2:1; ylabel(‘y’); y=1./(1+25*x.*x); hold on; xx=-1:0.02:1; p3=polyfit(x,y,3); p2=polyfit(x,y,2); yy=polyval(p3,xx); yy=polyval(p2,xx); plot(x,y,’o’,xx,yy); plot(x,y,’o’,xx,yy); hold off; xlabel(‘x’);

适当修改上述MATLAB 程序，也可以拟合其他你有兴趣的函数。实验要求： （1）将拟合的结果与拉格朗日插值及样条插值的结果比较。

（2）归纳总结数值实验结果，试定性地说明函数逼近各种方法的适用范围，及实际应用中选择方法应注意的问题。 实验过程及结果

（1）下图40为对龙格函数进行的不同阶多项式拟合。

图40 对龙格函数进行的不同阶多项式拟合

对于龙格函数，由于为偶函数，所以进行多项式拟合时，基函数只有偶次项，因此n （为偶数）阶拟合和n+1阶的拟合结果相同。

对比图40与拉格朗日插值结果和样条插值结果。首先比较多项式拟合的结果与拉格朗日插值结果，当离散点数目n 一定时，随着拟合阶数的增加，拟合结果逐渐向拉格朗日插值结果靠拢，当拟合阶数为n-1时，最佳拟合结果与拉格朗日插值的结果相同。并且，如果n 比较大，高次插值和拟合都会出现龙格现象。还有，阶数比较低时，同样阶数的拟合结果要比插值结果更趋向理论曲线。

比较拟合结果与样条插值结果，低阶的拟合结果要比样条插值结果好，但高阶时样条插值结果要比拟合结果好。

（2）从离散点的数据来源来说，当测量数据精确时，插值更合理一些；当测量数据误差比较显著时，拟合比较合理一些。因为插值函数一定过离散点，在离散点没有误差，而拟合不一定过离散点，存在误差。 另一方面，当离散数据点比较少且精确时用低次拉格朗日插值比较合理；当离散点数据少且存在一定误差时用低次拟合比较好；当数据点比较多且精确度高时，用样条插值比较好；当数据点比较多且存在一定误差时用拟合方法更合理。

一般情况，对实验数据进行拟合时, 函数模型通常已知, 仅需要拟合得到参数值。如果函数模型未知可以通过多次尝试寻找最佳函数模型。

实验3.2最小二乘拟合的经验公式和模型

1. （已知经验公式）：某类疾病发病率为y ‰和年龄段x (每五年为一段，例如0~5

bx

y =ae 岁为第一段，6~10岁为第二段……)之间有形如的经验关系，观测得到的

bx

y =ae （1）用最小二乘法确定模型中的参数a 和b (提示函数：lsqcurvefit ，

lsqnonlin ) 。

bx

y =ae （2）利用MATLAB 画出离散数据及拟合函数图形。

（3）利用MATLAB 画出离散点处的误差图，并计算相应的均方误差。

2．（最小二乘拟合模型未知） 某年美国轿车价格的调查资料如表，其中x i 表示轿车的使用年数，y i 表示相应的平均价格。实验要求：试分析用什么形式的曲线来拟合表中的数据，并预测使用1下面进行已知经验公式的拟合。

bx

y ae （1）确定模型中的参数a 和b ，Matlab 代码如下：

xd=[];

yd=[];

f=@(a,xd)a(1)*exp(a(2)*xd);

a0=[1 1];%a0为参数的预测值，这个值取得好坏直接影响到最后结果，可能需要多次尝试。

[a,resnorm,residual]=lsqcurvefit(f,a0,xd,yd);

其中a 即为参数，resnorm 为总的误差的平方和，residual 为各点的误差。 计算得参数a 为0.2369，参数b 为0.3090。 （2）离散数据及拟合函数的图像如图41。

图41 发病率随年龄段的变化图

（3）下面42绘制离散点处的误差图，并计算均方误差。计算得到均方误差为1.9857。

图42 离散点处发病率的拟合值与实际值的误差

2模型未知的最小二乘拟合。

下图43为汽车价格随使用年限的离散点分布情况，观察猜测汽车的平均价格随使用年数的变化模型为

y =abx

其中a,b 为参数。

图43 汽车价格随使用年限的离散点分布情况

下面进行拟合参数a,b 的值。拟合代码如下。

>> x=[1 2 3 4 5 6 7 8 9 10]; >> yd=[2615 1943 1494 1087 765 538 484 290 226 204]; >>f=@(a,xd)a(1)*power(a(2),xd); >>a0=[1 0.9];

>> [a,resnorm,residual]=lsqcurvefit(f,a0,xd,yd);

拟合得到a=3544.5，b=0.7。拟合的均方误差为32.5。下图44是汽车价格随使用年限的离散点分布情况和拟合曲线图。

图44 汽车价格随使用年限的离散点分布情况和拟合曲线图

从拟合结果来看，汽车使用4.5年后价格在922.5±32.5这个范围内。

实验3.3研究最佳平方逼近多项式的收敛性质

x

f (x ) =e 取函数，在[-1，1]上以勒让德多项式为基函数，对于n =0, 1, , 10构造

ε(x ) =f (x ) -p n (x )

最佳平方逼近多项式p n (x ) ，令n ，将εn (x ) ~x 的曲线画在一

个图上。

令E n =max f (x ) -p n (x ) ，画出E n ~n 的曲线。做出E n ~n 之间的最小二乘曲线，

-1≤x ≤1

能否提出关于收敛性的猜测。 数值计算过程和结果 首先，利用自编的利用勒让德多项式为基函数进行连续函数最佳平方逼近的函数Legendrefit ，进行函数f 的0到10阶最佳平方逼近。然后计算

εn (x ) =f (x ) -p n (x )

，计算结果如下图45、46所示。

图45 0到5阶逼近函数与理论函数偏差的图像

图46 6到10阶逼近函数与理论函数偏差的图像

然后，利用自编的利用勒让德多项式为基函数进行连续函数最佳平方逼近的函数Legendrefit ，进行函数f 的0到10阶最佳平方逼近。然后计算

E n =max f (x ) -p n (x ) ，计算结果如下图47所示，拟合结果如下图48。

-1≤x ≤1

图47 En随n 的变化趋势

图48 En随n 的变化趋势和Efit(n)模型的拟合结果

拟合曲线模型为f(x) = a*exp(b*x)其中a = 5.723，b =-1.309，置信度大于95%。

综合ε(x)~x和En~n的计算结果可以看出，随着n 的增大，ε(x)变小，En 也不断。可以猜测，f(x)的以勒让德多项式为基的最佳平方逼近函数可以收敛到原函数f(x)。这与书上理论相符：勒让德多项式为连续函数的最佳平方逼近，如果函数足够光滑，此时勒让德多项式一致收敛到被逼近函数。

思考题：（病态）考虑将[0,1]30等分节点，用多项式y =1+x + +x 5生成数据，再用polyfit 求其3次、5次、10次、15次拟合多项式，并分析误差产生的原因。

计算过程及结果 利用polyfit 和polyval 函数对上面问题进行3次、5次、10次、15次拟合多项式拟合。表9为各阶拟合值与理论值偏差的无穷范数，图49为函数f(x)的理论值与各阶拟合值的图像显示。

表9各阶拟合值与理论值偏差的无穷范数

图49 函数f(x)的理论值与各阶拟合值

从计算结果来看，从3阶到5阶拟合精度提高，然而从5阶到10精度不再提高，

从10阶到15阶精度反而下降。原因如下，从理论上计算，用5阶及5阶以上多项式拟合本问题时，最佳拟合结果都应该为y =1+x + +x 5这个结果，及幂次高于5的项的系数都应该为0。然而，由于数值计算的舍入误差，当用高于5阶的多项式拟合时，5次以上的项的系数并不为0，如下表10所示。这种数值计算中高阶项的舍入误差导致了高阶拟合中的病态性。本问题中计算区间为(0,1)高次项的系数被不断缩小，所以影响不大，如果计算区间大于1，此时高次项系数不为0所引起的影响将很大。

表10 多项式10阶拟合时各项系数

项 X 10 X 9 X 8 X 7 X 6 X 5 X 4 X 3 X 2 X 1 1

系数

-0.[**************] 0.[**************] -0.[**************] 0.[**************] -0.[**************] 1.[**************] 0.[**************] 0.[**************] 1.[**************] 0.[**************] 1.[**************]