《MATLAB 程序设计实践》课程考核一、编程实现以下科学计算算法，并举一例应用之。

（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“正交多项式最小二乘法拟合”正交多项式最小二乘法拟合原理正交多项式做最小二乘法拟合：不要求拟合函数y=f(x)经过所有点（x i,y i)，而只要求在给定点x i上残差δi=f(x i)-y i按照某种标准达到最小，通常采用欧式范数||δ||2作为衡量标准。

这就是最小二乘法拟合。

根据作为给定节点x0，x1，…x m及权函数ρ(x)>0,造出带权函数正交的多项式{P n（x）}。

注意n≤m,用递推公式表示P k（x），即()()()()()()()1101111,,(1,2,,1) k k k k kP xP x x P xP x P x P x k nααβ++-=⎧⎪=-⎨⎪=--=...-⎩这里的P k(x)是首项系数为1的k次多项式，根据P k(x)的正交性，得()()()()()()()()()()()()()()()()()()()()2i12i2i211i1x,,x,0,1,1,x,0,1,1,xmi k ik kik mk kk ii k kmk kk ik kik mk kk iix P xxP x P xaP x P xP xxP Pk nP PP xP Pk nP PP xρρρβρ=+==---=⎧⎪⎪==⎪⎪⎪==⋅⋅⋅-⎨⎪⎪===⋅⋅⋅-⎪⎪⎪⎩∑∑∑∑根据公式(1)和(2)逐步求P k（x）的同时，相应计算系数()()()2()(),(0,1,n(,)()mi j i k ik ik mk ki k iix x xf Pa kP Px xρϕϕρϕ=====⋅⋅⋅,∑∑）并逐步把*ka P k（x）累加到S（x）中去，最后就可得到所求的拟合函数曲线***0011n ny=S x=a P x+a P x++a P x⋅⋅⋅（）（）（）（）.流程图（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',...'X and Y vectors must be the same size.') end%检验X Y维数是否匹配x = x(:);y = y(:);if nargout > 2mu = [mean(x); std(x)];x = (x - mu(1))/mu(2);end%利用范德蒙德矩阵构造方程组系数矩阵V(:,n+1) = ones(length(x),1,class(x));for j = n:-1:1V(:,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.0e10if nargout > 2warning('MATLAB:polyfit:RepeatedPoints', ...'Polynomial is badly conditioned. Remove repeated data points.')elsewarning('MATLAB:polyfit:RepeatedPointsOrRescale', ...['Polynomial is badly conditioned. Remove repeated data points\n'...' or try centering and scaling as described in HELP POLYFIT.'])endendr = y - V*p;p = p.'; % 将多项式系数默认为行向量.5、运行流程图过程：clearx =[ 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000] y=[1.75 2.45 3.81 4.80 8.00 8.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,'*')hold onplot(x1,y1,'r')二、编程计算以下电路问题[例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('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），而且自由边界条件S ''(-3)=0,S ''(4)=0。

算法说明：三次样条也是分片三次插值函数，它是在给定的区间[a ，b]上的一个划分：a=x 0<x 1<…<x n =b,已知函数f(x)在xj 上的函数值为 f(x j )=y j ，（j=0,1,2,3...,n)如果存在分段函数[][][]⎪⎪⎩⎪⎪⎨⎧∈⋯⋯∈∈=-n n n x x x x S x x x x S x x x x S x S ,)(,)(,)()(1212101满足下述条件：（1）S(x)在每一个子区间 [ x j-1，x j ]，（j=0,1,2,3...,n)上是三次多项式； （2）S(x)在每一个内接点 （j=0,1,2,3...,n)具有直到二阶的连续导数； 则成为节点x 0，x 1，…，x n 上的三次样条函数。

若 S(x) 在节点 x 0，x 1，…，x n 上海满足插值条件：（3）S(x j)=y j （j=0,1,2,3...,n)则称S(x)为三次样条插值函数。

由（1）知，S(x)在每一个小区间[ x j-1，x j ]上是一三次多项式，若记为S j(x)，则可设：S j(x)=a j x3+b j x2+c j x+d j要确定函数S(x)的表达式，需确定4n个未知数{aj,bj ,cj, dj}（j=0,1,2,3...,n）。

由（2）知S(x)，S'(x)，S''(x)在内节点x1, x2,....,x n-1上连续，则Ⅰ型S(x j-0)=S(x j+0)Ⅱ型S'(x j-0)=S'(x j+0)Ⅲ型S''(x j-0)=S''(x j+0) j=0,1,2,3...,n-1可得3n-2个方程，又由条件（3）S(x j)=y j j=0,1,2,3...,n得n个方程，共得到4n-2个方程。

要确定4n个未知数，还差俩个方程。

通常在端点x0=a. x n=b处各附加一个条件称边界条件，常见有三种：（1）自然边界条件：S''(x0)=S''(x n)=0（2）固定边界条件：S'(x0)=f'(x0),S'(x)=f'(x n)（3）周期边界条件：S'(x0)=S'(x n)，S''(x0)=S''(x n)共4n个方程，可唯一的确定4n个未知数流程图：（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:1m(j)=a(j)*m(j+1)+b(j);endfor k=1:kmfor j=1:(n-1)if ((x(k)>=x0(j))&(x(k)<x0(j+1)))l(k)=j;endendendfor k=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 -2 1 4]y=[2 0 3 1]x0=[-3:0.15:4];y0=myspline(x,y,0,0,x0);plot(x0,y0)hold onplot(x,y,'or')legend('计算值','实验值')。