曲线拟合的最小二乘法讲解.docx

1、曲线拟合的最小二乘法讲解实验三函数逼近与曲线拟合、问题的提出：函数逼近是指“对函 数类A中给定的函数f(x),记作f(x)A,要求在另一类简的便于计算的函数类 B中求函数p(x)A,使p(x)与f (x)的误差在某中度量意义下最小”函数类A通常是区间a,b上的连续函数，记作 Ca,b,称为连续函数空间，而函数类B通常为n次多项式，有理函数或分段低次多项式等，函数逼近是数值分析的基础。主要内容有：(1)最佳一致逼近多项式(2)最佳平方逼近多项式 (3 )曲线拟合的最小二乘法1、 构造正交多项式；2、 构造最佳一致逼近；3、 构造最佳平方逼近多项式；4、 构造最小二乘法进行曲线拟合；5、 求出近似

2、解析表达式，打印出逼近曲线与拟合曲线，且打印出其在数据点上的偏差;6、 探讨新的方法比较结果。三、实验目的和意义:1、学习并掌握正交多项式的 MATLAB编程；4、 掌握曲线拟合的最小二乘法；5、 最小二乘法也可用于求解超定线形代数方程组；6、 探索拟合函数的选择与拟合精度之间的关系；四、算法步骤：an=0的n次多项式，r(x)为a,b上权函1、正交多项式序列的生成 ( X)o :设 ( X)是a,b上首项系数数，如果多项式序列 ( X)o：满足关系式(x)正交，称;:n (x)为a,b上带权(x)则称多项式序列 (X)o：为在a,b上带权的n次正交多项式。1 )输入函数 (x)和数据a,b;

3、2)分别求(xn, j(x),C j (x), j(x)的内积；. . n 2 (Xn,j(X), ,3)按公式；：o(X)=1, -(X) =Xn j j(X)计算；：n(X)，生成正交多项式;j鼻Wj(x),Wj(x)流程图：开始cz结束2、最佳一致逼近多项式f(x) Ca,b，若存在 R*(x) Hn 使得.：(f,P；En，则称 P； (x)是 f (x)在a,b上的最佳一致逼近多项式或最小偏差逼近多项式，简称最佳逼近多项式。现在我们所求的是最佳一次逼近多项式 R(x) = rx,其中f (a) f(X2) f (b) f (a) a X22 b a 21 )输入函数f (x)和数据a

4、, b ；2)计算：i和f (x)；3)求 X2和 f(X2)；4)按公式，计算I。；5)生成最佳一次逼近多项式;流程图:3、最佳平方逼近多项式对 f (x) Ca,b及Ca,b中的一个子集：二 span o(x), i(x)，n(x)，若存在* * 2 2 b 2S (x) 使II f (x) S (x) ll2 = m)nJI f(X) S(x) |2 = sm)%Ja P(x)f (x) - S(x) dx. 则称S*(x)是f(x)在子集 Ca,b中的最佳平方逼近函数。 k若取L(x)二x , T(x)三1, f (x) C0,1,则要在H n中求n次最佳平方逼近多项式* * * *

5、ns (x) =ao ax . a“x ,1 k * 1 1 k此时()(x)ar 0x jdrrn,(f(x), k(x)H。f(x)xd-dk右用H表示Gn=G (1, x,., x）对应的矩阵，既-他角）何0n)1-11/2 1/( n + 1)伴,0）（乳角）件Wn)1/21/3 1/(n +2)H =33 +39a9-件禹）（咒理1）伴 nWnLi i1/( n +1)1/(n +2) 1/(2 n+1)称为希尔伯特距阵，记a = (a0 , a1,.,an)T,d= (d,d1,.,dn)T,则Ha = d的解nak =a；(k =0,1, . .n)即为所求。平方误差为 |6(x

6、) |2斗| f(x)|；它 a：(x), f(x)；k01 )输入函数f (x)和数据a, b ；2)求 d =(d,d1,.,dn)T ；3) 解方程组 Ha=d，解出 a = (a0, a1,.,an)T ；4）生成最佳平方逼近多项式;流程图:4、曲线拟合的最小二乘法由已知的离散数据点选择与实验点误差最小的曲线S(x)二aoS(x) ai i(x) . - a(x)称为曲线拟合的最小二乘法。m m若记 (：j,k) f 异(X)：j(Xi)匚(Xi), (f, ：k) f：(Xi)f(X)：k(Xi)二 dk i=0 i=0n上式可改写为 v (Uj)aj二dk；(k-0,1,，n)这个

