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

1、曲线拟合的数值计算方法实验【摘要】实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合(curvefitting)是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。曲线直线化是曲线拟合的重要手段之一。对于某些非线性的资料可以通过简单的变量变换使之直线化，这样就可以按最小二乘法原理求出变换后变量的直线方程，在实际工作中常利用此直线方程绘制资料的标准工作曲线，同时根据需要可将此直线方程还原为曲线方程，实现对资料的曲线拟合。常用的曲线拟合有最小二乘法拟合、幕函数拟合、对数函数拟合、线性插值

2、、三次样条插值、端点约束。关键词曲线拟合、最小二乘法拟合、幕函数拟合、对数函数拟合、线性插值、三次样条插值、端点约束一、实验目的1 .掌握曲线拟合方式及其常用函数指数函数、幕函数、对数函数的拟合2 .掌握最小二乘法、线性插值、三次样条插值、端点约束等。3 .掌握实现曲线拟合的编程技巧。实验原理1 .曲线拟合曲线拟合是平面上离散点组所表示的坐标之间的函数关系的一种数据处理方法。用解析表达式逼近离散数据的一种方法。在科学实验或社会活动中，通过实验或观测得到量x与y的一组数据对(Xi,Yi)(i=1,2,.m),其中各Xi是彼此不同的。人们希望用一类与数据的背景材料规律相适应的解析表达式，y=f(x

3、,c)来反映量x与y之间的依赖关系，即在一定意义下最佳”地逼近或拟合已知数据。f(x,c)常称作拟合模型，式中C=(C1,C2,Cn)是一些待定参数。当C在f中线性出现时，称为线性模型，否则称为非线性模型。有许多衡量拟合优度的标准，最常用的一种做法是选择参数c使得拟合模型与实际观测值在各点的残差(或离差)ek=yf(fk,c)的加权平方和达到最小，此时所求曲线称作在加权最小二乘意义下对数据的拟合曲线。有许多求解拟合曲线的成功方法，对于线性模型一般通过建立和求解方程组来确定参数，从而求得拟合曲线。至于非线性模型，则要借助求解非线性方程组或用最优化方法求得所需参数才能得到拟合曲线，有时称之为非线性

4、最小二乘拟合。曲线拟合：贝塞尔曲线与路径转化时的误差。值越大，误差越大；值越小，越精确。2 .最小二乘法拟合:最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化嫡用最小二乘法来表达。函数曲线为：y=Ax+B其中系数满足下列的正规方程：NN、NN、NIZx2A+工xk旧=£xkykIkTJIkTJkTNN'xkANB八ykk耳k=13 .幕函数拟合:函数曲线为：设xk,y

5、k心有N个点，其中横坐标是确定的。最小二乘幕函数拟合曲线的系数A为：NNA=xMyk)/x；M)k1kT、4 .对数函数拟合:对数函数(lograrithmicfunction)的标准式形式为Y=ablnX(X0)b>0时，Y随X增大而增大，先快后慢；b<0时，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-

6、blnX5 .线性插值：在代数才S值中，为了提高插值多项式对函数的逼近程度一般是增加节点的个数，即提高多项式的次数，但这样做往往不能达到预想的效果。如下图所示：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现象。它表明通过增加节点来提高逼近程度是不适宜的，因而不采用高次多项式插值

7、。如果我们把以上的点用直线连接起来，显然比L7(x)要更逼近f(x)0这就是分段线性插值。而在实际应用中通常采用分段低次插值来提高近似程度，比如可用分段线性插值或分段三次埃尔米特插值来逼近已知函数，但它们的总体光滑性较差。为了克服这一缺点，一种全局化的分段插值方法一一三次样条插值成为比较理想的工具。数值方法实验报告6 .三次样条插值:设xk,ykW_q有N+1个点，其中a=x0<x1<x2<.<xN=b。如果存在N个三次多项式Sk(x)，系数为Sk0,Sk1,Sk2,Sk3满足如下性质:k'/k,OJk,1，k,2，k,3c/、C/、/、/、2/、3S(x)=S

