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

《MATLAB程序设计实践》课程考核一、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“正交多项式最小二乘法拟合”正交多项式最小二乘法拟合原理正交多项式做最小二乘法拟合：不要求拟合函数y=f(x)经过所有点（xi,yi)，而只要求在给定点xi上残差δi=f(xi)-yi按照某种标准达到最小，通常采用欧式范数||δ||2作为衡量标准。这就是最小二乘法拟合。根据作为给定节点x0，x1，…xm及权函数ρ(x)0,造出带权函数正交的多项式{Pn（x）}。注意n≤m,用递推公式表示Pk（x），即01101111,,(1,2,,1)kkkkkPxPxxPxPxPxPxkn这里的Pk(x)是首项系数为1的k次多项式，根据Pk(x)的正交性，得2i012i02i0211i10x,,x,0,1,1,x,0,1,1,xmikikkikmkkkiikkmkkkikkikmkkkiixPxxPxPxaPxPxPxxPPknPPPxPPknPPPx根据公式(1)和(2)逐步求Pk（x）的同时，相应计算系数020()(),(0,1,n(,)()mijikikikmkkikiixxxfPakPPxx）并逐步把*kaPk（x）累加到S（x）中去，最后就可得到所求的拟合函数曲线***0011nny=Sx=aPx+aPx++aPx（）（）（）（）.流程图（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',...'XandYvectorsmustbethesamesize.')end%检验XY维数是否匹配x=x(:);y=y(:);ifnargout2mu=[mean(x);std(x)];x=(x-mu(1))/mu(2);end%利用范德蒙德矩阵构造方程组系数矩阵V(:,n+1)=ones(length(x),1,class(x));forj=n:-1:1V(:,j)=x.*V(:,j+1);end%对矩阵进行QR分解以求得多项式系数值读取点集x,y和n数据size(x)=size(y)NY输出多项式系数P结束提示x,y维数不匹配利用x点集数据构造范德蒙德系数矩阵V调用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.0e10ifnargout2warning('MATLAB:polyfit:RepeatedPoints',...'Polynomialisbadlyconditioned.Removerepeateddatapoints.')elsewarning('MATLAB:polyfit:RepeatedPointsOrRescale',...['Polynomialisbadlyconditioned.Removerepeateddatapoints\n'...'ortrycenteringandscalingasdescribedinHELPPOLYFIT.'])endendr=y-V*p;p=p.';%将多项式系数默认为行向量.5、运行流程图过程：clearx=[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,'*')holdonplot(x1,y1,'r')调用mypolyfit.m进行运算根据所得到的多项式系数构造拟合函数结束输入变量x=[0.50001.00001.50002.00002.50003.0000]；y=[1.752.453.814.808.008.60]开始二、编程计算以下电路问题[例8-1-3]如图所示电路，已知R=5Ω，ωL=3Ω，C1=5Ω,Uc=1000,求I.R,I.C,I.和U.L，U.S，并画其相量图。理论分析：根据电路分析Z=R+j*(Xl-Xc)Ic=Uc/Z3;Z3=-j*XcIr=Ur/Z2=Uc/Z2;Z2=RI=Ir+IcUl=I*Z1;Z1=j*XLUs=Ul+Ur计算得Ir=2；Ic=2.00iI=2.00+2.00iUl=-6.00+6.00iUs=4.00+6.00iM文件clearR=5;XL=3;XC=5;UC=10;UR=UC;%为给定元件赋值Z1=j*XL;Z2=R;Z3=-j*XC;%定义各电抗disp('电流')IR=UR/Z2%计算IRIC=UC/Z3%计算ICI=IR+IC%计算Idisp('电压')UL=I*Z1%计算ULUS=UC+UL%计算USdisp('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)运行流程图：运行图开始读取已知数据构造电抗Z1，Z2，Z3计算Ir,Ic,I,Ul,UsIR=UR/Z2；IC=UC/Z3；I=IR+IC；UL=I*Z1；US=UC+UL结束输出Ir,Ic,I,Ul,Us并画出对应的相量图三、编程解决以下问题求自然三次样条曲线，经过点（-3,2），（-2,0），（1,3），（4,1），而且自由边界条件S''(-3)=0,S''(4)=0。算法说明：三次样条也是分片三次插值函数，它是在给定的区间[a，b]上的一个划分：a=x0x1…xn=b,已知函数f(x)在xj上的函数值为f(xj)=yj，（j=0,1,2,3...,n)如果存在分段函数nnnxxxxSxxxxSxxxxSxS,)(,)(,)()(1212101满足下述条件：（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）由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)进行循环，可得3n-2个方程S(xj)=yjj=0,1,2,3...,n得n个方程加上自由边界条件可得2个方程。共4n个方程，可唯一确定Sj(x)=ajx3+bjx2+cjx+dj方程中S(x)的表达式中需确定的4n个未知数{aj,bj,cj,dj}由确定的{aj,bj,cj,dj}得出方程的解Sj(x)开始储存结果（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个插值区间的长度endforj=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:1m(j)=a(j)*m(j+1)+b(j);endfork=1:kmforj=1:(n-1)if((x(k)=x0(j))&(x(k)x0(j+1)))l(k)=j;end用myspline的插值方法进行插值处理开始结束输入原始数据输出插值文件输出图形endendfork=1:kmsum=(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)holdonplot(x,y,'or')legend('计算