正交多项式最小二乘法拟合(共13页).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')二、 编程计算以下电路问题例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,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)=S''(xj+0) j=0,1,2,3.,n-1可得3n-2个方程，又由条件（3）S(xj)=yj j=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个未知数流程图：由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)=yj j=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) 开始 储存结果（1）用myspline的插值方法进行插值处理 开始 结束输入原始数据输出插值文件输出图形（2） 运行流程图 源代码：function s=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);for j=1:(n-1) h(j)=x0(j+1)-x0(j);%h(j)为第j个插值区间的长度endfor j=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);for j=(n-1):-1:1 m(j)=a(j)*m(j+1)+b(j);endfor k=1:km for j=1:(n-1) if (x(k)>=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(k)2-(x(k)-x0(l(k)3/h(l(k)3)*m(l(k)+1);end举例：x=-3 -2 1 4y=2 0 3 1x0=-3:0.15:4;y0=myspline(x,y,0,0,x0);plot(x0,y0)hold onplot(x,y,'or')legend('计算值','实验值')专心-专注-专业
