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

MATLAB 程序设计实践课程考核 一、编程实现以下科学计算算法，并举一例应用之。（参考书籍精通科学计算，王正林等著，电子工业出版社，年）“正交多项式最小二乘法拟合”正交多项式最小二乘法拟合原理正交多项式做最小二乘法拟合：不要求拟合函数y=f(x)经过所有点（xi,yi)，而只要求在给定点xi上残差i=f(xi)-yi按照某种标准达到最小，通常采用欧式范数|2作为衡量标准。这就是最小二乘法拟合。（1）根据作为给定节点x0，x1，xm及权函数(x)0,造出带权函数正交的多项式Pn（x）。注意nm,用递推公式表示Pk（x），即这里的Pk(x)是首项系数为1的k次多项式，根据Pk(x)的正交性，得（2）根据公式(1)和(2)逐步求Pk（x）的同时，相应计算系数并逐步把 Pk（x）累加到S（x）中去，最后就可得到所求的拟合函数曲线流程图开始读取点集x,y和n数据size(x)=size(y)N提示x,y维数不匹配Y利用x点集数据构造范德蒙德系数矩阵V调用QR分解函数求的多项式系数输出多项式系数P结束M文件function p = mypolyfit(x,y,n)%定义mypolyfit为最小二乘拟合函数%P = POLYFIT(X,Y,N)以计算以下多项式系数%P(1)*XN + P(2)*X(N-1) +.+ P(N)*X + P(N+1).if isequal(size(x),size(y) 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 pointsn . or try centering and scaling as described in HELP POLYFIT.) endendr = y - V*p;p = p.; % 将多项式系数默认为行向量.5、运行流程图 调用mypolyfit.m进行运算根据所得到的多项式系数构造拟合函数结束输入变量x=0.5000 1.00001.50002.00002.50003.0000；y=1.75 2.45 3.81 4.80 8.00 8.60开始过程：clearx = 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000y=1.75 2.45 3.81 4.80 8.00 8.60x1=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 onplot(x1,y1,r)2、 编程计算以下电路问题例8-1-3如图所示电路，已知R=5，L=3，=5,Uc=10,求R,C,和L，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(IR IC I UL US)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、 编程解决以下问题求自然三次样条曲线，经过点（-3,2），（-2,0），（1,3），（4,1），而且自由边界条件(-3)=0,(4)=0。算法说明：三次样条也是分片三次插值函数，它是在给定的区间a，b上的一个划分：a=x0x1=x0(j)&(x(k)x0(j+1) l(k)=j; end endendfor k=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(