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

《曲线拟合的数值计算方法试验》由会员分享，可在线阅读，更多相关《曲线拟合的数值计算方法试验（19页珍藏版）》请在装配图网上搜索。

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使得拟合模型与实际观测值在各点的残差(或离差)ek = yk-f (fk, c)的加权平方和到达最小，此时所求曲线称作在加权最小二乘意义下对数据的拟合曲线。 有许多求解拟合曲线的成功方法，对丁线 性模型一般通过建立和求解方程组来确定参数，从而求得拟合曲线。至丁非线性 模型，那么要借助求解非线

4、性方程组或用最优化方法求得所需参数才能得到拟合曲 线，有时称之为非线性最小二乘拟合。曲线拟合：贝塞尔曲线与路径转化时的误差。值越大，误差越大；值越小， 越精确。2. 最小二乘法拟合：最小二乘法(乂称最小平方法)是一种数学优化技术。它通过最小化误差的 平方和寻找数据的最正确函数匹配。利用最小二乘法可以简便地求得未知的数据， 并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用丁曲线拟合。其他一些优化问题也可通过最小化能量或最大化嫡用最小二乘法来 表达。函数曲线为：y=Ax+B其中系数满足以下的正规方程：,n *n * nIZ x2；A+ xk；B = Z xkykNN、xk A

5、 NB = yk3. 籍函数拟合：函数曲线为：设xk yk此有N个点，其中横坐标是确定的。最小二乘籍函数拟合曲线 的系数A为：NNA = C xk yk)/C x?M)k =1kT、4. 对数函数拟合：对数函数(lograrithmic function)的标准式形式为Y =a bln X(X 0)b0时，Y随X增X而增大，先快后慢；b0时，Y随X增X而减少，先快 后慢，见图12.4(c)、(d)。当以Y和lnX绘制的散点图呈直线趋势时，可考虑采 用对数函数描述Y与X之间的非线性关系，式中的b和a分别为斜率和截距。更一般的对数函数Y=a+bln(X+k)式中k为一常量，往往未知。(a) lnY

6、=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)也只能在很小的范

7、围内收敛， 这一现象称为Runge现象。它说明通过增加节点来提高逼近程度是不适宜的， 因而不采用高次多项式插值。如果我们把以上的点用直线连接起来，显然比L7(x)要更逼近f(x)。这就是分段线性插值。而在实际应用中通常采用分段低次插值来提高近似程度，比方可用分段线性插值或分段三次埃尔米特插值来逼近函数，但它们的总体光滑性较差。为了克服这一缺点，一种全局化的分段插值方法一一三次样条插值成为比拟理想的工 具。6. 三次样条插值:设 奴，yk汩工有N+1个点，其中a = x0入X2. fminsearch(fiveOne,1 1 1 1 1)ans =15.72251.371715.53591.27

8、6860.3579此时的所求值便为函数的待定系数。所以可得最小二乘曲线的表达式为：f (x) = 15.7225cos(1.3717x) 15.5359sin(1.2768) 60.3579然后进行相应的绘制图形便可以求出所要求出的结果。结果分析：通过最小二乘法多项式同样可以绘制出三角函数的曲线。并且程序通用性高，只要输入相应的插值便可以计算出结果， 结果系数的小数点有效位同时也可 以自己确定。3. P229 1实验描述：由题意可知，由丁这里涉及到了样条线的运算， 计算较为复杂。且要涉及到 画图的局部，所以此处采用 matlab计算较为方便快捷。而书本上给出了一个用来计算三次紧压样条曲线的可调

9、用函数，现在将其引用，并根据点得出相应的曲线。实验结果:代入后得出的结果如下所示: X=0 2 4 6 8; Y=0 40 160 300 480; S=csfit(X,Y,0,98)0.81258.3750-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四 个区间的插点。且多项式分别为S0(x) =0.8 1 2x5 8.3 7必S(x) =-2.4375(x-2)3 13.25(x-2)2 4

10、3.25(x-2) 40S2(x) =1.4375(x-4)3 -1.375(x-4)2 67(x-4) 1 6 0&(x) = -0.8125(x -6)3 7.25(x -6)2 78.75(x-6) 3 0 0在matlab中通过polyval作出相应的曲线，再用plot函数便可绘制图线。绘 制后图线如下：此时便可以直观的看到一个平滑的样条线结果分析:通过给定的一些数据点，便可以绘制出过这些点的相应的样条线，通过观察 能发现样条线的平滑度与你选择的样条线类型以及数据点的分布有一定关联。 不 仅仅是紧压样条线，相关其它的线也可以用同样的方法一一给出。4. P238 5实验描述：由题意可知，

11、题目是想求出所绘制数据点的三角多项式的逼近，三角多项式逼近的公式为：a0MTm (x) = (aj cos(jx) bj sin(jx)2 j此外，x为离散傅里叶级数，满足条件：2 j二xj = j, j = 0,1,.,N jN实验结果：而所给的点x为1,24的离散数值点，所以无法直接对其作出逼近公式，需 要进行尺度变换，将x点转换为：2j 二Xj = j, j = 0,1,.,24j 24(1) 通过matlab绘制出相关的三角多项式曲线，然后同样通过matlab的绘制点，将点绘制到这个曲线之中，具体的matlab代码如下：hold on;%保持图形不动，绘制新的图形入曲线中X=1 2 3

12、 4 5 6 7 8 9 10 11 12 13 14 15 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 60 60 5958;%数据点的y的取值plot(X,Y,?kk?%绘制出数据点(2) 然后便可以画出如下图的插值数据点。结果分析:三角多项式的曲线拟合度非常高，能很好的绘制出图像的具体形式而且曲线 平滑，但是它需要满足x届丁-pi到pi的区间内。5. P246 1实验描述：由题意可知，首先以a阶贝塞尔函数为带求函数。其求解格式如下：其中b为

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

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

15、也从数学层面解决了如何用计算机 绘制出光滑圆润的曲线，在一些设计软件中，例如 3dmax和maya的三维建模， Auto CAD的工程建模，贝塞尔运算对丁曲线曲面的设计有着举足轻重的作用。附件代码：1.P202 1#include#include#includevoid main()extern float afunction(float b);extern float bfunction(float z,float v);float *a,y,q,c;int i;a=new float5;定义动态数组aprintf(请在此输入x的各插值n);scanf(%f %f %f %f %f,&a0,

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

17、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 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=

18、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);%U(k)的值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-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)/

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

20、;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); 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()(ext

21、ern 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;

22、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 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;ix=X转换成元素的点乘除形式 y=Y转换成元素的点乘除形式 t=0:0.01:1;hold onplotx,y线matlab命令面板中输入如下指令%保持在一张图内绘制%在图形中绘制出所求得贝塞尔函数曲