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

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

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

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

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

5、，其中横坐标是确定的。最小二乘幕函数拟合曲线 的系数A为：NNA ( xJyk)/(xj)k 1k 14.对数函数拟合：对数函数(lograrithmic function)的标准式形式为Y a bl nX(X 0)b0时，Y随X增大而增大，先快后慢；b fminsearch(fiveOne,1 1 1 1 1) ans =此时的所求值便为函数的待定系数。所以可得最小二乘曲线的表达式为：f (x) 15.7225cos(1.3717x) 15.5359sin(1.2768) 60.3579然后进行相应的绘制图形便可以求出所要求出的结果。结果分析：通过最小二乘法多项式同样可以绘制出三角函数的曲线

6、。并且程序通用性 高，只要输入相应的插值便可以计算出结果， 结果系数的小数点有效位同时也可 以自己确定。1实验描述：由题意可知， 由于这里涉及到了样条线的运算， 计算较为复杂。 且要涉及到 画图的部分，所以此处采用matlab计算较为方便快捷。而书本上给出了一个用来计算三次紧压样条曲线的可调用函数，现在将其引 用，并根据已知点得出相应的曲线。实验结果：代入后得出的结果如下所示: X=0 2 4 6 8; Y=0 40 160 300 480; S=csfit(X,Y,0,98)S =0 0由结果可知此插值为了在区间0,8中有三个分别划分了0,2,2,4,4,6,6,8四个区间的插点。且多项式分

7、别为So(x)30.8125X8.375 x2Sdx)2.4375(x32)13.25(x2)243.25(x2)40S2(x)1.4375(x4)31.375(x4)267(x 4)160Sa(x)0.8125(x6)37.25(x6)278.75(x6)300在matlab中通过polyval作出相应的曲线，再用plot函数便可绘制图线 绘制后图线如下：此时便可以直观的看到一个平滑的样条线。结果分析:通过给定的一些数据点，便可以绘制出过这些点的相应的样条线，通过观察能发现样条线的平滑度与你选择的样条线类型以及数据点的分布有一定关联。不仅仅是紧压样条线，相关其它的线也可以用同样的方法一一给出

8、。实验描述：由题意可知，题目是想求出所绘制数据点的三角多项式的逼近, 逼近的公式为：a0MTM(X)-(ajcos(jx) bjSin(jx)2j i此外，X为离散傅里叶级数，满足条件：实验结果：而所给的点X为1,24的离散数值点，所以无法直接对其作出逼近公式，需 要进行尺度变换，将X点转换为：2jXj, j0,1,.,24j24(1)通过matlab绘制出相关的三角多项式曲线，然后同样通过matlab的 绘制点，将点绘制到这个曲线之中，具体的matlab代码如下：hold on; %保持图形不动，绘制新的图形入曲线中X=1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

9、16 17 18 19 20 21 22 2324;%数据点的X的取值Y=58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 63 63 62 61 6060 59 58;%数据点的y的取值plot(X,Y,xk)%绘制出数据点(2)然后便可以画出如图所示的插值数据点三角多项式Xj2jN，j0,1,., N结果分析：三角多项式的曲线拟合度非常高，能很好的绘制出图像的具体形式而且曲线 平滑，但是它需要满足x属于-pi到pi的区间内。1实验描述：由题意可知，首先以33阶贝塞尔函数为带求函数。其求解格式如下：P(t)Bi,3(t)其中b为伯恩斯坦多项式，其

10、定义如下：3i 3 iBi,3(t)ti(1 t)3 ii用matlab先设置一个参数t,然后再根据公式，通过所给的数据点将伯恩斯 坦多项式，以及x,y的关于参数t的多项式求出来。然后再把t设置为精度足够 大的单位序列，绘制图线即可得出贝塞尔的效果。实验结果：其matlab运行后结果如下：【以第一组控制点为例】 X=0 1 1 3X =0113 Y=0 2 1 0Y =0 2 1 0 fiveTwo(X,Y)x =3*t*(t - 1)A2 - 3*tA2*(t - 1) + 3*tA36*t*(t - 1)A2 - 3*tA2*(t - 1)且绘制出第一组控制点所在位置以及三阶贝塞尔曲线如下

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

12、学层面解决了如何用计算机 绘制出光滑圆润的曲线，在一些设计软件中，例如3dmax和maya的三维建模，Auto CAD的工程建模，贝塞尔运算对于曲线曲面的设计有着举足轻重的作用。附件（代码）#include#include #include void main()extern float afunction(float b);extern float bfunction(float z,float v);float *a,y,q,c;int i;a=new float5;cos(i*B)+C.*si n( i*D)+E-y(i)A2+z;涵数的线性最小二乘法的多项式end之后在Matlab的命

13、令窗口输入如下语句：fminsearch(fiveOne,1 1 1 1 1)(绘制图形方程组) x1=0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24; %数值点的x取值 y1=58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 6363 62 61 60 60%数值点%绘制出24个数值点的图形%绘制出此函数的二乘法后的函数1(三次紧压样条线计算函数)function S=csfit(x,y,dx0,dxn) N=length(x)-1;H=diff(x);D=diff

14、(y)./H;59 58;的x取值 plot(x1,y1,x) x=linspace(0,25,100); hold on plot(x,y,k)曲线%d(k)的值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);for k=2:N-1%U(k)的值temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-

15、1);end M(N)=U(N-1)/B(N-1);for k=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;的值for k=0:N-1S(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1);S(k+1,2)=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(曲线绘制函数程序) x1=0:.01:2; y1=polyval(S(1,

16、:),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); plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,x)并且标注数据点的位置。%三次紧压约束中S(0)%三次紧压约束中S(k)%S(k,3)的值%S(k,2)的值%S(k, 1)的值%S(k,0)的值%第三段样条线%第四段样条线%绘制连接成完整曲线，#include#include#includevoid main()exte

17、rn float afunction(float X,float Y,int temp);extern float bfunction(float X,float Y,int temp);int i;float a24=1,2,3,4,5,6,7,8,9,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=new float7,d=new float7;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;float