曲线拟合的数值计算方法实验教材.docx

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

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

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

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

5、lograrithmic function）的标准式形式为 b0时，Y随X增大而增大，先快后慢；b fminsearch(fiveOne,1 1 1 1 1)ans = 15.7225 1.3717 15.5359 1.2768 60.3579此时的所求值便为函数的待定系数。所以可得最小二乘曲线的表达式为：然后进行相应的绘制图形便可以求出所要求出的结果。结果分析：通过最小二乘法多项式同样可以绘制出三角函数的曲线。并且程序通用性高，只要输入相应的插值便可以计算出结果，结果系数的小数点有效位同时也可以自己确定。3.P229 1实验描述：由题意可知，由于这里涉及到了样条线的运算，计算较为复杂。且要涉

6、及到画图的部分，所以此处采用matlab计算较为方便快捷。而书本上给出了一个用来计算三次紧压样条曲线的可调用函数，现在将其引用，并根据已知点得出相应的曲线。实验结果：代入后得出的结果如下所示： X=0 2 4 6 8; Y=0 40 160 300 480; S=csfit(X,Y,0,98)S = 0.8125 8.3750 0 0 -2.4375 13.2500 43.2500 40.0000 1.4375 -1.3750 67.0000 160.0000 -0.8125 7.2500 78.7500 300.0000由结果可知此插值为在区间0,8中有三个分别划分了0,2,2,4,4,6,

7、6,8四个区间的插点。 且多项式分别为 在matlab中通过polyval作出相应的曲线，再用plot函数便可绘制图线。绘制后图线如下：此时便可以直观的看到一个平滑的样条线。结果分析：通过给定的一些数据点，便可以绘制出过这些点的相应的样条线，通过观察能发现样条线的平滑度与你选择的样条线类型以及数据点的分布有一定关联。不仅仅是紧压样条线，相关其它的线也可以用同样的方法一一给出。4.P238 5实验描述：由题意可知，题目是想求出所绘制数据点的三角多项式的逼近，三角多项式逼近的公式为： 此外，x为离散傅里叶级数，满足条件： 实验结果：而所给的点x为1,24的离散数值点，所以无法直接对其作出逼近公式，

8、需要进行尺度变换，将x点转换为: （1）通过matlab绘制出相关的三角多项式曲线，然后同样通过matlab的绘制点，将点绘制到这个曲线之中，具体的matlab代码如下： hold on; %保持图形不动，绘制新的图形入曲线中X=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的取值Y=58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 63 63 62 61 60 60 59 58; %数据点的y的取值plot(X,Y,xk) %绘制出数据点（2）然后便可以

9、画出如图所示的插值数据点。 结果分析：三角多项式的曲线拟合度非常高，能很好的绘制出图像的具体形式而且曲线平滑，但是它需要满足x属于-pi到pi的区间内。 5.P246 1实验描述：由题意可知，首先以3阶贝塞尔函数为带求函数。其求解格式如下：其中b为伯恩斯坦多项式，其定义如下： 用matlab先设置一个参数t,然后再根据公式，通过所给的数据点将伯恩斯坦多项式，以及x,y的关于参数t的多项式求出来。然后再把t设置为精度足够大的单位序列，绘制图线即可得出贝塞尔的效果。实验结果：其matlab运行后结果如下：【以第一组控制点为例】 X=0 1 1 3X = 0 1 1 3 Y=0 2 1 0Y = 0

10、 2 1 0 fiveTwo(X,Y)x = 3*t*(t - 1)2 - 3*t2*(t - 1) + 3*t3 y = 6*t*(t - 1)2 - 3*t2*(t - 1) 且绘制出第一组控制点所在位置以及三阶贝塞尔曲线如下所示：同理，第2,3组控制点所作的图形如下所示： 结果分析：可以通过Matlab函数绘制出三阶的贝塞尔函数。只要给出控制点，便可自动绘制出所控制的三阶贝塞尔函数以及控制点的位标。由此可以观察到贝塞尔与控制点的约束作用，以及所求得贝塞尔函数是个相对平滑的曲线。5、实验结论曲线拟合对于实际的工程以及理论推导问题都有着重大的作用。在具体的实验，数据分析上，往往我们只有巨量的

11、已知离散点，想从离散点中得出函数表达式则就需要曲线拟合进行，根据不同的要求，我们可以选择最小二乘法的幂函数拟合或者是一次函数，二次函数拟合，抑或是精度非常高的傅里叶变换的三角函数拟合。同时在建模方面，贝塞尔函数的引用也从数学层面解决了如何用计算机绘制出光滑圆润的曲线，在一些设计软件中，例如3dmax和maya的三维建模，Auto CAD的工程建模，贝塞尔运算对于曲线曲面的设计有着举足轻重的作用。附件（代码）:1.P202 1#include#include#includevoid main() extern float afunction(float b); extern float bfun

12、ction(float z,float v); float *a,y,q,c; int i; a=new float5; /定义动态数组a printf(请在此输入x的各插值n); scanf(%f %f %f %f %f,&a0,&a1,&a2,&a3,&a4); /手动输入x各值 q=afunction(a); /对x值进行求和 a=new float5; /重新为a进行数据分配 printf(请在此输入y的各插值n); scanf(%f %f %f %f %f,&a0,&a1,&a2,&a3,&a4); /手动输入y各值 c=afunction(a); /对y值进行求和 y=bfunct

13、ion(q,c); /计算出k的系数的大小 printf(此数据的最小二乘法曲线表达式为n Y=%f*xn,y); return;float afunction(float b) float x=0; int i; for(i=0;ifminsearch(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 63 63 62 61 60

14、60 59 58; %数值点的x取值 plot(x1,y1,x) %绘制出24个数值点的图形 x=linspace(0,25,100); hold on plot(x,y,k) %绘制出此函数的二乘法后的函数曲线3.P229 1 （三次紧压样条线计算函数）function S=csfit(x,y,dx0,dxn)N=length(x)-1;H=diff(x);D=diff(y)./H; %d(k)的值A=H(2:N-1);B=2*(H(1:N-1)+H(2:N);C=H(2:N);U=6*diff(D); %U(k)的值B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1)-dx

15、0);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(dxn-D(N); for k=2:N-1 temp=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);for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2; %三次紧压约束中S(0)的值M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2; %三次紧压约束中S(k)的值for k=0:N-1

16、 S(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1); %S(k,3)的值 S(k+1,2)=M(k+1)/2; %S(k,2)的值 S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2)/6; %S(k,1)的值 S(k+1,4)=y(k+1); %S(k,0)的值End（曲线绘制函数程序） 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);

17、 %第三段样条线 x4=6:.01:8; y4=polyval(S(4,:),x4-X(4); %第四段样条线 plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,x) %绘制连接成完整曲线，并且标注数据点的位置。 #include#include#includevoid main() extern 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

18、,15,16,17,18,19,20,21,22,23,24; float b24=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,

19、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 afunction(float X,float Y,int temp) float A=0; int i; for(i=0;i=24;i+) A=Yi*cos(temp*Xi)/6+A; return A;float bfunction(float X,float Y,int temp) float A=0; int i; for(i=0;i=24;i+) A=Yi