曲线拟合的最小二乘法讲解.docx
《曲线拟合的最小二乘法讲解.docx》由会员分享，可在线阅读，更多相关《曲线拟合的最小二乘法讲解.docx（17页珍藏版）》请在冰豆网上搜索。
曲线拟合的最小二乘法讲解
实验三函数逼近与曲线拟合
、问题的提出：
函数逼近是指“对函数类A中给定的函数f（x）,记作f（x）・A,要求在另一类简
的便于计算的函数类B中求函数p（x）・A,使p（x）与f（x）的误差在某中度量意义
下最小”函数类A通常是区间[a,b]上的连续函数，记作C[a,b],称为连续函数空间，
而函数类B通常为n次多项式，有理函数或分段低次多项式等，函数逼近是数值分析
的基础。
主要内容有：
（1）最佳一致逼近多项式
（2）最佳平方逼近多项式
（3）曲线拟合的最小二乘法
1、构造正交多项式；
2、构造最佳一致逼近；
3、构造最佳平方逼近多项式；
4、构造最小二乘法进行曲线拟合；
5、求出近似解析表达式，打印出逼近曲线与拟合曲线，且打印出其在数据点上的偏差;
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;
2）分别求（xn,j（x））,Cj（x）,j（x））的内积；
..n2（Xn,®j（X））,,
3）按公式①；：
o（X）=1,-（X）=Xnjj（X）计算；：
n（X），生成正交多项式;
j鼻Wj（x）,Wj（x））
流程图：
开始
cz>
结束
2、最佳一致逼近多项式
f（x）C[a,b]，若存在R*（x）Hn使得.：
（f,P；^En，则称P；（x）是f（x）在[a,b]
上的最佳一致逼近多项式或最小偏差逼近多项式，简称最佳逼近多项式。
现在我们所求的是最佳一次逼近多项式R（x）=「rx,其中
f（a）f（X2）f（b）f（a）aX2
2b—a2
1）输入函数f（x）和数据a,b；
2）计算：
i和f（x）；
3）求X2和f（X2）；
4）按公式②，计算I。
；
5）生成最佳一次逼近多项式;
流程图:
3、最佳平方逼近多项式
对f（x）C[a,b]及C[a,b]中的一个子集：
二span{o（x）,i（x），…，’n（x）}，若存在
**22b2
S（x）"•使IIf（x）—S（x）ll2=m）nJIf（X）—S（x）||2=sm）%JaP（x）[f（x）-S（x）]dx.则称S*（x）是f（x）在子集〉C[a,b]中的最佳平方逼近函数。
k
若取L（x）二x,T（x）三1,f（x）C[0,1],则要在Hn中求n次最佳平方逼近多项式
****n
s（x）=aoa〔x...a“x,
1k*11k
此时（）（x）「ar0xjd^rrn,（f（x）,k（x）H。
f（x）xd-dk
右用H表示Gn
=G（1,x,...,x
）对应的矩阵，
既
-
他角）…
何0®n）1
-
-1
1/2…
1/（n+1）]
伴,®0）
（乳角）…
件Wn）
1/2
1/3…
1/（n+2）
H=
3
3+
3
—
9
a
9
-
件禹）
（咒理1）…
伴nWnL
ii
1/（n+1）
1/（n+2）…
1/（2n+1）
称为希尔伯特距阵，记
a=（a0,a1,...
an）T,d
=（d°,d1,...,dn）T,则
Ha=d的解
n
ak=a；（k=0,1,..n）即为所求。
平方误差为||6（x）||2斗|f（x）||；它a：
®（x）,f（x））；
k^0
1）输入函数f（x）和数据a,b；
2）求d=（d°,d1,...,dn）T；
3）解方程组Ha=d，解出a=（a0,a1,...,an）T；
4）生成最佳平方逼近多项式;
流程图:
4、曲线拟合的最小二乘法
由已知的离散数据点选择与实验点误差最小的曲线
S（x）二aoS（x）•aii（x）...-a^\（x）称为曲线拟合的最小二乘法。
mm
若记（：
：
j,「k）f异（X）：
j（Xi）匚（Xi）,（f,：
：
k）f』：
（Xi）f（X）：
k（Xi）二dki=0i=0
n
上式可改写为v（Uj）aj二dk；（k-0,1,…，n）这个方程成为法方程，可写成距阵形式j=o
Ga二d.
其中a=（ao,ai,...,an）T,d=（do,di,...,dn）T,
Co」0）（—1）…Co，n）
|件,％）伴1,®1）…件件）
333
fn,®0）件円1）…（咋件）一
m
它的平方误差为:
22
||「.||2八，（Xi）[S（Xi）—f（Xi）].
i兰
零矩阵A；
|B|=
2）然后开始构造Gram矩阵伽皿）（仍皿）
@j®）=丫叭坞）例佃）佻（珂）
7（在下面程序里我们把克莱姆矩阵用A来表示）
3）然后求列矩阵b，因为Aa=b，所以求a=A\b；（d就是列矩阵b）4）然后找对应数据的最小二乘拟合方程和画出它的图像;
5）在m文件里制好以上规定的程序后，在matlab的命令窗口输入数组x和数组y及所
选择的拟合多项式次数m,然后执行就可以得到曲线二乘拟合方程和它的图像。
流程图:
开始输入数组*尤和旳橫数m+
LuuCJfAdh«*^EMXaAdir
制牆零矩阵A*
.*
执行循环语句far1=0:
m,forjpOiin^专n
sm«a求得直矩阵:
克麴矩阵〉a°+屮扈吨伝.仪少：
i*_「翹原聽>b（j+l）=sumte^i*y）^
求a=A\b;*J
求尸輕俎）；^（l）:
0.01：
x（end）-'越=pgjy嗖血用】）；，
所找到的方程〔多项式）变换尸paly2£tr氯捏）;心
用pht画出所找到的罢项武图像•
边慝〔a1,a2八b；岳材.匕‘naarkErth:
匕2D”
五、Matlab程序源代码：
1、正交多项式序列的生成
functiongouzaozhengjiaoduoxiangshi
数值例题：
1）当区间为［-1,1］，权函数r（x）=1时，由｛1,x,...,xn,...｝正交化得到的多项式就成为勒让
1dn
德多项式，它的表达式为P0（x）=1,Pn（x）n-比｛（x2-1）n｝；
2n!
dx
递推公式为（n1）Pn1（x）=（2n1）xPn（x）-nPn4（x）,（n=1,2,...）,
由F>（x）=1,P（x）二x，利用
（1）就可推出
2
P2（x）=（3x-1）/2,
P3（x）=（5x3-3x）/2,
42
p4（x）=（35x-30x3）/8,
由序列｛1,x,…，xn,...｝正交化得到的多项式
就是切比雪夫多项式，它可表示为Tn（x）=cos（narccosx）,|x|空1.
若令x=cost，贝UTn（x）二cosnv,0「“：
：
二.
可推出
递推关系IxTxTndTnW，2，...），
To（X）—1,Ti（X）=X.
2
T2（x=2x-1）,
=4x3-3x,
T3（x）
2、最佳一致逼近多项式
functionyicibijin
_•-
数值例题：
（1）求f（x）=--1•x在［0,1］上的最佳一次逼近多项式
结果为:
R（x）=0.9550.414x；
误差限为
max|.1x?
-R（x）|乞0.045.0知
3、最佳平方逼近多项式
functionpingfang
数值例题：
设f（x）二.1x2，求［0,1］上的一次最佳平方逼近多项式；结果为：
S（X）=0.934+0.426X,||6（x）血=0.0026；
4、曲线拟合的最小二乘法：
functionp=nihe（x,y,m）
数值实例：
例1:
下面给定的是乌鲁木齐最近1个月早晨7:
00左右（新疆时间）的天气预报所得到的
温度数据表，按照数据找出任意次曲线拟合方程和它的图像•
（2006年10月26~11月26）
天数
1
2
3
4
5
6
7
8
9
10
温度
9
10
11
12
13
14
13
12
11
9
天数
11
12
13
14
15
16
17
18
19
20
温度
10
11
12
13
14
12
11
10
9
8
天数
21
22
23
24
25
26
27
28
29
30
温度
7
8
9
1
9
7
6
5
3
1
1
解：
x=1:
30
y=[91011121314131211910111213141211109878911876531]
m=9；（m是任意真整数,但是不要取的太大）
运行p=nihe（x,y,m）
结果为:
p=
图形为:
或者任意取部分数据也可得到相应得多项式和它的图像
m=15时
x=1:
30
y=[91011121314131211910111213141211109878911876531]运行p=nihe（x,y,m）
结果为：
p=
图形为：
如果:
m=3
x=1:
30
y=[91011121314131211910111213141211109878911876531]运行p=nihe（x,y,m）
结果为：
p=
图形为：
（m=length（x））跟X的数
如果写M文件的时候我们把想得到的多项式的次数直接定义为
量一样取得最小二乘法的曲线拟合方程，这样上术问题中的m=30,这时候它的程序代码为：
functionp=nn（x,y）
赋值x,y:
x=1:
30
y=[91011121314131211910111213141211109878911876531]；运行
p=nn（x,y）
解为：
P=
图像为：
结果分析：
以上结果可以看到用最小二乘拟合来求解问题时，有时候他的结果很接近实际情况，有
时候跟实际情况里的太远，因为所求得多项式次数太小时数据点之间差别很大•次数最大是
误差最小但是有时后不符合实际情况，所以用最小二乘法时次数要取合适一点
或采用：
functionp=duoxiangxi（x,y,m）
数值例题：
（1）合成纤维抽丝工段中第一导丝盘的速度对丝的质量有很大的影响，第
一丝盘的速度和电流周波有重要关系。
下面是一组实例数据：
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
其中x代表电流周波，y代表第一导丝盘的速度冃的：
研究第一导丝盘速度y与电流周波x的关系。
解：
根据所给数据，在matlab窗口运行:
x=[49.250.049.349.049.049.549.849.950.250.2]
y=[16.717.016.816.616.716.816.917.017.017.1]
plot（x,y,'p','markersize',15）
17.1
17.05
17
☆
1B.95
IS.S
1B.75
16.65
49.4
☆
50
SQ.2
90.4
S1（x）=a
从图中看到各点在一条直线附近，故可选择线性函数作拟合曲线，既令
故，根据曲线拟合的最小二乘法，公式（5.5），（5.6），即
10
（:
j,-k）八"j「k；
i=Q
（K=091，ooooo）
10
dk=（f,A）八讦「k；
（K=0,1,000o）
i£
这里，m=10,n=1,「二1，l（x）二1,I（x）二X,
扌巴x=[49.250.049.349.049.049.549.849.950.250.2]；
y=[16.717.016.816.616.716.816.917.017.017.1]；
代入公式（5.5），（5.6），
（%%），（申0,叫）=（久％），严汽），严0,f），（申1,f）分别可得:
10
（；：
0,；：
0）10；
i=1
1010
（：
0,1）=（1,：
0）八％=496.1，（：
1,1）Xi2=24614
i=1i=1
（dfj二'fi=168.6,（仆fj二、Xifi=8364.9
i4i4
得方程组
10ao496.1ai=168.6
496.1a024614a1=8364.9
解得：
a1=0.3389,a0=0.0490。
于是所求拟合曲线为：
S=0.33886x+0.048969
102
误差平方和：
|■■：
■（x）2=二i[s（Xi）-f（Xi）]=0.0469
i二
MATLAB上运行的结果：
x=[49.250.049.349.049.049.549.849.950.250.2];
y=[16.717.016.816.616.716.816.917.017.017.1];
m=1;
S=
0.33886x+0.04896
P=
0.33890.0490
程序运行的结果中S（x）=0.33886x+0.048969为所求拟合曲线的方程
下面看离散点x，f=y和所得的拟合曲线的逼近情况：
在matlab窗口输入
x=[49.250.049.349.049.049.549.849.950.250.2];
f=[16.717.016.816.616.716.816.917.017.017.1];
x仁49:
0.01:
50.4;
s='0.3389*x+0.0490';
s=inline（s）;
s1=s（x1）;
图二
误差平方和：
在matlab窗口运行以下程序
f=[16.717.016.816.616.716.816.917.017.017.1];t=0;
x=[49.250.049.349.049.049.549.849.950.250.2];
S=0.3389*x+0.0490;
n=length（x）;
fori=1:
n
t=t+（f（i）-S（i）F2;
end,t
结果：
t=0.0469
（2）已知一组实验数据如下，求它的拟合曲线。
Xi
1
2
3
4
5
fi
4
4.5
6
8
8.5
2
1
3
1
1
结果为：
S；（x）=2.771.13x；
如图所示：
e
-■+1
T.-S
F
£&
-
-
$
-
45
«■
丹
1
1U1
1
1
1
-
1n
25
1A
J
A
升级会员 