线性最小二乘法拟合.doc
MATLAB程序设计实践课程考核一 实现以下科学计算算法，并举一例应用之。（参考书籍精通MATLAB科学计算，王正林等著，电子工业出版社，2009年）“线性最小二乘法拟合”解答：一， 最小二乘法拟合；方法介绍在科学实验的统计方法研究中，往往要从一组实验数据（xi，yi）中寻找出自变量x和因变量y之间的关系yf（x）。由于观测数据往往不准确，因此并不要求yf（x）经过所有的点（xi，yi），而只要求在给定的点xi上误差ðf（xi）yi按照某种标准达到最小，通常采用欧氏范数ð2作为误差量度的标准。这就是所谓的最小二乘法。MATLAB实现在MATLAB中实现最小二乘法拟合通常可以采用如下途径：利用polyfit函数进行多项式拟合。polyfit函数源程序：function p,S,mu = polyfit(x,y,n)%POLYFIT Fit polynomial to data.% P = POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of% degree N that fits the data Y best in a least-squares sense. P is a% row vector of length N+1 containing the polynomial coefficients in% descending powers, P(1)*XN + P(2)*X(N-1) +.+ P(N)*X + P(N+1).% P,S = POLYFIT(X,Y,N) returns the polynomial coefficients P and a% structure S for use with POLYVAL to obtain error estimates for% predictions. S contains fields for the triangular factor (R) from a QR% decomposition of the Vandermonde matrix of X, the degrees of freedom% (df), and the norm of the residuals (normr). If the data Y are random,% an estimate of the covariance matrix of P is (Rinv*Rinv')*normr2/df,% where Rinv is the inverse of R.% P,S,MU = POLYFIT(X,Y,N) finds the coefficients of a polynomial in% XHAT = (X-MU(1)/MU(2) where MU(1) = MEAN(X) and MU(2) = STD(X). This% centering and scaling transformation improves the numerical properties% of both the polynomial and the fitting algorithm.% Warning messages result if N is >= length(X), if X has repeated, or% nearly repeated, points, or if X might need centering and scaling.% Class support for inputs X,Y:% float: double, single% See also POLY, POLYVAL, ROOTS.% Copyright 1984-2004 The MathWorks, Inc.% $Revision: 5.17.4.5 $ $Date: 2004/07/05 17:01:37 $% The regression problem is formulated in matrix format as:% y = V*p or% 3 2% y = x x x 1 p3% p2% p1% p0% where the vector p contains the coefficients to be found. For a% 7th order polynomial, matrix V would be:% V = x.7 x.6 x.5 x.4 x.3 x.2 x ones(size(x);if isequal(size(x),size(y) error('MATLAB:polyfit:XYSizeMismatch',. 'X and Y vectors must be the same size.')endx = x(:);y = y(:);if nargout > 2 mu = mean(x); std(x); x = (x - mu(1)/mu(2);end% Construct Vandermonde matrix.V(:,n+1) = ones(length(x),1,class(x);for j = n:-1:1 V(:,j) = x.*V(:,j+1);end% Solve least squares problem.Q,R = qr(V,0);ws = warning('off','all'); p = R(Q'*y); % Same as p = Vy;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.' % Polynomial coefficients are row vectors by convention.% S is a structure containing three elements: the triangular factor from a% QR decomposition of the Vandermonde matrix, the degrees of freedom and% the norm of the residuals.S.R = R;S.df = length(y) - (n+1);S.normr = norm(r);流程图：例题设yspan1,x,x2,用最小二乘法拟合如表所示数据。X0.51.01.52.02.53.0y1.752.453.814.808.008.60开始用polyfit功能函数进行拟合。流程图输入数据调用polyfit函数拟合终止建立m文件源代码：function task1()x=0.5:0.5:3.0;y=1.75 2.45 3.81 4.80 8.00 8.60;a=polyfit(x,y,2)x1=0.5:0.05:3.0y1=a(3)+a(2)*x1+a(1)*x1.2plot(x,y,'*')hold onplot(x1,y1,'-r')legend('实验值','拟合值')end在matlab命令窗口中输入：task1()运行结果：task1()a = 0.4900 1.2501 0.8560x1 = Columns 1 through 9 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 Columns 10 through 18 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 1.3000 1.3500 Columns 19 through 27 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 Columns 28 through 36 1.8500 1.9000 1.9500 2.0000 2.0500 2.1000 2.1500 2.2000 2.2500 Columns 37 through 45 2.3000 2.3500 2.4000 2.4500 2.5000 2.5500 2.6000 2.6500 2.7000 Columns 46 through 51 2.7500 2.8000 2.8500 2.9000 2.9500 3.0000y1 = Columns 1 through 9 1.6036 1.6918 1.7825 1.8756 1.9712 2.0692 2.1697 2.2726 2.3780 Columns 10 through 18 2.4859 2.5961 2.7089 2.8241 2.9417 3.0618 3.1843 3.3093 3.4367 Columns 19 through 27 3.5666 3.6989 3.8337 3.9709 4.1106 4.2528 4.3973 4.5444 4.6939 Columns 28 through 36 4.8458 5.0002 5.1570 5.3163 5.4780 5.6422 5.8088 5.9779 6.1494 Columns 37 through 45 6.3234 6.4999 6.6787 6.8601 7.0439 7.2301 7.4188 7.6099 7.8035 Columns 46 through 51 7.9995 8.1980 8.3989 8.6023 8.8081 9.0164二， 编程解决以下科学计算问题1. 按图8-1-1所示的梯形直流电路，问（1） 如Us=10V，求Ubc，I7，Ude（2） 如Ude=4V，求Ubc，I7，Us解答：此电路中设节点电压为变量,共有ua,ub,uc,ud四个。将各支路电流均用这些电压来表示.对b点和c点列出节点电流方程：I1=I3+I7，I3=I4+I6，将上述各值代入化简，得对问题（1）给出Us时，这两个方程中只有两个未知数Ub和Uc，可以写成矩阵方程其中：a11=1/(r1+r2)+1/r3+1/r7;a12=-1/r3;a13=1/(r1+r2);a21=1/r3;a22=-1/r3-1/(r4+r5)-1/r6; a23=0并用矩阵左除解出Ub和Uc再由即可得到问题（1）的解。对问（2），可以推出类似的矩阵表达式，只是输入输出量不同开始流程图：输入Ude输入Us代入矩阵方程代入矩阵方程Ubc=u(1)-u(2)I7=u(1)/r7Us=(r1+r2)*(u(1)*a11+u(2)*a12)Ubc=u(1)-u(2)Ude=u(2)*r5/(r4+r5)I7=u(1)/r7源程序代码：r1=2;r2=4;r3=4;r4=4;r5=2;r6=12;r7=12;%为元件赋值%将各系数矩阵赋值a11=1/(r1+r2)+1/r3+1/r7;a12=-1/r3;a13=1/(r1+r2);a21=1/r3;a22=-1/r3-1/(r4+r5)-1/r6; a23=0; Us=input('Us='); % 输入解(1)的已知条件A1=a11,a12;a21,a22 % 列出系数矩阵A1u=A1a13*Us;0% u=ub;ucUbc=u(1)-u(2),Ude=u(2)*r5/(r4+r5),I7=u(1)/r7% 解(1)Ude=input('Ude='),% 输入解(2)的已知条件A2=A1,a13;a23;0,r5/(r4+r5),0 % 列出系数矩阵A2u=A20;0;Ude,%u=ub;uc;usUbc=u(1)-u(2),I7=u(1)/r7,Us=(r1+r2)*(u(1)*a11+u(2)*a12)% 解(2)运行结果：输入Us=10 时， Ubc = 2.2222, Ude = 0.7407,I7 = 0.3704输入Ude=4 时， Ubc = 12.0000, I7 = 2.0000, Us = 54.0000 如下：task2Us=10A1 = 0.5000 -0.2500 0.2500 -0.5000u = 4.4444 2.2222Ubc = 2.2222Ude = 0.7407I7 = 0.3704Ude=4Ude = 4A2 = 0.5000 -0.2500 0.1667 0.2500 -0.5000 0 0 0.3333 0u = 24.0000 12.0000 -54.0000Ubc = 12.0000I7 = 2.0000Us = 54.0000>>2. （1） 差商表：源代码：x=4.0 5.0 6.0 7.0 8.0; y=2.000 2.3607 2.44949 2.64575 2.82843; n=length(x); newton=x',y'for j=2:nfor i=n:-1:1if i>=jy(i)=(y(i)-y(i-1)./(x(i)-x(i-j+1);elsey(i)=0;endendnewton=newton,y'enddisp('下三角状的牛顿差商表如下：')newton运行结果:下三角状的牛顿差商表如下：newton = 4.0000 2.0000 0 0 0 0 5.0000 2.3607 0.3607 0 0 0 6.0000 2.4495 0.0888 -0.1360 0 0 7.0000 2.6458 0.1963 0.0537 0.0632 0 8.0000 2.8284 0.1827 -0.0068 -0.0202 -0.0209>>(2)写出牛顿多项式Nk（x），k=1,2,3,4.(3)在给定值x处求牛顿多项式Nk（x），k=1,2,3,4的值。(4)比较（3）中的结果与实际函数值。（2）,(3),(4)编为一个程序，（2）（3）以结果表示出来，（4）以图形表示出来流程图：源程序代码： %输入参数中x,y为元素个数相等的向量，t为待估计的点，可以为数字或向量。%输出参数中p2为所求得的牛顿插值多项式，z为利用多项式所得的t的函数值。x=4 5 6 7 8;y=2 2.3607 2.44949 2.64575 2.82843;t=4.5,7.5;n=length(x);chaS(1)=y(1);for i=2:n x1=x;y1=y; x1(i+1:n)=; y1(i+1:n)=; n1=length(x1); s1=0; for j=1:n1 t1=1; for k=1:n1 if k=j continue; else t1=t1*(x1(j)-x1(k); end end s1=s1+y1(j)/t1; end chaS(i)=s1;endb(1,:)=zeros(1,n-1) chaS(1);cl=cell(1,n-1);for i=2:n u1=1; for j=1:i-1 u1=conv(u1,1 -x(j); cli-1=u1; end cli-1=chaS(i)*cli-1; b(i,:)=zeros(1,n-i),cli-1;endp2=b(1,:);for q=2:n p2=p2+b(q,:); disp(strcat('牛顿',num2str(q-1),'次多项式')vpa(poly2sym(p2),6)%显示第（2）问牛顿多项式的结果，vpa函数用于精确到小数后4位if length(t)=1 rm=0; for i=1:n rm=rm+p2(i)*t(n-i); end z=rm;else k1=length(t); rm=zeros(1,k1); for j=1:k1 for i=1:n rm(j)=rm(j)+p2(i)*t(j)(n-i); end z=rm; endenddisp(strcat('牛顿',num2str(q-1),'次多项式的值')z %显示第（3）问给定x值处牛顿多项式的值subplot(2,2,(q-1)x0=0:0.01:10;plot(x0,sqrt(x0),'-',t,z,'*')title(strcat('牛顿',num2str(q-1),'次多项式与实际函数值比较')%作图比较（3）中结果与实际函数值end运行结果：牛顿1次多项式 ans = .*x+. 牛顿1次多项式的值z = 2.1803 3.2625牛顿2次多项式 ans = -.*x2+1.58430*x-2.16190 牛顿2次多项式的值z = 2.2143 2.0728牛顿3次多项式 ans = .e-1*x3-1.08440*x2+6.26332*x-9.74950 牛顿3次多项式的值z = 2.2380 2.9027牛顿4次多项式 ans = -.e-1*x4+.*x3-4.81678*x2+19.5664*x-27.2646 牛顿4次多项式的值z = 2.2576 2.7659>>
