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

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

《正交多项式最小二乘法拟合.docx》由会员分享，可在线阅读，更多相关《正交多项式最小二乘法拟合.docx（16页珍藏版）》请在冰豆网上搜索。

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

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

一、编程实现以下科学计算算法，并举一例应用之。

（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）

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

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

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

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

这就是最小二乘法拟合。

根据作为给定节点x0，x1，…xm及权函数ρ（x）>0,造出带权函数正交的多项式{Pn（x）}。

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

这里的Pk（x）是首项系数为1的k次多项式，根据Pk（x）的正交性，得

根据公式

（1）和

（2）逐步求Pk（x）的同时，相应计算系数

并逐步把

Pk（x）累加到S（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））

error（'MATLAB:

polyfit:

XYSizeMismatch',...

'XandYvectorsmustbethesamesize.'）

end

%检验XY维数是否匹配

x=x（:

）;

y=y（:

）;

ifnargout>2

mu=[mean（x）;std（x）];

x=（x-mu

（1））/mu

（2）;

end

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

V（:

n+1）=ones（length（x）,1,class（x））;

forj=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）;

ifsize（R,2）>size（R,1）

warning（'MATLAB:

polyfit:

PolyNotUnique',...

'Polynomialisnotunique;degree>=numberofdatapoints.'）

elseifcondest（R）>1.0e10

ifnargout>2

warning（'MATLAB:

polyfit:

RepeatedPoints',...

'Polynomialisbadlyconditioned.Removerepeateddatapoints.'）

else

warning（'MATLAB:

polyfit:

RepeatedPointsOrRescale',...

['Polynomialisbadlyconditioned.Removerepeateddatapoints\n'...

'ortrycenteringandscalingasdescribedinHELPPOLYFIT.']）

end

end

r=y-V*p;

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

5、运行流程图

过程：

clear

x=[0.50001.00001.50002.00002.50003.0000]

