最小二乘法实验报告

【实验目的】：观察最小二乘多项式的数值不稳定现象

【实验内容】：1 在[-1,1]区间上取n=20个等距节点，计算出以相应节点上x e 的值做为数据样本，以21,,,,l x x x 为基函数作出3,5,7,9,11,13,15l =次的最小二乘多项式，画出

2ln(())cond A ～l 之间的曲线，其中A 是确定最小二乘多项式的系数矩阵。计算出不同阶

最小二乘多项式给出的最小误差2

1

()(())

n

i

i

i l y x y σ==

-∑

2 在[-1,1]区间上取n=20个等距节点，计算出以相应节点上x e 的值做为数据样本，以121,(),(),()l p x p x p x 为基函数作出3,5,7,9,11,13,15l =次的最小二乘多项式，其中，

()i p x 是勒让德多项式。画出2ln(())cond A ～l 之间的曲线，其中A 是确定最小二乘多项式

的系数矩阵。计算出不同阶最小二乘多项式给出的最小误差2

1

()(())

n

i

i

i l y x y σ==

-∑，把结

果与1比较

【实验步骤及结果】：在[-1,1]区间上取n=20个等距节点，步长h=2/19，计算出以相应节点上x

e 的值做为数据样本，数据如表格1。

表格 1 数据样本值

（1）

以21,,,,l

x x x 为基函数拟合x e

在matlab 中编写函数lsmex （x, y, l ）,生成最小二乘法的系数矩阵A 、右端向量d ，求

出系数a =[a 0,a 1,a 2,…,a l ]T =A −1d ，得不同阶数下的最小二乘多项式y l x = a i x l

l i =0=a T X ，其中，X =[1,x ,x 2,…,x l ]T 。计算系数a 的结果如下：

0.9955486716942

0.99757900166

0.540351495874469

0.176998749075551

②l=5，a=

1.000038552490093

1.000020023273578

0.499246808308251

0.1697111425869

0.043753810227387

0.0086818882249

③l=7，a=

0.999999832859033

0.999999916830632

0.500005765944358

0.166667853917488

0.0416331725214

0.008328850726684

0.0014385132655

0.0002045690792

④l=9，a=

1.000000000409680

1.000000000193609

0.499999977559871

0.166666662393254

0.041666859391586

0.008333359277458

0.001388319252092

0.000198349327547

0.000025478122502

0.000002822438546

⑤l=11，a=

1.000000000008706

1.000000000002729

0.500000000050966

0.166666666671517

0.041666666004638

0.008333333185874

0.0013881754987

0.000198411522433

0.000024795348762

0.000002756249160

0.000000281543798

0.000000025378540 ⑥ l=13，a = 0.999999999770700 0.999999998584826 0.500000000002793 0.166666666744277 0.041666666993633 0.008333334699273 0.0013884421398 0.000198394060135 0.0000247998846 0.000002682209015 0.000000272755242 0.000000014901161 0.000000003426662 0 ⑦ l=15，a = 1.000000006592917 1.000000012805685 0.500000000015796 0.1666666183140 0.041666662936834 0.008333325386047 0.001388842803919

0.000197768211365

0.000025072928111 0.000001430511475 0.000000863215519 0.000001907348633 0.000000217004981 -0.000002384185791 0.000000011572638

0.000000238418579

计算出不同阶最小二乘多项式的误差并比较得到最小误差，最后计算cond (A )2，

绘出2ln(())cond A ～l 之间的曲线如图1，拟合误差与阶数的关系曲线如图2。

表格 2 不同阶数下的系数矩阵条件数和拟合误差值

最小误差σ 11

1.191342440055042e-018

图 1 系数矩阵A 的谱条件数和阶数的关系

图 2 拟合误差与阶数的关系

（2）以121,(),(),()l p x p x p x 为基函数拟合x

e

在matlab 中编写函数lsmexlgd （x, y, l ）,生成最小二乘法的系数矩阵A 、右端向量d ，求出系数a =[a 0,a 1,a 2,…,a l ]T =A −1d ，得不同阶数下的最小二乘多项式y l x =

a i p l (x )l

i =0=a T X ，其中，X =[1,p 1(x ),p 2(x ),…,p l (x )]T 。计算系数a 的结果如下：

①l =3，a =

1.175666152008432 1.103778251135020 0.360234330582979

0.070799499630219 ② l =5，a = 1.1752049173049 1.103639099368339 0.357833382811630 0.070457461559185 0.010000870909119 0.001102461996817

③ l =7，a =

1.175201209747137 1.1036383263283 0.357814********* 0.070455639883299 0.0099652650863 0.001099594726534 0.000099637282407 0.000007629616011 ④

l =9，a =

1.175201193697940 1.103638323523287 0.357814********* 0.070455633687156 0.009965128554402 0.001099586149963 0.000099454781535 0.000007620559433 0.000000506790883 0.000000029722041

⑤ l =11，a =

1.1752011933971 1.103638323514353 0.3578143508158 0.070455633668543 0.009965128150024 0.001099586127269 0.000099454340203 0.000007620541351

0.0000005072467 0.000000029718138 0.000000001560386

0.000000000074146 ⑥

l =13，a =

1.1752011933802 1.103638323514325 0.3578143507374 0.070455633668488 0.009965128148872 0.001099586127206 0.000099454339116 0.000007620541308 0.0000005071976 0.000000029718141 0.000000001560885 0.000000000074202 0.000000000003220 0.000000000000132 ⑦ l =15，a = 1.1752011933800

1.103638323514317 0.3578143507367 0.070455633668470

0.009965128148861

0.001099586127186 0.000099454339106 0.000007620541295 0.0000005071972 0.000000029718142 0.0000000015600 0.000000000074217 0.000000000003233 0.000000000000154

