实验一-数据拟合

1实验一数据拟合一、实验目的和要求了解最小二乘法；了解数据拟合；了解Matlab做拟合的基本方式；了解拟合在数值计算中的应用。二、背景知识1.最小二乘法与数据拟合数据拟合问题的提法是，已知一组数据miyxii,,3,2,1,,，寻求一个拟合函数xfy，使将𝑥𝑖代入得到的f(𝑥𝑖)在某种准则下与𝑦𝑖最为接近（i=1,⋯,m）。在几何上，可以看作是找一条曲线xfy，与平面上已知的m个点(𝑥𝑖,𝑦𝑖)在某种准则下最为接近。常用的拟合函数有多项式函数、指数函数等。最常用的准则是使f(𝑥𝑖)与𝑦𝑖之差的平方和最小，称为最小二乘准则。根据最小二乘准则计算得到的待定系数的解称为最小二乘解，相应的方法称为最小二乘法。2.多项式拟合拟合函数为多项式函数mininiiinxaxaxaayaaaF12221010,,时，极小化上述函数得到的naaa,,10为最小二乘解。这个问题非常特殊，可以推导出解析解的表达式。利用多元函数极值的必要条件，得到线性方程组000122101221012210miininiinimiininiiimiininiiyxaxaxaaxyxaxaxaaxyxaxaxaa即2miinimiiimiinminiminiminiminimiimiiminimiiyxyxyaaaxxxxxxxx11110121111112111m，解此方程组可得到唯一一组naaa,,10的取值。3.matlab中的多项式拟合函数在matlab软件中，也可以借助专用的多项式拟合命令函数polyfit，它的基本使用格式为：函数：p=polyfit(x,y,n)功能：用n阶多项式拟合数据列(x,y)。说明：参数(x,y)是数据列的坐标，参数n是多项式的阶。返回值p是多项式从高阶到低阶的系数4.多项式求值通过数据拟合得到的只是多项式的系数，如果要对这个多项式求值，则调用函数polyval即可，这个函数的基本使用格式为：函数：polyval(p,x)功能：求多项式p在x处的函数值。x是数组的话返回函数值数组。说明：参数p是多项式的系数数组，从高阶到低阶。三、例题例1下表是纽约黑鲈的体重W与体长L之间的一组测量数据，试拟合出体重W与体长L之间的函数关系。表1黑鲈的体重W与体长L测量数据W(ozs)171723274149L(in)12.5012.6314.1314.5017.2517.75（1）所给数据之间关系的散点图代码：L=[12.5,12.63,14.13,14.50,17.25,17.75];W=[17,17,23,27,41,49];3plot(L,W,'r+')xlabel('体长');ylabel('体重');legend('体重与体长关系散点图')运行结果：图2-1。（2）选择3阶多项式，先数据拟合代码：p=polyfit(L,W,3)运行结果：0.42932-18.709274.49-1329.8。即：320.4293218.709274.491329.8WLLL（3）拟合曲线与原始数据对比图代码：plot(L,W,'r+',L,polyval(p,L))xlabel('体长');ylabel('体重');运行结果：图2-2。图2-1图2-2例2(斐波那契数列)某人养了一对兔，一个月后生育了一对小兔。假设小兔一个月后就可以长大成熟，而每对成熟的兔每月都将生育一对小兔，且假设兔子不会死亡。问：一年后共有多少对兔子？（答案：144对）这个问题，最早由意大利数学家斐波那契(Fibonacci)，于1202年在其著作《珠算原理》中提出。根据问题的假设，兔子的总数目是如下数列：1，1，2，3，5，8，13，21，34，55，89，144，233，…这个数列通常被称为斐波那契(Fibonacci)数列，我们会问：第n个月后，总共有多少对兔子？即Fibonacci数列的第n项是多少？这就需要我们探索Fibonacci数列的通项公式。根据问题的描述，我们知道第n+2个月后兔子的对数，等于第n+1个月后兔子的对数（表示原来就有的老兔子对数），加上第n个4月后兔子的对数（表示生育出来的新兔子对数）。这样就得到关于Fibonacci数列的一个递推公式：21nnnFFF有了这个递推公式，使用数学方法就能够得到这个数列的通项公式如下：{[(15)2][(15)2]}5nnnF这个公式是法国数学家比内(Binet)早在1843年发现的，称为比内公式。有了这个公式后，第n个月后兔子的对数，就是计算nF。请利用matlab软件强大的数据可视化功能，将Fibonacci数列显示成平面曲线的形式，更直观地观察其变化规律，并试着推导Fibonacci数列的通项公式。实验过程本试验将Fibonacci数列的有限项，看成是待处理的数据。首先利用matlab软件的可视化功能，将这些数据显示在平面坐标系中，观察其图形类似什么曲线，结论是：指数函数的曲线。进一步，利用指数函数与对数函数的互逆关系，将原有数据取对数，再观察其曲线形状是否类似直线，以验证原来的观察是否正确。通过观察到的目标函数，再利用matlab软件的数据拟合功能，得到Fibonacci数列通项公式的近似关系。最后，从近似关系出发，推导出来Fibonacci数列的通项公式。1.观察数列的大概函数关系为了研究Fibonacci数列的变化规律，我们取此数列的前有限项来观察。利用Matlab软件的数据可视化功能，将这些数据显示在平面坐标系中，观察其中蕴涵的函数关系。为了调用方便，设计了一个函数，具体代码如下：functionfib1(n)%显示Fibonacci数列前n项fn=[1,1];%将数列的前两项放到数组fn中fori=3:nfn=[fn,fn(i-2)+fn(i-1)];%将第i项添加到数组fn中endplot(fn)%画曲线将这个文件保存在Matlab的当前目录下，文件名fib1.m。在这个函数里，n是参数，即自变量，在调用时根据需要赋具体的值，调用方式与系统内部函数5相同。（1）选择n=30，调用上述函数画图：代码：fib1(30);legend('n=30');运行结果：图2-1-5。（2）选择n=50，调用上述函数画图：代码：fib1(50);legend('n=50');运行结果：图2-1-6。图2-1-5n=30图2-1-6n=50（3）选择n=500，调用上述函数画图：代码：fib1(500);legend('n=500');运行结果：图2-1-7。（4）选择n=1000，调用上述函数画图：代码：fib1(1000);legend('n=1000');运行结果：图2-1-8。图2-1-7n=500图2-1-8n=1000（5）初步结论6经过反复观察，觉得曲线上升的速率非常快，远远快过n增大的速率。而常见的函数中，上升速率最快的应该是指数函数，故初步猜测Fibonacci数列通项公式的函数关系，应该近似为指数函数。读者可以取更多的参数n的值，反复观察寻找规律。2.进一步验证上述结论经过上一步的观察，觉得这些数据应该是指数函数的形式。为了进一步验证这个结论是否正确，可以利用指数函数与对数函数的互逆关系。如果这些数据确实是指数函数的形式，则经过取对数后应该是一个线性关系，即一阶多项式，从图形上看应该象一条直线。因此，再利用Matlab软件的数据可视化功能，将这些数据取对数后显示在平面坐标系中，观察它是否象一条直线。为了调用方便，设计了一个函数，具体代码如下：functionfib2(n)%显示取对数后的前n项fn=[1,1];%将数列的前两项放到数组fn中fori=3:nfn=[fn,fn(i-2)+fn(i-1)];%将第i项添加到数组fn中endfn=log(fn);%将原来的数据取对数plot(fn)%画曲线将这个文件保存为名fib2.m。（1）选择n=30，调用上述函数画图：代码：fib2(30);legend('n=30');运行结果：（2）选择n=50，调用上述函数画图：代码：fib2(50);legend('n=50');运行结果：（3）选择n=500，调用上述函数画图：7代码：fib2(500);legend('n=500');运行结果：图2-1-11。（4）选择n=1000，调用上述函数画图：代码：fib2(1000);legend('n=1000');运行结果：图2-1-12。图2-1-11n=500图2-1-12n=1000（5）确切的结论经过反复观察，以上四条曲线的确非常直。故n比较大以后，Fibonacci数列取对数后的函数曲线，可以近似看成直线，因而验证了上一步的猜测是基本正确的，Fibonacci数列的函数关系，应该就是为指数函数。3.获取数列的近似函数关系经过以上第一步的观察，确定Fibonacci数列的数据是指数函数的关系，第二步验证了第一步得到的结论，因此我们认为Fibonacci数列的数据关系就是指数函数，取对数后就是线性函数，即一阶多项式。利用Matlab软件的数据拟合功能，通过取对数后的数据，用一阶多项式拟合出它的函数关系式，可以得到Fibonacci数列通项公式的一个近似表达式。为了调用方便，设计了一个函数，具体代码如下：functiony=fib3(n)%针对取对数后的数据，做线性拟合fn=[1,1];%将数列的前两项放到数组fn中fori=3:nfn=[fn,fn(i-2)+fn(i-1)];%将第i项添加到数组fn中endxn=1:n;%定义横坐标8fn=log(fn);%将原来的数据取对数，竖坐标y=polyfit(xn,fn,1);%一阶多项式拟合将这个文件保存为fib3.m。在这个函数里，y是因变量，用于将拟合结果传到函数外。（1）选择n=30，调用上述函数做拟合：代码：p1=fib3(30)运行结果：p1=0.4799-0.7768。结论：取前30项做拟合，得到：log()0.7768+0.4799nnF（2）选择n=50，调用上述函数做拟合：代码：p2=fib3(50)运行结果：p2=0.4807-0.7881。结论：取前50项做拟合，得到：log()0.7881+0.4807nnF（3）选择n=500，调用上述函数做拟合：代码：p3=fib3(500)运行结果：p3=0.4812-0.8031。结论：取前500项做拟合，得到：log()0.8031+0.4812nnF（4）选择n=1000，调用上述函数做拟合：代码：p4=fib3(1000)运行结果：p4=0.4812-0.8039。结论：取前1000项做拟合，得到：log()0.8039+0.4812nnF（5）最终结论从上面四次拟合的结果来看，拟合系数基本趋于稳定，故我们选定第四个作为Fibonacci数列取对数后的近似函数关系，即：log()0.8039+0.4812nnF对上面的表达式，做指数函数运算来还原：代码：p5=exp(p4)运行结果：p5=1.61800.4476。结论：Fibonacci数列的近似函数关系是：n0.44761.6180nF4.观察拟合数据与原始数据的吻合程度9经过第三步的拟合，得到了Fibonacci数列近似的通项公式，为了观察其吻合程度，我们将Fibonacci数列的拟合数据与原始数据的图形显示出来，进行对比观察。为了调用方便，设计了一个函数，具体代码如下：functiony=fib4(n)%显示拟合数据与原始数据的前n项fn1=[];%装拟合数据的数组fori=1:nfn1=[fn1,0.4476*1.618^i];%将第i项添加到数组endfn2