曲线拟合的数值计算方法实验

《数值方法》实验报告1曲线拟合的数值计算方法实验郑发进2012042020022【摘要】实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合（curvefitting）是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。曲线直线化是曲线拟合的重要手段之一。对于某些非线性的资料可以通过简单的变量变换使之直线化，这样就可以按最小二乘法原理求出变换后变量的直线方程，在实际工作中常利用此直线方程绘制资料的标准工作曲线，同时根据需要可将此直线方程还原为曲线方程，实现对资料的曲线拟合。常用的曲线拟合有最小二乘法拟合、幂函数拟合、对数函数拟合、线性插值、三次样条插值、端点约束。关键词曲线拟合、最小二乘法拟合、幂函数拟合、对数函数拟合、线性插值、三次样条插值、端点约束一、实验目的1.掌握曲线拟合方式及其常用函数指数函数、幂函数、对数函数的拟合。2.掌握最小二乘法、线性插值、三次样条插值、端点约束等。3.掌握实现曲线拟合的编程技巧。二、实验原理1.曲线拟合曲线拟合是平面上离散点组所表示的坐标之间的函数关系的一种数据处理方法。用解析表达式逼近离散数据的一种方法。在科学实验或社会活动中，通过实验或观测得到量x与y的一组数据对（Xi，Yi)（i=1，2，...m），其中各Xi是彼此不同的。人们希望用一类与数据的背景材料规律相适应的解析表达式，y=f(x，c）来反映量x与y之间的依赖关系，即在一定意义下“最佳”地逼近或拟合已知数据。f(x，c)常称作拟合模型，式中c=(c1，c2，…cn)是一些待定参数。《数值方法》实验报告2当c在f中线性出现时，称为线性模型，否则称为非线性模型。有许多衡量拟合优度的标准，最常用的一种做法是选择参数c使得拟合模型与实际观测值在各点的残差(或离差)，c）-f（fyekkk的加权平方和达到最小，此时所求曲线称作在加权最小二乘意义下对数据的拟合曲线。有许多求解拟合曲线的成功方法，对于线性模型一般通过建立和求解方程组来确定参数，从而求得拟合曲线。至于非线性模型，则要借助求解非线性方程组或用最优化方法求得所需参数才能得到拟合曲线，有时称之为非线性最小二乘拟合。曲线拟合：贝塞尔曲线与路径转化时的误差。值越大，误差越大；值越小，越精确。2.最小二乘法拟合：最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。函数曲线为：y=Ax+B其中系数满足下列的正规方程：NkkkNkkNkkyxBxAx1112NkkNkkyNBAx113.幂函数拟合：函数曲线为：设Νkk，kyx1有N个点，其中横坐标是确定的。最小二乘幂函数拟合曲线的系数A为：)/()(121NkMkNkkMkxyxA、4.对数函数拟合：《数值方法》实验报告3对数函数（lograrithmicfunction）的标准式形式为)0(lnaYXXbb0时，Y随X增大而增大，先快后慢；b0时，Y随X增大而减少，先快后慢，见图12.4(c)、(d)。当以Y和lnX绘制的散点图呈直线趋势时，可考虑采用对数函数描述Y与X之间的非线性关系，式中的b和a分别为斜率和截距。更一般的对数函数Y=a+bln(X+k)式中k为一常量，往往未知。(a)lnY=lna+bX(b)lnY=lna-bX(c)Y=a+blnX(d)Y=a-blnX5.线性插值：在代数插值中,为了提高插值多项式对函数的逼近程度一般是增加节点的个数，即提高多项式的次数，但这样做往往不能达到预想的效果。如下图所示：f(x)=1/(1+x2)如果在区间[-5，5]上取7个等距节点:xk=5*(k/3-1)(k=0,1,2,...,6),由lagrange插值公式可得到f(x)的次L7(x)。如图所示：L7(x)仅在区间的中部能较好的逼近函数f(x),在其它部位差异较大，而且越接近端点，逼近效果越差。可以证明，当节点无限加密时，Ln(x)也只能在很小的范围内收敛，这一现象称为Runge现象。它表明通过增加节点来提高逼近程度是不适宜的，因而不采用高次多项式插值。如果我们把以上的点用直线连接起来，显然比L7(x)要更逼近f(x)。这就是分段线性插值。而在实际应用中通常采用分段低次插值来提高近似程度，比如可用分段线性插值或分段三次埃尔米特插值来逼近已知函数，但它们的总体光滑性较差。为了克服这一缺点，一种全局化的分段插值方法——三次样条插值成为比较理想的工具。《数值方法》实验报告46.三次样条插值：设Nkkk，yx0有N+1个点，其中bx...xxxaN210。如果存在N个三次多项式x）Sk（，系数为3210k，k，k，k，，S，S，SS满足如下性质：33,22,1,0,)()()()()(kkkkkkkkxxsxxsxxssxSxS1,....,1,0],,[1Nkxxxkk)('')('')(')(')()()(111111111kkkkkkkkkkkxSxSxSxSxSxSyxS2,....1,02,....1,02,....1,0,....1,0NkNkNkNk则成函数S(x)为三次样条函数。7.端点约束：紧压样条：存在唯一的三次样条曲线，其一阶导数的边界条件是：NdbSdaS)(',)('0natural样条：存在唯一的三次样条曲线，它的自由边界条件是：0)('',0)(''bSaS外推样条：存在唯一的三次样条曲线，其中通过对点x1和x2进行外推得到aS，同时通过对点X(n-1)和X(N-2)进行外推得到bS。端点曲率调整：存在唯一的三次样条曲线，其中二阶导数的边界条件aS和bS是确定的。抛物线终结样条：存在唯一的三次样条曲线，其中二阶在区间[X0，X1]内0xS）（，而在[Xn-1,Xn]内0xS）（。《数值方法》实验报告5三、实验内容1.P2021胡克定律指出F=kx，其中F是拉伸弹簧的拉力（单位为盎司），x为拉伸的长度（单位为英寸）。根据下列试验数据，求解拉伸常量k的近似值。（a）（b）2.P2151洛杉矶（美国城市）郊区11月8日的温度记录入下表所示，其中共有24个数据点。（a）根据例5.5中的处理过程（使用fmins命令），对给定的数据求解最小二乘曲线EDxCBxAxf)sin()cos()(。（b）求。)(2fE（c）在同一坐标系下画出这些点集和(a)得出的最小二乘曲线。xkFk0.23.60.47.30.610.90.814.51.018.2XkFk0.25.30.410.60.615.90.821.21.026.4时间，p.m.温度时间，a.m.温度1661582662583653584644585635576636577627578618589609601060106411591167午夜58正午68《数值方法》实验报告63.P2291一个轿车在时间Tk时经过的距离dk，如下表所示。使用程序5.3，并根据一阶导数边界条件98)8(',0)0('SS,求这些数据的三次紧压样条插值。时间，tk02468距离，dk0401603004804.P2385美国洛杉矶郊区11月9日的温度（华氏温度）如表5.10所示。采用24小时制。（a）求三角多项式)(7xT（b）在同一坐标系下，画出图)(7xT和24个数据点。（c）使用本地的温度情况重新求解问题(a)和问题(b)。时间,p.m温度时间,a.m温度1661582662583653584644585635576636577627578618589609601060106411591167午夜58正午685.P2461《数值方法》实验报告7编写Matlab程序，生成并绘制组合贝塞尔曲线。利用该程序生成和绘制过3个控制点集{（0,0），（1,2），（1,1），（3,0）}，{（3,0），（4，-1），（5，-2），（6,1），（7,0）}，{（7,0）,(4,-3),(2,-1),(0,0)}的贝塞尔曲线。四、实验结果及分析1.P2021实验描述：由题意可知，此题需要用最小二乘法进行计算，因为已知函数的5个插点并且知道它们的x,y的值。且函数的表达式为F=kx，所以只需用方程中NkkNkkyAx11)(便可计算出k的数值。定义一个动态数组a,用来依次存取x和y的插值。其中x,y的插值通过键盘手动输入并赋予给a中的元素。然后通过相应的求和和基本运算便可以得到相应插值下的最小二乘法的函数表达式。（其中k保留小数点后4位）具体函数运行效果如下：（a）请在此输入x的各插值0.20.40.60.81.0请在此输入y的各插值3.67.310.914.58.2此函数的最小二乘法曲线表达式为Y=14.8333*x请按任意键继续。。。（b）请在此输入x的各插值0.20.40.60.81.0请在此输入y的各插值5.310.615.921.226.4此函数的最小二乘法曲线表达式为Y=26.4667*x请按任意键继续。。。结果分析：易得，最小二乘法多项式计算可以很好的做出较为精确的最小二乘法拟合曲线，并且程序通用性高，只要输入相应的插值便可以计算出结果，结果系数的小数点有效位同时也可以自己确定。《数值方法》实验报告82.P2151实验描述：给出的最小二乘曲线表达式为：EDxCBxAxf)sin()cos()(其中变量有5个，而给出的数据点有24个，在C语言中可以用牛顿-拉夫森算法迭代计算分别得出5个变量的值，但是方法繁琐，且迭代计算量庞大，因此这里采用Matlab进行相关的计算，调用fminsearch函数，求得当5个参量都为1附近时候多项式的最小值，此时便可求出此5个参变量的值.然后继续通过Matlab,将得到的公式以及各点，用plot函数，便可以求得。实验结果：运行matlab结果如下：fminsearch('fiveOne',[11111])ans=15.72251.371715.53591.276860.3579此时的所求值便为函数的待定系数。所以可得最小二乘曲线的表达式为：3579.60)2768.1sin(5359.15)3717.1cos(7225.15)(xxf然后进行相应的绘制图形便可以求出所要求出的结果。结果分析：通过最小二乘法多项式同样可以绘制出三角函数的曲线。并且程序通用性高，只要输入相应的插值便可以计算出结果，结果系数的小数点有效位同时也可以自己确定。3.P2291实验描述：由题意可知，由于这里涉及到了样条线的运算，计算较为复杂。且要涉及到画图的部分，所以此处采用matlab计算较为方便快捷。而书本上给出了一个用来计算三次紧压样条曲线的可调用函数，现在将其引用，并根据已知点得出相应的曲线。《数值方法》实验报告9实验结果：代入后得出的结果如下所示：X=[02468];Y=[040160300480];S=csfit(X,Y,0,98)S=0.81258.375000-2.437513.250043.250040.00001.4375-1.375067.0000160.0000-0.81257.250078.7500300.0000由结果可知此插值为在区间[0,8]中有三个分别划分了[0,2],[2,4],[4,6],[6,8]四个区间的插点。且多项式分别为230375.88125.0)(xxxS40)2(25.43)2(25.13)2(4375.2)(231xxxxS160)4(67)4(375.1)4(4375.1)(232