0.000000000000015

0.000000000000020

计算出不同阶最小二乘多项式的误差并比较得到最小误差，最后计算cond (A )2，绘出

2ln(())cond A ～l 之间的曲线如图3，拟合误差与阶数的关系曲线如图4。

表格 3 不同阶数下的系数矩阵条件数和拟合误差值

最小误差σ13 1.528726152656812e-029

图 3 系数矩阵A的谱条件数和阶数的关系

图 4 拟合误差与阶数的关系

【对比与分析】图5为两种方法的系数矩阵条件数与阶数的关系曲线对比，从图看出，用勒让德多项式拟合原函数（方法一）时系数矩阵条件数大大小于以21,,,,l x x x 为基函数拟合原函数（方法二）时系数矩阵的条件数。系数矩阵条件数越大，系数的线性方程组的解误差越大。当l =15时，方法一求出的系数a 出现负值。因此 ，从系数矩阵条件数考虑 ，方法一优于方法二。

图 5 两种方法系数矩阵条件数与阶数的关系曲线对比

图 6 两种方法拟合误差与阶数的关系曲线对比

方法一

方法二

方法一

方法二

图6为两种方法的最小拟合误差与阶数关系曲线对比。在阶数 l ≤8时，两种方法的最小拟合误差相同，当l >8时，方法一的最小拟合误差远远大于方法二的最小拟合误差，且误差逐渐增大，出现误差不稳定现象。对于方法二，当l >13时误差开始增大，出现不稳定现象。因此，从最小拟合误差考虑，方法二优于方法一。

综上，在[-1,1]上拟合函数x e 时，以勒让德多项式为基函数比以21,,,,l x x x 为基函数更好。

附函数源代码：

（1）以21,,,,l x x x 为基函数拟合

function lsmex(x,y,l)

%给定一组数据，以span{1,x,x^2,x^3,...,x^l;l=0,1,2,3,...}基函数，选取不同的l

%用最小二乘法拟合exp （x ）,确定相应的最小二乘法系数矩阵A ，画出ln(Cond(A))-l 的 %曲线，计算最小误差delta ，绘制拟合误差与阶数l 的关系曲线。 for m = 1:7 t(m) = l(m)+1 C = ones(20,t(m)) for i = 1:20

for j = 1:t(m)

C(i,j)=x(i)^(j-1)

end end

A = C'*C %生成最小二乘法系数矩阵A d = C'*y' %生成右端向量d a = inv(A)*d %计算系数a %求拟合值yy yy0 = ones(20,t(m)) for i = 1:20

for j = 1:t(m) yy0(i,j) = x(i)^(j-1) end end

yy = yy0*a %求误差delta delta = 0 for k = 1:20

delta = delta +(yy(k)-y(k))^2 end

deltaA(m) = delta

condA(m) = cond(A,2) %计算A 的条件数 end

%求最小误差mindelta mindelta = min(deltaA)

figure(1)

plot(l,log(condA),'k-.o','LineWidth',2)

grid

title('系数矩阵A的谱条件数与阶数的关系')

xlabel('阶数l')

ylabel('ln(cond(A))')

%绘制拟合误差曲线

figure(2)

plot(l,log(deltaA),'k-.p','LineWidth',2)

grid

title('拟合误差与阶数的关系')

xlabel('阶数l')

ylabel('ln(delta)')

end

（2）以勒让德多项式为基函数拟合

function lsmexlgd(x,y,l)

%给定一组数据，以span{1,p1(x),p2(x),p3(x),...,pl(x);l=0,1,2,3,...}为基函数

%（其中p(x)为Legendre多项式），

%选取不同的l用最小二乘法拟合exp（x）,确定相应的最小二乘法系数矩阵A，画出%ln(Cond(A))-l的曲线，并计算最小误差delta，绘制拟合误差与阶数l的曲线。

for m = 1:7

t(m) = l(m)+1

for i = 1:20

b = LegendreP(x(i),t(m))

for j = 1:t(m)

C(i,j) = b(j)

end

end

A = C'*C %生成最小二乘法系数矩阵A

d = C'*y' %生成右端向量d

a = A\\d %计算系数a

condA(m) = cond(A,2)

%求拟合值yy

yy = C*a

%求误差delta

delta = 0

for k = 1:20

delta = delta +(yy(k)-y(k))^2

end

deltaA(m) = delta

condA(m) = cond(A,2) %计算A的条件数

end

%求最小误差mindeltamindelta = min(deltaA)

%绘制ln(Cond(A))-l的曲线

figure(1)

plot(l,log(condA),'k-.o','LineWidth',2)

grid

title('系数矩阵A的谱条件数与阶数的关系')

xlabel('阶数l')

ylabel('ln(cond(A))')

%绘制拟合误差曲线

figure(2)

plot(l,log(deltaA),'k-.p','LineWidth',2)

grid

title('拟合误差与阶数的关系')

xlabel('阶数l')

ylabel('ln(delta)')

end

其中LegendreP为勒让德多项式函数，源代码如下：

function LegendreMATRIX=LegendreP(X,N)

if N<2

fprintf('N smaller than needed parameter\\n') ;%可能缺少一个跳出指令end

LegendreMATRIX=ones(size(X,2),N);%生成Legendre多项式矩阵

if N>=2

LegendreMATRIX(:,2)=X';

end

if N>=3

for i=3:N

LegendreMATRIX(:,i)=[(2*(i-1)-1)*X(1,:)'.*LegendreMATRIX(:,i-1)-(i-2) *LegendreMATRIX(:,i-2)]/(i-1);

%注意此时的阶次与矩阵序列号的关系，i与i-1的区别

end

end

end