正交多项式最小二乘法拟合.docx
《正交多项式最小二乘法拟合.docx》由会员分享，可在线阅读，更多相关《正交多项式最小二乘法拟合.docx（15页珍藏版）》请在冰豆网上搜索。
正交多项式最小二乘法拟合
《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（'计算值','实验值'）
升级会员 