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

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

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

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

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

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

不要求拟合函数y=f(x)经过所有点（x i ,y i ),而只要求在给定点x i 上残差δi=f(x i )-y i 按照某种标准达到最小,通常采用欧式范数||δ||2作为衡量标准。这就是最小二

乘法拟合。

根据作为给定节点x 0,x 1,…x m 及权函数ρ(x)>0,造出带权函数正交的多项式{P n （x ）}。注意n ≤m,用递推公式表示P k （x ）,即

()()()()()()()01101

111,,

(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)的正交性,得

()()()()()()()()()()()()()()()()()()()()2i 012i 02i 0211i 10x ,,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 ）的同时,相应计算系数

()()()020()(),(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 ⋅⋅⋅（）（）（）（）.

流程图

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))

（2）

（1）

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、运行流程图

过程：

clear

x =[ 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000]

y=[1.75 2.45 3.81 4.80 8.00 8.60]

x1=0.5:0.05:3.0;

p=mypolyfit(x,y,2)

y1=p(3)+p(2)*x1+p(1)*x1.^2;

plot(x,y,'*')

hold on

plot(x1,y1,'r')

二、编程计算以下电路问题

[例8-1-3]如图所示电路,已知R=5Ω,ωL=3Ω,C 1ω=5Ω,Uc=100∠0,求I .R ,I .C ,I .和U .L ,U .

S ,并画其相量图。