毕业论文 曲线拟合.doc

目 录摘要1前言21 问题提出32 插值介绍42.1拉格朗日公式求解42.1.1 算法分析52.1.2 程序设计52.1.3 计算结果82.1.4 计算量分析82.1.5 图形绘制与分析82.2 逐次线性插值公式求解92.2.1 算法分析102.2.2 程序设计102.2.3 计算结果122.2.4 计算量分析122.3 牛顿插值公式求解132.3.1 算法分析132.3.2 程序设计142.3.3 计算结果152.3.4 计算量分析162.3.5 图形绘制与分析163 拟合介绍173.1 最小二乘拟合求解173.1.1 算法分析183.1.2 程序设计183.1.3 计算结果213.1.4 计算量分析213.1.5 图形绘制与分析23结论25参考文献26致谢27附录28附1 MATLAB绘图源程序28附2 拉格朗日插值求解系数的具体流程图29摘 要曲线拟合是应用数学里数值计算方法中的一个小的分支，本文先给出湖水污染中氯磷浓度的几对离散数据，然后用不同的方法解决这个问题，具体是大学期间所学的拉格朗日插值、逐次线性插值、牛顿插值与最小二乘拟合法。同时，对各种方法进行了系统的算法分析、程序设计，并对各种方法的计算量进行了分析，绘制了不同方法下的各种图形，针对计算量与图形效果，讨论了各种方法的优点、缺点，从而得出了什么情况下选用何种方法更好。关键词：插值拟合；数值计算方法；应用数学；程序设计。AbstractCurve approach is an important part of the computing technology of number value in applied mathematics. In this thesis the author firstly gives a few pairs of dispersed data of the phosphorus density of chlorine in the pollutions of the lake, then work out the problem with different ways, they are Lagrange Interpolation, linear Interpolation, Newton Interpolation and minimum to be approach, which I learned in the university. At the same time, I Proceeded the algorithm analysis, procedure design by system amount of these different methods, analysis the calculating amount of these different methods and design all kinds of graphics under the different methods, Aim at calculating amount and The result of graphics, discusses the advantages and the disadvantages of the different methods, and then work out which method is the best way under the certain circumstance.Key word: Interpolation approach; Value computational method; applied mathematics; Procedure design.前 言在解决一个现实工程的问题时，由于计算机的高速、大容量、多功能、又为现代科学技术的发展提供了最优、最快的新途径，一般可按四个阶段进行。一、工程问题数学化（数学建模），二、数学问题数值化（算法与分析），三、数值问题机器化（程序设计），四、科学试验。任何一个工程问题，任何一项能产生社会和经济效益的科研课题，必须经过上述4个阶段，第一阶段是根本，模型是否反映工程实际问题的需要，取决于当代高级科技人员的理论水平；第二阶段是桥梁，所构造的算法是否能真正取代原数学模型，它的有效性和可实现性，取决于从事计算数学的研究人员的理论水平和创新能力；第三阶段是检验前两阶段工作的有效性和可行性的方法和手段；第四阶段检验该项研究成果是否有实际应用价值，同时也检验了前三阶段研究工作的可靠性，也是能否见到经济效益和社会效益的关键。在科学研究和其他许多实际问题中，常常有函数不便于处理或计算的情形。有时候函数关系没有明显的解析表达式，需要根据实验观测或其他方法来确定与自变量的某些值相对应的函数值；有时函数虽有解析表达式，但是使用很不方便。因此，希望对这些问题中的函数建立一个简单的便于计算处理的近似表达式，即用一个简单的函数来近似代替这些不便处理的函数，用近似数代替准确数。我们把观测得来的离散数据重新构造函数的方法叫做离散问题解析化。离散问题解析化，一般有两种途径：其一，插值法，包括多项式插值，分段低次多项式插值（样条插值等），逐步线性插值等；其二，逼近法，最小二乘逼近，最优一致逼近，小波逼近等。本文介绍大学期间所学的几种方法，对同一个问题用不同方法进行计算，并对结果进行分析，方法的选择是本文要讨论的内容。1 问题提出河流水质的主要污染指标为总磷、氨、氮、汞、挥发酚和铅，影响湖泊水质的主要指标是总氮、总磷及耗氧有机物。湖泊中氮磷及其耗氧有机物含量高处在越高的富营养状态。湖水中氯和磷的浓度之间有一定的关系，下面表1是一组不同湖水中的测量值。 表1湖泊磷浓度（p）mg/m3氯浓度（c）mg/m312345674.58.05.539.019.517.521.00.82.01.211.04.43.85.5当氯浓度为15mg/m3时，磷的浓度是多少？试建立用磷浓度预报氯浓度的经验公式，并试分析河流中氯磷浓度关系。2 插值介绍设函数在区间上给出了一系列的函数值或者给出一张函数表 表2x这里 b选择一个函数，满足作为函数的近似表达式，称这样的问题为插值问题，满足关系式的称为的插值函数，称为被插值函数，点称为插值节点。区间称为插值区间。这样在给点计算插值函数的值作为被插值函数的近似值，这一过程称为插值，点称为插值点。所以说所谓插值就是根据已知点的函数值求其余点的函数值，也不是依据被插值函数给出的函数表“插出”所要点的函数值。由于插值函数的选择不同，就产生不同类型的插值，若为代数多项式就是代数插值，若为三角多项式就称为三角多项式插值，若为有理函数就称为有理函数插值等，选用不同的插值函数，近似的效率也不同。2.1 拉格朗日公式求解拉格朗日插值是最常用的插值公式，以下仅介绍其公式，构造方法及证明请查阅相关资料，它的一次插值构造方法是：已知函数y=f(x)在节点x0,x1有函数值y0=f(x0)，y1=f(x1)。如图1所示， y=P(x)x1x0y1y0y=f(x) 图1这种插值叫线性插值，P(x)就是由(x0,y0)与(x1,y1)构造出来的插值函数，这个新构造出来的函数就是为了替代原来我们不一定知道的f(x)或者难以用简单曲线描绘的f(x)。为了便于推广，记我们就有了，同理有n次拉格朗日插值多项式则有2.1.1 算法分析（采用自然语言描述）1） 输入插值节点的总数n，并保证2n15；2） 输入每个插值点，在确认输入无误后，输入1确认输入的插值点正确；3） 利用辅助数组xx 、xy 计算出插值节点函数的系数数组xi ，并且将计算后的插值函数输出；4） 输入所求点的x值，确定它在输入的插值点的范围内，利用循环调用函数，函数作用为计算；5） 输出计算结果。2.1.2 程序设计float max(float array ,int m) /*找出输入点x的最大值*/float p=array0; int i; for(i=0;im;i+) if(parrayi) p=arrayi; return(p); float min(float array ,int m) /*找出输入点x的最小值*/float p=array0; int i; for(i=0;iarrayi) p=arrayi; return(p); float lk(float array ,int j,int m,float a) /*求出lk值*/ float p=1; int i; for(i=0;im;i+) if(j!=i)p=p*(a-arrayi)/(arrayj-arrayi); return(p); #define N 15main()int n,i,j,k; /*n个插值节点*/ float xN,yN,xiN+1,xxN+1,xyN+1; float a,x1,y1=0; /*所求节点XY值*/ doprintf(please input n,2nN)printf(too bign); else if(n2)printf(too smalln); while (nN); /*保证值不太多或太少*/ dofor(i=0;in;i+) printf(please input x%d,y%dn,i,i); scanf(%f,%f &xi,&yi); /*输入已知点*/ for(i=0;in;i+) printf(%f,%fn,xi,yi); /*输出已知点*/ printf(OK? right input 1, eles input other number); scanf(%d,&i); while(i!=1); /*确认输入各点并无误*/ for(j=0;jN+1;j+) /*系数数组清零*/ xij=0; for(k=0;kn;k+) for(xx0=1,xy0=0,j=1;jN+1;j+) /*辅助数组清零*/ xxj=0; xyj=0; for(i=0;in;i+) /*计算x-xi的系数*/ for(j=0;jN;j+) if(i!=k) xyj+1=xxj*(-xi); for(j=0;jN+1;j+) if(i!=k) xxj=xxj+xyj; for(a=1,i=0;in;i+) /*计算1/xk-xi的值*/ if(i!=k) a=a*1/(xk-xi); a=a*yk; /*计算yk/xk-xi的值*/ for(j=0;jN+1;j+) xxj=xxj*a; for(j=0;jN+1;j+) /*计算插值多项式系数*/ xij=xij+xxj; printf(P(x)=); for(j=0;jn;j+) /*控制和输出系数项*/ printf(%.2fx%d,xij,n-j-1); if(j!=n-1) printf(+); doprintf(nplease input x, min(x)x=max(x,n)|x1=min(x,n); /*保证所求点在输入点范围内*/ for(i=0;in;i+) y1=y1+lk(x,i,n,x1)*yi; /*调用lk子函数求所求点值*/ printf(x=%f,y=%fn,x1,y1); /*输出所求点拉格朗日插值*/ getch(); 2.1.3 计算结果在同一点即时y点的值，按照插值节点就所求点近的原则依次输入了，就拉格朗日不同的次数对求值，以下是利用程序得到的结果： 2.1.4 计算量分析在计算时先进行加减法计算量为。而乘除法计算量为；再计算时，加减法的计算量为，乘除法的计算量为；由于计算机作一次乘除法所需时远大于作一次加减法需的时间，因此估计一个算法时，一般只需考虑乘除法次数；拉格朗日插值公式计算量为；节点每增加一点，公式计算复杂度增加的乘除计算量为，n为增加一点后的总结点数。2.1.5 图形绘制与分析二点一次插值前面六个点均均匀分布在线的附近，最后一点偏离插值函数较远；三点二次插值插值点大致分布在曲线附近，构造的效果应该与现实情况更接近一些；三次、四次、五次构造函数前五六点大致接近测量值，但与第七点偏差很大，误差的原因是：高次插值函数系数精确度不够高，因为电脑计算精度的原因，只要将上面的程序显示精度调高我们就可以得到更精确的图。图2为数学软件MATLAB绘制的图形：图22.2 逐次线性插值公式求解拉格朗日插值多项式计算函数值，精度不够增加节点时，因基函数和第一个节点有关，原来的数据都不能利用，要重新计算，为使计算有“承袭性”，埃特金Aitken算法和内维尔Neville算法解决了这个问题。已知三个节点及其函数值,先过其中两个点,作一次插值；再过两个点,作一次插值；当给定，都是函数值，因此将,也看成是两个点，再进行线性插值，有，将和代入并化简可得三点二次抛物插值。同理多次插值当时，得到，当时，得到，当时，得到。反复执行这一算式，可以逐步构造出如下插值顺序，按这个顺序可求出某点的插值。埃特金算法的计算顺序： 内维尔算法的计算顺序：2.2.1 算法分析1） 输入插值节点的总数n，并保证2n10；2） 输入每个插值点，将赋给第一列赋给第二列在确认输入无误后，输入1确认输入的插值点正确；3） 输入所求点的x值，确定它在输入的插值点的范围内；4） 埃特金算法：计算第三列的值，接着计算第四列以后的值，输出埃特金算法的矩阵；5） 内维尔算法：计算第三列的值，接着计算第四列以后的值，输出内维尔算法的矩阵。2.2.2 程序设计#define N 10float max(float array N,int j) /*返回输入节点x最大值*/ float p=array00; int i; for(i=1;ij;i+) if(parrayi0) p=arrayi0; return(p); float min(float array N,int j) /*返回输入节点x最小值*/ float p=array00; int i; for(i=1;iarrayi0) p=arrayi0; return(p); main()int n,i,j=0; float aNN; float x; doprintf(please input n,2nN)printf(too bign); else if(n2)printf(too smalln); while (nN); /*保证值不太多或太少*/ dofor(i=0;in;i+) printf(please input x%d,y%dn,i,i); scanf(%f,%f,&aij,&aij+1); /*输入各节点数值*/ for(i=0;in;i+) printf(%f,%fn,aij,aij+1); /*输出各节点数值*/ printf(OK? right, input 1,eles input other number:); scanf(%d,&i); while(i!=1); /*确定输入无误*/ doprintf(please input x, min xx=max(a,n)|x=min(a,n); /*确定所求点在输入点范围内*/ for(i=1;in;i+) /*计算第三列数值*/ ai2=(x-ai0)*a01/(a00-ai0)+(x-a00)*ai1/(ai0-a00); for(j=3;jn+1;j+) /*计算第四列以后的数值*/ for(i=j-1;in;i+) aij=(x-ai0)*aj-2j-1/(aj-20-ai0)+(x-aj-20)*aij-1/(ai0-aj-20); printf(Aitkenn); for(i=0;in;i+) for(j=0;ji+2;j+) printf(%12.5f,aij); /*输出埃特金算法形式*/ printf(n); for(i=1;in;i+) /*计算第三列数值*/ ai2=(x-ai0)*ai-11/(ai-10-ai0)+(x-ai-10)*ai1/(ai0-ai-10); for(j=3;jn+1;j+) /*计算第四列以后的数值*/ for(i=j-1;in;i+) aij=(x-ai0)*ai-1j-1/(ai+1-j0-ai0)+(x-ai+1-j0)*aij-1/(ai0-ai+1-j0); printf(Nevillen); for(i=0;in;i+) for(j=0;j1) printf(x-%.2f),arrayj-20); xxi(j-1,array); float max(float array N,int j) /*返回输入节点x最大值*/ float p=array00; int i; for(i=1;ij;i+) if(parrayi0)p=arrayi0; return(p); float min(float array N,int j) /*返回输入节点x最小值*/ float p=array00; int i; for(i=1;iarrayi0)p=arrayi0; return(p); main()int n,i,j=0; float aNN; float x,y; doprintf(please input n,1nN)printf(too bign); else if(n2)printf(too smalln); while (nN); /*保证值不太多或太少*/ dofor(i=0;in;i+) printf(please input x%d,y%dn,i,i); scanf(%f,%f,&aij,&aij+1); /*输入各节点数值*/ for(i=0;in;i+) printf(%f,%fn,aij,aij+1); /*输出各节点数值*/ printf(OK? right, input 1,eles input other number:); scanf(%d,&i); while(i!=1); /*确定输入无误*/ for(i=1;in;i+) /*计算第三列数值*/ ai2=(ai1-ai-11)/(ai0-ai-10); for(j=3;jn+1;j+) /*计算第四列以后的数值*/ for(i=j-1;in;i+) aij=(aij-1-ai-1j-1)/(ai0-ai+1-j0); printf(Newtonn); /*输出牛顿差商表*/ for(i=0;in;i+) for(j=0;ji+2;j+) printf(%10.5f,aij); printf(n); printf(P(x)=); for(i=1;in+1;i+) /*输出牛顿插值多项式*/ printf(%.6f,ai-1i); xxi(i,a); if(i!=n)printf(+); doprintf(nplease input x, min xx=max(a,n)|x=min(a,n); /*确定所求点在输入点范围内*/ for(y=0,i=1;in+1;i+) /*计算牛顿插值*/ y=y+ai-1i*xzi(i,a,x); printf(y=%f,y); getch( ); /*暂停,方便观察*/2.3.3 计算结果；2.3.4 计算量分析计算差商表的计算量与逐次相比，计算量少一半，共进行次乘除运算，在计算时，所执行的乘除法次数为；综上，牛顿插值所做的乘除法计算量为；每增加一点计算量增加。2.3.5 图形绘制与分析图4牛顿插值与拉格朗日插值情况大体相同。3 拟合介绍似合又叫逼近法，逼近法与插值法类似，构造一分析性能好，又便于应用的函数去取代，或离散点，不同之处在于逼近函数不一定过离散点，这些离散点尽可能在的附近，附近的程度取决于各种逼迫标准，而这些标准一般是根据实际问题的需要提出的。3.1最小二乘拟合求解已知有离散点，要求在某特定函数类|（例如多项式）中找出一个函数作为的近似函数，使得在上的误差按某种度量标准为最小，这就是拟合问题，也称为出线拟合。记向量，要求残差按某种度量标准最小，即要求向量的某种范数最小。例如，可要求1-范数或-范数即最小。可是计算不方便，因此通常要求的2-范数即为最小。这种要求误差平方和最小的拟合的最小二乘法。直线拟合设已知数据点，分布大致为一条直线。作拟合直线，该直线不是通过所有的，而是使残差平方和为最小，即使为最小。确定直线参数的方法是使达到极限的参数应满足即成立 解得值，则有为直线拟合的函数。多项式拟合有时所给数据点用直线拟合并不合适，这时可考虑用多项式拟合。对于给定的一组数据寻求作次多项式，使总误差为最小。正则方程组解得，则有为多项式拟合的函数。3.1.1 算法分析1） 输入拟合所用节点的总数n，并保证2n15；2） 输入每个拟合点，构造正则方程组的系数矩阵，输出构成后的系数矩阵；3） 用高斯约当消去法计算，输出计算后的系数矩阵；4） 计算出的值，并将拟合后的函数写出；5） 输入所求点的x值；6） 计算所求点的y值，并输出。3.1.2 程序设计#define N 15#define M 10float zhz(float array,int n,int i) float p=1; if(n=0) p=1; else p=p*arrayi*zhz(array,n-1,i); return(p); float zhx(float arrax,int n,int m) float p=0; int i; for(i=0;im;i+) p=p+zhz(arrax,n,i); return(p);float zhy(float arrax,float array,int n,int m)float p=0; int i; for(i=0;im;i+) p=p+arrayi*zhz(arrax,n,i); return(p); float zhi(int j,float x)float p=1; if(j=0) p=1; else if(j!=0) p=x*zhi(j-1,x); return(p);main()int n,m,i,j,h,a,b; float xN,yN,zM+2M+2; float x1,y1; /*所求节点XY值*/ doprintf(please input n,2nN)printf(too bign); else if(n2)printf(too smalln); while (nN); /*保证值不太多或太少*/ dofor(i=0;in;i+) printf(please input x%d,y%dn,i,i); scanf(%f,%f,&xi,&yi); /*输入已知点*/ for(i=0;in;i+) printf(%f,%fn,xi,yi); /*输出已知点*/ printf(OK? right input 1, eles input other number); scanf(%d,&i); while(i!=1); /*确认输入各点并无误*/ doprintf(please input m 0m%d,M); /*输入多项式拟合的次数*/ scanf(%d,&m); if(mM) printf(Erraor,input again!); while(mM); /*确定次数不太高*/ for(j=0;jm+1;j+) /*给线性方程组左边矩阵项赋值*/ for(i=0;im+1;i+) if(i=0&j=0) zij=n; else zij=zhx(x,i+j,n); for(i=0;im+1;i+) /*给线性方程组右边项赋值*/ zim+1=zhy(x,y,i,n); printf(hang lie shi wein); /*输出赋值后的多项式系数项*/ for(i=0;im+1;i+) for(j=0;jm+2;j+) printf(%12.2f,zij); printf(n); for(h=0;hm+1;h+) /*用高斯消去法计算矩阵*/ for(i=0;im+1;i+) if(i!=h) for(j=0;jm+2;j+) zm+1j=zhj*zih/zhh; for(j=0;jm+2;j+) zij=zij-zm+1j; printf(ji suan hou de hang lie shin); /*输出计算后的多项式系数项*/ for(i=0;im+1;i+) for(j=0;jm+2;j+) printf(%12.2f,zij); printf(n); for(i=0;im+1;i+) /*计算a0到an的值*/ zm+1i=zim+1/zii; printf(P(x)=); /*输出生成后的n次函数*/ for(i=0;im+1;i+) printf(%.10fx%d,zm+1i,i); if(i!=m) printf(+); printf(nplease input xn); /*输入所求点对应的y值*/ scanf(%f,&x1); for(y1=0,i=0;im+1;i+) y1=y1+zm+1i*zhi(i,x1); printf(x=%f,y=%.6fn,x1,y1); getch(); 3.1.3 计算结果在插值中随着插值节点的增加，次数也随之增加。先用这种方法观察拟合的函数。；在拟合中，我们可以用多个点来拟合，以下是用七点来进行拟合的；3.1.4 计算量分析在生成正则方程组中，设输入个拟合离散点，生成的拟合曲线为次，先对的系数进行计算和赋值，再用高斯约当消去法变换，使之成为对角矩阵，再计算的值。在对的系数进行计算时，先计算右边一列时的乘法计算量为0；的乘法计算量为；的乘法计算量为；总计右边一列计算量为；将左边的矩阵按从右上到左下角这种划分方式来划分列.，左上角为第一列，依次推算右下角为列。第1列，乘除计算量为0；第2列，乘除计算量也为0；第3列，乘除计算量为；第m+1列，乘除法计算量为；第m+2列，乘除法计算量为；第2m+1列，乘除法计算量为；计算左边矩阵，总计计算量为：在这里高斯约当消去法实际运行的乘除次数为；教材9上的资料是；求的计算量为；将所求点代入后，计算量为；最小二乘曲线拟合计算量由赋值到行列式，计算系数和代值求解构成，总计计算量为：当时构成的类拟于插值的函数（两点一次、三点二次），其计算量为。3.1.5 图形绘制与分析图5图5是类拟与插值函数的图形（二点一次拟合，三点二次拟合），构造而成的效果大致与插值函数相同，三点二次效果最好。而图6是最小二乘拟合由MATLAB绘制的图形，从图上看到，二次、三次、四次的逼近效果都不错，五次逼近函数在附近达到了最高峰，六次函数效果相比前面低次要差许多，但它相比前面的插值函数要好，它的整个函数图像均落在图框内部。图6结论就本文中的例题来看，从计算量的角度分析，拉格朗日插值和逐次线性插值计算量几乎相同，而牛顿插值的计算量略大于以上两种方法，它们插值的效果大致相同，二次插值曲线为最好的情况，插值构造的曲线均在这个区域内。构造函数次数选择和节点选择同样重要，高次函数构造所要求的系数精度相当高，而且随着次数的增加远离所求点的数值偏差越大（越远离的点）。而曲线拟合的情况有些不同，类似于插值所构造的拟合函数（两点一次，三点二次，）跟插值函数大致相似，三次以上逼近的效果不理想，而最小二乘拟合所构造的函数，二次、三次、四次效果都比较理想，但其计算量却比插值计算大许多。综上所述，通过对湖水氯磷浓度估测和计算，得出当我们只要求一个或几个相当接近的点的数值时，我们只需构造简单的插值函数来进行计算；如果要求对所涉及到的对象有系统的了解时，拟合函数比插值函数计算量大