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

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

这就是最小二乘法拟合。

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

注意n ≤m,用递推公式表示P k （x ）,即()()()()()()()01101111,,(1,2,,1)k k k k k P x P x x P x P x P x P x k n ααβ++-=⎧⎪=-⎨⎪=--=...-⎩ 这里的P k (x)是首项系数为1的k 次多项式,根据P k (x)的正交性,得()()()()()()()()()()()()()()()()()()()()2i 012i 02i 0211i 10x ,,x ,0,1,1,x ,0,1,1,x m i k i k k i k m k k k i i k k m k k k i k k i k m k k k i i x P x xP x P x a P x P x P x xP P k n P P P x P P k n P P P x ρρρβρ=+==---=⎧⎪⎪==⎪⎪⎪==⋅⋅⋅-⎨⎪⎪===⋅⋅⋅-⎪⎪⎪⎩∑∑∑∑ 根据公式(1)和(2)逐步求P k （x ）的同时,相应计算系数()()()020()(),(0,1,n (,)()m i j i k i k i k m k k ik i i x x x f P a k P P x x ρϕϕρϕ=====⋅⋅⋅,∑∑） 并逐步把*k a P k （x ）累加到S （x ）中去,最后就可得到所求的拟合函数曲线***0011n n y=S x =a P x +a P x ++a P x ⋅⋅⋅（）（）（）（）.流程图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))（2）（1）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Ω,C 1ω=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=x0<x1<…<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 x 2+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)在内节点x 1, x 2,....,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 个未知数,还差俩个方程。

通常在端点x 0=a. x n =b 处各附加一个条件称边界条件,常见有三种：（1）自然边界条件：S''(x 0)=S''(x n )=0（2）固定边界条件：S'(x 0)=f'(x 0),S'(x)=f'(x n )（3）周期边界条件：S'(x 0)=S'(x n ),S''(x 0)=S''(x n )共4n 个方程,可唯一的确定4n 个未知数流程图：（2） 运行流程图。