正交多项式最小二乘法拟合Word文档下载推荐.docx
《正交多项式最小二乘法拟合Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《正交多项式最小二乘法拟合Word文档下载推荐.docx（15页珍藏版）》请在冰豆网上搜索。
MATLAB:
polyfit:
XYSizeMismatch'
...
'
XandYvectorsmustbethesamesize.'
）
end
%检验XY维数是否匹配
x=x（:
）;
y=y（:
ifnargout>
2
mu=[mean（x）;
std（x）];
x=（x-mu
（1））/mu
（2）;
%利用范德蒙德矩阵构造方程组系数矩阵
V（:
n+1）=ones（length（x）,1,class（x））;
forj=n:
-1:
1
V（:
j）=x.*V（:
j+1）;
%对矩阵进行QR分解以求得多项式系数值
[Q,R]=qr（V,0）;
ws=warning（'
off'
'
all'
p=R\（Q'
*y）;
warning（ws）;
ifsize（R,2）>
size（R,1）
warning（'
PolyNotUnique'
...
Polynomialisnotunique;
degree>
=numberofdatapoints.'
elseifcondest（R）>
1.0e10
ifnargout>
RepeatedPoints'
Polynomialisbadlyconditioned.Removerepeateddatapoints.'
else
RepeatedPointsOrRescale'
['
Polynomialisbadlyconditioned.Removerepeateddatapoints\n'
...
ortrycenteringandscalingasdescribedinHELPPOLYFIT.'
]）
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
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
电压'
UL=I*Z1%计算UL
US=UC+UL%计算US
IRICIULUS'
幅值'
disp（abs（[IR,IC,I,UL,US]））
相角'
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<
x1<
…<
xn=b,已知函数f（x）在xj上的函数值为f（xj）=yj，（j=0,1,2,3...,n）如果存在分段函数
满足下述条件：
（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）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）固定边界条件：
（x0）=f'
（x0）,S'
（x）=f'
（xn）
（3）周期边界条件：
（xn），S'
共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个插值区间的长度
forj=2:
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）:
m（j）=a（j）*m（j+1）+b（j）;
fork=1:
km
forj=1:
if（（x（k）>
=x0（j））&
（x（k）<
x0（j+1）））
l（k）=j;
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）;
举例：
x=[-3-214]
y=[2031]
x0=[-3:
0.15:
4];
y0=myspline（x,y,0,0,x0）;
plot（x0,y0）
holdon
or'
legend（'
计算值'
实验值'
升级会员 