8、k(x)=Sk,0,Sk,i(x-xk).Sk,2(x-xk)-Sk,3(x-xk)xxk,xki,k=0,1,.,N-1k=0,1,.Nk=0,1,.N-2k=0,1,.N-2k=0,1,.N-2S(xJ=ykS(xk1)-Sk1(xk1)S'(xk1)-S'k1(xk1)S''(xk1)=S''k1(xk1)则成函数S(x)为三次样条函数。7 .端点约束：紧压样条：存在唯一的三次样条曲线，其一阶导数的边界条件是：S'(a)=d0,S(b)=dNnatural样条：存在唯一的三次样条曲线，它的自由边界条件是：S''(a)

9、=0,S''(b)=0外推样条：存在唯一的三次样条曲线，其中通过对点x1和x2进行外推得到S“(a),同时通过对点X(n-1)和X(N-2)进行外推得到SKb)端点曲率调整：存在唯一的三次样条曲线，其中二阶导数的边界条件S”(a)和S”(b是确定的。抛物线终结样条：存在唯一的三次样条曲线，其中二阶在区间X0,X1内S"(x)三0,而在Xn-1,Xn内S”(x)三0。三、实验内容1.P2021X为拉伸胡克定律指出F=kx,其中F是拉伸弹簧的拉力（单位为盎司），的长度（单位为英寸）。根据下列试验数据，求解拉伸常量k的近似值xkFk0.23.610.47.310.610.9

10、0.814.511.018.2(b)XkFk0.215.30.410.60.615.90.8121.21.026.42.P2151洛杉矶（美国城市）郊区11月8日的温度记录入下表所示，其中共有24个数据点。（a）根据例5.5中的处理过程（使用fmins命令），对给定的数据求解最小二乘曲线f(x)=Acos(Bx)+Csin(Dx)+E。（b）求Ez（f）。（c）在同一坐标系下画出这些点集和（a）得出的最小二乘曲线。时间，p.m.温度时间，a.m.温度1P6615812662583653584r644581563557663657762757861858960960106010641159116

11、7午夜58正午68一个轿车在时间Tk时经过的距离dk,如下表所示。使用程序5.3,并根据阶导数边界条件S'(0)=0,S'(8)=98,求这些数据的三次紧压样条插值。时间，tk02468距离，dk0401603004804.P2385美国洛杉矶郊区11月9日的温度(华氏温度)如表5.10所示。采用24小时制。(a)求三角多项式T7(x)(b)在同一坐标系下，画出图T7(x)和24个数据点。(c)使用本地的温度情况重新求解问题(a)和问题(b)。时间,p.m温度时间,a.m温度1661582662583653584644585635576636577627578618589609

12、601060106411591167午夜58正午685.P2461编写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,所以只需用方程中NN(ZXk)A=£yk便可计算出k的数值。k工k4定义一个动态数组a【】,用来依次存取x和

13、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请按任意键继续。结果分析：易得，最小二乘法多项式计算可以很好的做出较为精

14、确的最小二乘法拟合曲线，并且程序通用性高，只要输入相应的插值便可以计算出结果，结果系数的小数点有效位同时也可以自己确定。实验描述：给出的最小二乘曲线表达式为：f(x)=Acos(Bx)Csin(Dx)E其中变量有5个，而给出的数据点有24个，在C语言中可以用牛顿-拉夫森算法迭代计算分别得出5个变量的值，但是方法繁琐，且迭代计算量庞大，因此这里采用Matlab进行相关的计算，调用fminsearch函数，求得当5个参量都为1附近时候多项式的最小值，此时便可求出此5个参变量的值.然后继续通过Matlab,将得到的公式以及各点，用plot函数，便可以求得。实验结果：运行matlab结果如下：>

15、>fminsearch('fiveOne',11111)ans=15.72251.371715.53591.276860.3579此时的所求值便为函数的待定系数。所以可得最小二乘曲线的表达式为：f(x)=15.7225cos(1.3717x)15.5359sin(1.2768)60.3579然后进行相应的绘制图形便可以求出所要求出的结果。结果分析：通过最小二乘法多项式同样可以绘制出三角函数的曲线。并且程序通用性高，只要输入相应的插值便可以计算出结果，结果系数的小数点有效位同时也可以自己确定。3.P2291实验描述：由题意可知，由于这里涉及到了样条线的运算，计算较为复杂。且

16、要涉及到画图的部分，所以此处采用matlab计算较为方便快捷。而书本上给出了一个用来计算三次紧压样条曲线的可调用函数，现在将其引用，并根据已知点得出相应的曲线。实验结果:代入后得出的结果如下所示:> >X=02468;> >Y=040160300480;> >S=csfit(X,Y,0,98)S=0043.250040.000067.0000160.000078.7500300.00000.81258.3750-2.437513.25001.4375-1.3750-0.81257.2500由结果可知此插值为在区间0,8中有三个分别划分了0,2,2,4,4,6

17、,6,8四个区间的插点。且多项式分别为S0(x)=0.812x58.37X232Si(x)=-2.4375(x-2)13.25(x-2)43.25(x-2)40S2(x)=1.4375(x-4)3-1.375(x-4)267(x-4)160S3(x)=-0.8125(x-6)37.25(x-6)278.75(x-6)300平滑的样条线在matlab中通过polyval作出相应的曲线，再用plot函数便可绘制图线。绘制后图线如下：此时便可以直观的看到一个结果分析:通过给定的一些数据点，便可以绘制出过这些点的相应的样条线，通过观察能发现样条线的平滑度与你选择的样条线类型以及数据点的分布有一定关联。

18、不仅仅是紧压样条线，相关其它的线也可以用同样的方法一一给出。4.P2385实验描述：由题意可知，题目是想求出所绘制数据点的三角多项式的逼近，三角多项式逼近的公式为：a0MTm(x)=-'(ajcos(jx)bjsin(jx)2jd止匕外，x为离散傅里叶级数，满足条件：2j二xj=-二j,j=0,1,.,NjN实验结果：而所给的点x为1,24的离散数值点，所以无法直接对其作出逼近公式，需要进行尺度变换，将x点转换为：2j二Xj一j,j=0,1,.,24j24(1)通过matlab绘制出相关的三角多项式曲线，然后同样通过matlab的绘制点，将点绘制到这个曲线之中，具体的matlab代码如

19、下：>>holdon;%保持图形不动，绘制新的图形入曲线中>>X=123456789101112131415161718192021222324;%数据点的x的取值>>Y=585858585757575860646768666665646363626160605958;%数据点的y的取值>>plot(X,Y,?kk?%绘制出数据点(2)然后便可以画出如图所示的插值数据点。结果分析:三角多项式的曲线拟合度非常高，能很好的绘制出图像的具体形式而且曲线平滑，但是它需要满足x属于-pi到pi的区间内。5.P2461实验描述：由题意可知，首先以回阶贝塞尔函

20、数为带求函数。其求解格式如下：P(t)="Bi,3(t)其中b为伯恩斯坦多项式，i4定义如下：Bi,3(t)=.ti(1-t)3"U)用matlab先设置一个参数t,然后再根据公式，通过所给的数据点将伯恩斯坦多项式，以及x,y的关于参数t的多项式求出来。然后再把t设置为精度足够大的单位序列，绘制图线即可得出贝塞尔的效果。实验结果：其matlab运行后结果如下：【以第一组控制点为例】>>X=0113X=0113>>Y=0210Y=0210>>fiveTwo(X,Y)x=3*t*(t-1)A2-3*tA2*(t-1)+3*tA36*t*(t-

21、1卜2-3中2*。-1)且绘制出第一组控制点所在位置以及三阶贝塞尔曲线如下所示:同理，第2,3组控制点所作的图形如下所示:结果分析：可以通过Matlab函数绘制出三阶的贝塞尔函数。只要给出控制点，便可自动绘制出所控制的三阶贝塞尔函数以及控制点的位标。由此可以观察到贝塞尔与控制点的约束作用，以及所求得贝塞尔函数是个相对平滑的曲线。五、实验结论曲线拟合对于实际的工程以及理论推导问题都有着重大的作用。在具体的实验，数据分析上，往往我们只有巨量的已知离散点，想从离散点中得出函数表达式则就需要曲线拟合进行，根据不同的要求，我们可以选择最小二乘法的幕函数拟合或者是一次函数，二次函数拟合，抑或是精度非常高的

22、傅里叶变换的三角函数拟合。同时在建模方面，贝塞尔函数的引用也从数学层面解决了如何用计算机绘制出光滑圆润的曲线，在一些设计软件中，例如3dmax和maya的三维建模，AutoCAD的工程建模，贝塞尔运算对于曲线曲面的设计有着举足轻重的作用。附件（代码）1.P2021#include<stdio.h>#include<math.h>#include<stdlib.h>voidmain()(externfloatafunction(floatb);externfloatbfunction(floatz,floatv);float*a,y,q,c;inti;a=ne

23、wfloat5;/定义动态数组aprintf("请在此输入x的各插值n");scanf("%f%f%f%f%f",&a0,&a1,&a2,&a3,&a4);手动输入x各值q=afunction(a);/对x值进行求a=newfloat5;/重新为a口进行数据分printf("请在此输入y的各插值n");scanf("%f%f%f%f%f",&a0,&a1,&a2,&a3,&a4);手动输入y各值c=afunction(a);/对y值进行

24、求和y=bfunction(q,c);/计算出k的系数的大小printf("此数据的最小二乘法曲线表达式为nY=%f*xn",y);return;floatafunction(floatb口)(floatx=0;inti;for(i=0;i<5;i+)x=bi+x;returnx;floatbfunction(floatz,floatv)(floatx;x=v/z;returnx;)2.P2151functionz=fiveOne(u,y)A=u(1);B=u(2);C=u(3);D=u(4);E=u(5);Z=0;fori=1:1:24z=(A.*cos(i*B)+

25、C.*sin(i*D)+E-y(i)A2+z;%函数的线性最小二乘法的多项式end之后在Matlab的命令窗口输入如下语句：>>fminsearch(,fiveOne?,11111)(绘制图形方程组)>>x1=0123456789101112131415161718192021222324;%数值点的x取值>>y1=585858585757575860646768666665646363626160605958;%数值点的x取值>>plot(x1,y1,'x')%绘制出24个数值点的图形>>x=linspace(0,2

26、5,100);>>holdon>>plot(x,y,'k')%绘制出此函数的二乘法后的函数曲线3.P2291(三次紧压样条线计算函数)functionS=csfit(x,y,dx0,dxn)N=length(x)-1;H=diff(x);%d(k)的值D=diff(y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N);C=H(2:N);U=6*diff(D);B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1)-dx0);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(dxn-D(N);f

27、ork=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);endM(N)=U(N-1)/B(N-1);fork=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;值%U(k)的值%三次紧压约束中S"(0)的值%三次紧压约束中S"(k)的fork=0:N-1S(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1);S(k+1,2)

28、=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2)/6;S(k+1,4)=y(k+1);End%S(k,3)的值%S(k,2)的值%S(k,1)的值%S(k,0)的值(曲线绘制函数程序)>>x1=0:.01:2;y1=polyval(S(1,:),x1-X(1);>>x2=2:.01:4;y2=polyval(S(2,:),x2-X(2);>>x3=4:.01:6;y3=polyval(S(3,:),x3-X(3);>>x4=6:.01:8;y4=polyval(S(4,:),x4-X(4);>

29、;>plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,'x')注数据点的位置。%第一段样条线%第二段样条线%第三段样条线%第四段样条线%绘制连接成完整曲线，并且标#include<stdio.h>#include<stdlib.h>#include<math.h>voidmain()externfloatafunction(floatX,floatY口,inttemp);externfloatbfunction(floatX口,floatY,inttemp);inti;floata24=1,2,3,4,5,6,7,8,9

30、,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24;floatb24=58,58,58,58,57,57,57,58,60,64,67,68,66,66,65,64,63,63,62,61,60,60,59,58);float*c;float*d;c=newfloat7,d=newfloat7;for(i=0;i<=24;i+)ci=afunction(a,b,i);di=bfunction(a,b,i);)printf("由所给的数据点可以得出n三角多项式T7(x)中a的系数分别为n");printf("%f%f%f%f%f%f%fn",c0,c1,c2,c3,c4,c5,c6);printf("由所给的数据点可以得出n三角多项式T7(x)中b的系数分别为n");printf("%f%f%f%f%f%f%fn",d0,d1,d2,d3,d4,d5,d6);system("pause");return;)floatafunction(flo