正交多项式最小二乘法拟合

《MATLAB 程序设计实践》课程考核

一、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）

“正交多项式最小二乘法拟合”

正交多项式最小二乘法拟合原理

正交多项式做最小二乘法拟合：

不要求拟合函数y=f(x)经过所有点（x i,y i)，而只要求在给定点x i上残差δi=f(x i)-y i按照某种标准达到最小，通常采用欧式范数||δ||2作为衡量标准。这就是最小二乘法拟合。根据作为给定节点x0，x1，…x m及权函数ρ(x)>0,造出带权函数正交的多项式{P n（x）}。注

意n≤m,用递推公式表示P k（x），即

()

()()()

()()()

110

111

1,

,

(1,2,,1) k k k k k

P x

P x x P x

P x P x P x k n

α

αβ

++-

=

⎧

⎪

=-

⎨

⎪=--=...-

⎩

这里的P k(x)是首项系数为1的k次多项式，根据P k(x)的正交性，得

()()

()()

()()

()

()()

()

()()

()()

()

()

()

()

()

()

2

i

1

2

i

2

i

211

i1

x

,

,

x

,

0,1,1

,

x

,

0,1,1

,

x

m

i k i

k k

i

k m

k k

k i

i k k

m

k k

k i

k k

i

k m

k k

k i

i

x P x

xP x P x

a

P x P x

P x

xP P

k n

P P

P x

P P

k n

P P

P x

ρ

ρ

ρ

β

ρ

=

+

=

=

--

-

=

⎧

⎪

⎪==

⎪

⎪⎪

==⋅⋅⋅-

⎨

⎪

⎪

===⋅⋅⋅-

⎪

⎪

⎪⎩

∑

∑

∑

∑

根据公式(1)和(2)逐步求P k（x）的同时，相应计算系数

()()

()

2

()()

,

(0,1,n

(,)

()

m

i j i k i

k i

k m

k k

i k i

i

x x x

f P

a k

P P

x x

ρϕϕ

ρϕ

=

=

===⋅⋅⋅,

∑

∑

）

并逐步把*

k

a P k（x）累加到S（x）中去，最后就可得到所求的拟合函数曲线

***

0011n n

y=S x=a P x+a P x++a P x

⋅⋅⋅

（）（）（）（）.

流程图

（2）

（1）

M文件

function [p] = mypolyfit(x,y,n)

%定义mypolyfit为最小二乘拟合函数

%P = POLYFIT(X,Y,N)以计算以下多项式系数

%P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).

if ~isequal(size(x),size(y))

error('MATLAB:polyfit:XYSizeMismatch',...

'X and Y vectors must be the same size.') end

%检验X Y维数是否匹配

x = x(:);

y = y(:);

if nargout > 2

mu = [mean(x); std(x)];

x = (x - mu(1))/mu(2);

end

%利用范德蒙德矩阵构造方程组系数矩阵

V(:,n+1) = ones(length(x),1,class(x));

for j = n:-1:1

V(:,j) = x.*V(:,j+1);

end

% 对矩阵进行QR分解以求得多项式系数值

[Q,R] = qr(V,0);

ws = warning('off','all');

p = R\(Q'*y);

warning(ws);

if size(R,2) > size(R,1)

warning('MATLAB:polyfit:PolyNotUnique', ...

'Polynomial is not unique; degree >= number of data points.') elseif condest(R) > 1.0e10

if nargout > 2

warning('MATLAB:polyfit:RepeatedPoints', ...

'Polynomial is badly conditioned. Remove repeated data points.')

else

warning('MATLAB:polyfit:RepeatedPointsOrRescale', ...

['Polynomial is badly conditioned. Remove repeated data points\n'...

' or try centering and scaling as described in HELP POLYFIT.'])

end

end

r = y - V*p;

p = p.'; % 将多项式系数默认为行向量.

5、运行流程图