7、方程成为法方程，可写成距阵形式 j =oGa 二 d.其中 a =(ao,ai,.,an)T,d =(do,di,.,dn)T,Co0）（1）Co，n）|件,）伴1,1）件件）3 3 3fn,0）件円1）（咋件）一m它的平方误差为:2 2|.|2 八，(Xi)S(Xi) f(Xi).i兰零矩阵A ；|B|=2）然后开始构造 Gram矩阵 伽皿）（仍皿）j） = 丫叭坞）例佃）佻（珂）7 （在下面程序里我们把克莱姆矩阵用 A来表示）3）然后求列矩阵 b，因为Aa=b，所以求 a=Ab ； （ d就是列矩阵b） 4）然后找对应数据的最小二乘拟合方程和画出它的图像;5）在m文件里制好以上规定的程序后

8、，在 matlab的命令窗口 输入数组x和数组y及所选择的拟合多项式次数 m,然后执行就可以得到曲线二乘拟合方程和它的图像。流程图:开始 输入数组*尤和旳橫数m+L uuCJfAdh*EM X a Adir制牆零矩阵A*. * 执行循环语句far 1=0 :m, for jpOiin专 nsma求得直矩阵:克麴矩阵a+屮扈吨伝.仪少： i * _ 翹原聽 b(j+l)=sumtei *y) 求 a=Ab;*J求尸輕俎)； (l):0.01：x(end)- 越=pgjy嗖血用】)；，所找到的方程多项式)变换尸paly2tr氯捏);心用pht画出所找到的罢项武图像边慝a 1, a2 八 b ；岳材

9、.匕naarkErth:匕 2 D”五、Matlab程序源代码：1、正交多项式序列的生成function gouzaozhe ngjiaoduoxia ngshi数值例题：1 )当区间为-1,1，权函数r(x)=1时，由1, x,.,xn,.正交化得到的多项式就成为勒让1 d n德多项式，它的表达式为 P0(x)=1, Pn(x) n-比(x2-1)n；2 n! dx递推公式为(n 1)Pn1(x)=(2n 1)xPn(x)-nPn4(x), (n =1,2,.),由F (x) =1,P(x)二x，利用(1 )就可推出2P2(x)=(3x -1)/2,P3(x) =(5x3 -3x)/2,4

10、2p4(x)=(35x -30x 3)/8,由序列1,x,，xn,.正交化得到的多项式就是切比雪夫多项式，它可表示为 Tn（x） =cos（n arccosx）, | x |空1.若令 x = cost，贝U Tn(x)二 cosnv,0“：二.可推出递推关系 IxTxTndTnW，2，.），To（X）1 ,Ti（X）= X.2T2(x=2x -1),= 4x3 -3x,T3(x)2、最佳一致逼近多项式function yicibijin_ -数值例题：（1 ）求f（x） =-1 x在0, 1上的最佳一次逼近多项式结果为:R(x) =0.955 0.414x ；误差限为max | . 1 x?

11、 -R(x) |乞 0.045. 0知3、 最佳平方逼近多项式fun ctio n pingfang数值例题：设f(x)二.1 x2，求0,1上的一次最佳平方逼近多项式；结果为：S(X)=0.934 +0.426X , | 6(x)血=0.0026 ；4、 曲线拟合的最小二乘法：fun cti on p=n ihe(x,y,m)数值实例：例1:下面给定的是乌鲁木齐最近 1个月早晨7:00左右(新疆时间)的天气预报所得到的温度数据表，按照数据找出任意次曲线拟合方程和它的图像 （2006 年 10 月 2611 月 26）天数12345678910温度910111213141312119天数111

12、21314151617181920温度101112131412111098天数21222324252627282930温度78919765311解：x=1:30y=9 10 11 12 13 14 13 12 11 9 10 11 12 13 14 12 11 10 9 8 7 8 9 11 8 7 6 5 3 1m=9； ( m 是任意真整数 , 但是不要取的太大)运行 p=nihe(x,y,m)结果为 :p =图形为 :或者任意取部分数据也可得到相应得多项式和它的图像m=15 时x=1:30y=9 10 11 12 13 14 13 12 11 9 10 11 12 13 14 12 11