y=[1.752.453.814.808.008.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,'*'）

holdon

plot（x1,y1,'r'）

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

[例8-1-3]如图所示电路，已知R=5Ω，ωL=3Ω，

=5Ω,Uc=10

求

R,

C,

和

L，

S，并画其相量图。

理论分析：

根据电路分析Z=R+j*（Xl-Xc）

Ic=Uc/Z3;Z3=-j*Xc

Ir=Ur/Z2=Uc/Z2;Z2=R

I=Ir+Ic

Ul=I*Z1;Z1=j*XL

Us=Ul+Ur

计算得

Ir=2；

Ic=2.00i

I=2.00+2.00i

Ul=-6.00+6.00i

Us=4.00+6.00i

M文件

clear

R=5;XL=3;XC=5;UC=10;UR=UC;%为给定元件赋值

Z1=j*XL;Z2=R;Z3=-j*XC;%定义各电抗

disp（'电流'）

IR=UR/Z2%计算IR

IC=UC/Z3%计算IC

I=IR+IC%计算I

disp（'电压'）

UL=I*Z1%计算UL

US=UC+UL%计算US

disp（'IRICIULUS'）

disp（'幅值'）;disp（abs（[IR,IC,I,UL,US]））

disp（'相角'）;disp（angle（[IR,IC,I,UL,US]）*180/pi）

ph=compass（[IR,IC,I,UL,US]）;%分别画出IR,IC,I,UL,US相量图

set（ph,'linewidth',3）

运行流程图：

运行图

三、编程解决以下问题

求自然三次样条曲线，经过点（-3,2），（-2,0），（1,3），（4,1），而且自由边界条件

（-3）=0,

（4）=0。

算法说明：

三次样条也是分片三次插值函数，它是在给定的区间[a，b]上的一个划分：

a=x0

满足下述条件：

（1）S（x）在每一个子区间[xj-1，xj]，（j=0,1,2,3...,n）上是三次多项式；

（2）S（x）在每一个内接点（j=0,1,2,3...,n）具有直到二阶的连续导数；

则成为节点x0，x1，…，xn上的三次样条函数。

若S（x）在节点x0，x1，…，xn上海满足插值条件：

（3）S（xj）=yj（j=0,1,2,3...,n）

则称S（x）为三次样条插值函数。

由

（1）知，S（x）在每一个小区间[xj-1，xj]上是一三次多项式，若记为Sj（x），则可设：

Sj（x）=ajx3+bjx2+cjx+dj

要确定函数S（x）的表达式，需确定4n个未知数{aj,bj,cj,dj}（j=0,1,2,3...,n）。

由

（2）知S（x），S'（x），S''（x）在内节点x1,x2,....,xn-1上连续，则

Ⅰ型S（xj-0）=S（xj+0）

Ⅱ型S'（xj-0）=S'（xj+0）

Ⅲ型S''（xj-0）=S''（xj+0）j=0,1,2,3...,n-1

可得3n-2个方程，又由条件（3）S（xj）=yjj=0,1,2,3...,n

得n个方程，共得到4n-2个方程。

要确定4n个未知数，还差俩个方程。

通常在端点x0=a.xn=b处各附加一个条件

称边界条件，常见有三种：

（1）自然边界条件：

S''（x0）=S''（xn）=0

（2）固定边界条件：

S'（x0）=f'（x0）,S'（x）=f'（xn）

（3）周期边界条件：

S'（x0）=S'（xn），S''（x0）=S''（xn）

共4n个方程，可唯一的确定4n个未知数

流程图：

（1）

（2）运行流程图

源代码：

functions=myspline2（x0,y0,y21,y2n,x）

%s=myspline（x0,y0,y21,y2n,x）

%x0,y0为已知插值点,y21,y2n为二型边界条件

n=length（x0）;

km=length（x）;

a

（1）=-0.5;

b

（1）=3*（y0

（2）-y0

（1））/（2*（x0

（2）-x0

（1）））;

forj=1:

（n-1）

h（j）=x0（j+1）-x0（j）;%h（j）为第j个插值区间的长度

end

forj=2:

（n-1）

alpha（j）=h（j-1）/（h（j-1）+h（j））;

beta（j）=3*（（1-alpha（j））*（y0（j）-y0（j-1））/h（j-1）+alpha（j）*（y0（j+1）-y0（j））/h（j））;

a（j）=-alpha（j）/（2+（1-alpha（j））*a（j-1））;

b（j）=（beta（j）-（1-alpha（j））*b（j-1））/（2+（1-alpha（j））*a（j-1））;

end%alpha（j）,beta（j）为插值基函数

m（n）=（3*（y0（n）-y0（n-1））/h（n-1）+y2n*h（n-1）/2-b（n-1））/（2+a（n-1））;

forj=（n-1）:

-1:

1

m（j）=a（j）*m（j+1）+b（j）;

end

fork=1:

km

forj=1:

（n-1）

if（（x（k）>=x0（j））&（x（k）

l（k）=j;

end

end

end

fork=1:

km

sum=（3*（x0（l（k）+1）-x（k））^2/h（l（k））^2-2*（x0（l（k）+1）-x（k））^3/h（l（k））^3）*y0（l（k））;

sum=sum+（3*（x（k）-x0（l（k）））^2/h（l（k））^2-2*（x（k）-x0（l（k）））^3/h（l（k））^3）*y0（l（k）+1）;

sum=sum+h（l（k））*（（x0（l（k）+1）-x（k））^2/h（l（k））^2-（x0（l（k）+1）-x（k））^3/h（l（k））^3）*m（l（k））;

s（k）=sum-h（l（k））*（（x（k）-x0（l（k）））^2/h（l（k））^2-（x（k）-x0（l（k）））^3/h（l（k））^3）*m（l（k）+1）;

end

举例：

x=[-3-214]

y=[2031]

x0=[-3:

0.15:

4];

y0=myspline（x,y,0,0,x0）;

plot（x0,y0）

holdon

plot（x,y,'or'）

legend（'计算值','实验值'）