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

实验三 函数逼近与曲线拟合一、问题的提出：函数逼近是指“对函数类A 中给定的函数)(x f ,记作A x f ∈)(,要求在另一类简的便于计算的函数类B 中求函数A x p ∈)(,使 )(x p 与)(x f 的误差在某中度量意义下最小”。

函数类A 通常是区间],[b a 上的连续函数，记作],[b a C ,称为连续函数空间，而函数类B 通常为n 次多项式，有理函数或分段低次多项式等，函数逼近是数值分析的基础。

主要内容有：（1）最佳一致逼近多项式（2）最佳平方逼近多项式（3）曲线拟合的最小二乘法二、实验要求：1、构造正交多项式；2、构造最佳一致逼近；3、构造最佳平方逼近多项式；4、构造最小二乘法进行曲线拟合；5、求出近似解析表达式，打印出逼近曲线与拟合曲线，且打印出其在数据点上的偏差；6、探讨新的方法比较结果。

三、实验目的和意义:1、学习并掌握正交多项式的MATLAB 编程；2、学习并掌握最佳一致逼近的MATLAB 实验及精度比较；3、学习并掌握最佳平方逼近多项式的MATLAB 实验及精度比较；4、掌握曲线拟合的最小二乘法；5、最小二乘法也可用于求解超定线形代数方程组；6、 探索拟合函数的选择与拟合精度之间的关系；四、 算法步骤：1、正交多项式序列的生成{n ϕ（x ）}∞0：设n ϕ（x ）是],[b a 上首项系数a ≠n 0的n 次多项式，)(x ρ为],[b a 上权函数，如果多项式序列{n ϕ（x ）}∞0满足关系式⎩⎨⎧=>≠==⎰.,0,,0)()()()(),(k j A k j x d x x x kk j bak j ϕϕρϕϕ则称多项式序列{n ϕ（x ）}∞0为在],[b a 上带权)(x ρ正交，称n ϕ（x ）为],[b a 上带权)(x ρ 的n 次正交多项式。

1）输入函数)(x ρ和数据b a ,;2）分别求))(),(()),(,(x x x x j j j nϕϕϕ的内积； 3）按公式①)())(),(())(,()(,1)(10x x x x x x x x j n j j jj n nn ϕϕϕϕϕϕ∑-=-==计算)(x n ϕ，生成正交多项式；流程图：开始否是结束2、 最佳一致逼近多项式],[)(b a C x f ∈，若存在n n H x P ∈)(*使得n n E P f =∆),(*，则称)(*x P n 是)(x f 在],[b a 上的最佳一致逼近多项式或最小偏差逼近多项式，简称最佳逼近多项式。

现在我们所求的是最佳一次逼近多项式x x P 101)(αα+=,其中)()()(21x f ab b f a f '=--=α ①2)()(2)()(220x a a b a f b f x f a f +--+=α ② 1）输入函数)(x f 和数据b a ,； 2）计算1α和)(x f '； 3）求2x 和)(2x f ； 4）按公式②，计算0α； 5）生成最佳一次逼近多项式；流程图：3、最佳平方逼近多项式对],[)(b a C x f ∈及],[b a C 中的一个子集)}(),...,(),({10x x x span n ϕϕϕϕ=，若存在.)(*ϕ∈x S 使||.)]()()[(min||)()(||min ||)()(2)(22)(22*dx x S x f x x S x f x S x f bax S x S ⎰-=-=-∈∈ρϕϕ则称)(*x S 是)(x f 在子集],[b a C ⊂ϕ中的最佳平方逼近函数。

若取],1,0[)(,1)(,)(C x f x x x kk ∈≡=ρϕ则要在n H 中求n 次最佳平方逼近多项式,...)(**1*0*nn x a x a a x s +++=此时,11))(),((1++==⎰+j k dx x x x j k k j ϕϕk k k d dx x x f x x f ≡=⎰10)())(),((ϕ若用H 表示),...,,1(nn x x G G =对应的矩阵，既⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡+++++=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)12/(1)2/(1)1/(1)2/(13/12/1)1/(12/11),(),(),(),(),(),(),(),(),(11110101000n n n n n H n n n n n n nϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ称为希尔伯特距阵，记,),...,,(,),...,,(1010Tn T n d d d d a a a a ==则d Ha =的解),...,1,0(*n k a a kk ==即为所求。

平方误差为))(),((||)(||||)(||0*2222x f x a x f x k nk k ϕδ∑=-=；1）输入函数)(x f 和数据b a ,；2）求Tn d d d d ),...,,(10=；3）解方程组 d Ha =，解出 Tn a a a a ),...,,(10=；4）生成最佳平方逼近多项式；流程图：4、曲线拟合的最小二乘法由已知的离散数据点选择与实验点误差最小的曲线)(...)()()(1100x a x a x a x S n n ϕϕϕ+++= 称为曲线拟合的最小二乘法。

若记 ),()()(),(0i k i jmi ik j x x x ϕϕωϕϕ∑==k i k i mi i k d x x f x f ≡=∑=)()()(),(0ϕωϕ上式可改写为),...,1,0(;),(n k d a k j noj j k-=∑=ϕϕ这个方程成为法方程，可写成距阵形式.d Ga =其中,),...,,(,),...,,(1010Tn T n d d d d a a a a ==⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=),(),(),()(),(),(),(),(),(101110101000n n n n n n G ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ 。

它的平方误差为：.)]()([)(||||222iimi ix f x S x -=∑=ωδ1） 按照最小二乘法的性质构造Gram 矩阵G ，并求解Ga=d ；构造的时候首先构造一个零矩阵A ；2）然后开始构造Gram 矩阵（在下面程序里我们把克莱姆矩阵用A 来表示）3）然后求列矩阵b ，因为Aa=b ，所以求 a=A\b ；（d 就是列矩阵b ）; 4）然后找对应数据的最小二乘拟合方程和画出它的图像；5）在m 文件里制好以上规定的程序后，在matlab 的命令窗口 输入数组x 和数组y 及所选择的拟合多项式次数 m ，然后执行就可以得到曲线二乘拟合方程和它的图像。

流程图：五、Matlab 程序源代码：1、正交多项式序列的生成function gouzaozhengjiaoduoxiangshi 数值例题：1）当区间为]1,1[-，权函数1)(≡x ρ时，由,...},...,,1{nx x 正交化得到的多项式就成为勒让德多项式，它的表达式为1)(0=x P ， })1{(!21)(2nnn nn x dxd n x P -=； 递推公式为,...),2,1(),()()12()()1(11=-+=+-+n x nP x xP n x P n n n n 由x x P x P ==)(,1)(10，利用（1）就可推出,8/)33035()(,2/)35()(,2/)13()(2443322+-=-=-=x x x p x x x P x x P…y 12)当权函数21)(x x -=ρ，区间为]1,1[-时，由序列,...},...,,1{nx x 正交化得到的多项式就是切比雪夫多项式，它可表示为.1||),arccos cos()(≤=x x n x T n 若令θcos =x ，则.0,cos )(πθθ≤≤=n x T n递推关系.)(,1)(,...),2,1(,)(2)(1011x x T x T n T x xT x T n n n ===-=-+ 可推出,188)(,34)(),12(2443322+-=-=-=x x x T x x x T x x T2、 最佳一致逼近多项式function yicibijin数值例题：（1）求21)(x x f +=在[0，1]上的最佳一次逼近多项式 ；结果为：x x P 414.0955.0)(1+=；误差限为.045.0|)(1|max 1210≤-+≤≤x P x x3、最佳平方逼近多项式function pingfang● 数值例题：设21)(x x f +=，求]1,0[上的一次最佳平方逼近多项式；结果为：x x S 426.0934.0)(*1+=，0026.0||)(||22=x δ；4、曲线拟合的最小二乘法：function p=nihe(x,y,m)● 数值实例：例1:下面给定的是乌鲁木齐最近1个月早晨7:00左右(新疆时间)的天气预报所得到的温度数据表 ,按照数据找出任意次曲线拟合方程和它的图像.(2006年10月26~11月26)解: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 1]m=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 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文件的时候我们把想得到的多项式的次数直接定义为(m=length(x))跟X的数量一样取得最小二乘法的曲线拟合方程,这样上术问题中的m=30,这时候它的程序代码为:f unction p=nn(x,y)赋值x,y :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 1]；运行p=nn(x,y)解为:p =图像为:结果分析:以上结果可以看到用最小二乘拟合来求解问题时,有时候他的结果很接近实际情况,有时候跟实际情况里的太远,因为所求得多项式次数太小时数据点之间差别很大.次数最大是误差最小但是有时后不符合实际情况，所以用最小二乘法时次数要取合适一点.或采用：function p=duoxiangxi(x,y,m)数值例题：（1）合成纤维抽丝工段中第一导丝盘的速度对丝的质量有很大的影响，第一丝盘的速度和电流周波有重要关系。