13、 10 9 8 7 8 9 11 8 7 6 5 3 1 运行 p=nihe(x,y,m)结果为：p =图形为：如果:m=3x=1:30y=9 10 11 12 13 14 13 12 11 9 10 11 12 13 14 12 11 10 9 8 7 8 9 11 8 7 6 5 3 1 运行 p=nihe(x,y,m)结果为：p =图形为：(m=length(x) 跟 X 的数如果写 M 文件的时候我们把想得到的多项式的次数直接定义为量一样取得最小二乘法的曲线拟合方程 ，这样上术问题中的 m=30,这时候它的程序代码为：function p=nn( x,y)赋值x,y :x=1:30y=

14、9 10 11 12 13 14 13 12 11 9 10 11 12 13 14 12 11 10 9 8 7 8 9 11 8 7 6 5 3 1 ； 运行p=nn (x,y)解为：P =图像为：结果分析：以上结果可以看到用最小二乘拟合来求解问题时 ，有时候他的结果很接近实际情况 ，有时候跟实际情况里的太远，因为所求得多项式次数太小时数据点之间差别很大 次数最大是误差最小但是有时后不符合实际情况，所以用最小二乘法时次数要取合适一点或采用：fun ctio n p=duoxia ngxi(x,y,m)数值例题：(1)合成纤维抽丝工段中第一导丝盘的速度对丝的质量有很大的影响，第一丝盘的速度和

15、电流周波有重要关系。下面是一组实例数据：X49.250.049.349.049.049.549.849.950.250.2y16.717.016.816.616.716.816.917.017.017.1其中x代表电流周波，y代表第一导丝盘的速度 冃的：研究第一导丝盘速度 y与电流周波x的关系。解：根据所给数据，在 matlab窗口运行:x=49.2 50.0 49.3 49.0 49.0 49.5 49.8 49.9 50.2 50.2y=16.7 17.0 16.8 16.6 16.7 16.8 16.9 17.0 17.0 17.1plot（x,y,p,markersize,15）17.

16、117.05171B.95IS.S1B.7516.6549.450SQ.290.4S1 （x）= a从图中看到各点在一条直线附近，故可选择线性函数作拟合曲线，既令故，根据曲线拟合的最小二乘法，公式（ 5.5）， （ 5.6），即10（:j, - k）八jk ；i =Q（K=09 1，ooooo）10dk=（ f, A）八讦k；（K=0, 1 , 0 0 0 o）i 这里，m=10,n=1,二 1， l（x）二 1, I （x）二 X,扌巴 x=49.2 50.0 49.3 49.0 49.0 49.5 49.8 49.9 50.2 50.2；y=16.7 17.0 16.8 16.6 16.7

17、 16.8 16.9 17.0 17.0 17.1；代入公式（5.5），（ 5.6），（%），（申0,叫）=（久），严汽），严0, f），（申1, f ）分别可得:10（；：0, ；：0） 10 ；i =110 10（：0, 1）=（ 1, ：0）八 =496.1， （ ：1, 1）Xi2 = 24614i =1 i =1(d fj 二 fi =168.6 , (仆 fj 二、Xi fi =8364.9i 4 i 4得方程组10 ao 496.1ai =168.6496.1a0 24614a1 =8364.9解得：a1 =0.3389 , a0= 0.0490。于是所求拟合曲线为：S = 0.

18、33886 x + 0.04896910 2误差平方和：| ： (x) 2=二 i s(Xi) - f (Xi) =0.0469i 二MATLAB上运行的结果：x=49.2 50.0 49.3 49.0 49.0 49.5 49.8 49.9 50.2 50.2;y=16.7 17.0 16.8 16.6 16.7 16.8 16.9 17.0 17.0 17.1;m=1;S =0.33886 x + 0.04896P =0.3389 0.0490程序运行的结果中 S( x)=0.33886 x + 0.048969为所求拟合曲线的方程下面看离散点x，f=y和所得的拟合曲线的逼近情况：在 ma

19、tlab窗口输入x=49.2 50.0 49.3 49.0 49.0 49.5 49.8 49.9 50.2 50.2;f=16.7 17.0 16.8 16.6 16.7 16.8 16.9 17.0 17.0 17.1;x仁49:0.01:50.4;s = 0.3389*x+0.0490;s=in li ne(s);s1=s(x1);图二误差平方和：在matlab窗口运行以下程序f=16.7 17.0 16.8 16.6 16.7 16.8 16.9 17.0 17.0 17.1; t=0;x=49.2 50.0 49.3 49.0 49.0 49.5 49.8 49.9 50.2 50.2;S=0.3389*x+0.0490;n=len gth(x);for i=1: nt=t+(f(i)-S(i)F2;en d,t结果：t=0.0469(2)已知一组实验数据如下，求它的拟合曲线。Xi12345fi44.5688.521311结果为：S；(x) =2.77 1.13x ；如图所示：e- + 1T.-SF&-$-4 5丹11 U 1111-1 n2 51AJA