数值分析考试复习题 附录：matlab函数 .doc

模式识别_周晓勇_Z201102021附录：函数%function p=newpoly(x,y,n,x0)%p=newpoly(x,y,n)%(x,y)为节点坐标%p为所求的的差值多项式的系数向量%n为牛顿插值次数,n=1时为线性插值,n=2时为二次插值.N=length(x);d=x;y;if nargin=4 for i=1:N %根据插值误差估计式(根据距x0的远近)重新排列插值节点 for j=i+1:N if abs(d(1,i)-x0)abs(d(1,j)-x0) t=d(:,i); d(:,i)=d(:,j); d(:,j)=t; end end endenddq=zeros(n+1,n+2);dq(:,1)=d(1,1:n+1);dq(:,2)=d(2,1:n+1);for i=2:n+1 for j=3:i+1 dq(i,j)=(dq(i,j-1)-dq(i-1,j-1)/(dq(i,1)-dq(i-j+2,1);%差商表 endends=diag(dq(2:end,3:end);p=dq(1,2);t=1;for i=1:n t=conv(t,1 -dq(i,1); p=polysum(p,s(i).*t);%函数见polysum.mend%function x,k,flag=SOR(A,b,delta,w,step)%函数格式：x,k,flag=SOR(A,b,delta,w,step)%A为方程组的系数矩阵%b为方程组右端项%delta为精度要求%step为最大迭代次数,缺省值为100%w为松弛因子,w1时为超松弛法,w=1时为高斯-塞德尔迭代%x为方程组的解%k为迭代次数%flag为收敛标志if nargin5 step=100;endif nargin4 w=1;endif nargin3 ep=1e-5;endn=length(A);k=0;x=zeros(n,1);flag=1;while 1 y=x; for i=1:n z=b(i); for j=1:n if j=i z=z-A(i,j)*x(j); end end if abs(A(i,i)1e-10|k=step flag=0; return; end z=z/A(i,i); x(i)=(1-w)*x(i)+w*z; end if norm(y-x,inf)delta break; end k=k+1;end%function y=GaussLegendre(func,a,b,n)%y=GaussLegendre(func,a,b,n)%func:被积函数表达式,例如x.2.*cos(x)%a:积分下限%b:积分上限%n:节点数xTable=. 0.0 NaN 0 0 0 0; 0.5773503 -0.5773503 0 0 0 0; 0.7745967 -0.7745967 0.0 0 0 0; 0.8611363 -0.8611363 0.3399810 -0.3399810 0 0; 0.9061798 -0.9061798 0.5384693 -0.5384693 0.0 0; 0.9324685 -0.9324685 0.6612094 -0.6612094 0.2386192 -0.2386192;ATable=. 2.0 0 0 0 0 0; 1.0 1.0 0 0 0 0; 0.5555556 0.5555556 0.8888889 0 0 0; 0.3478548 0.3478548 0.6521452 0.6521452 0 0; 0.2369269 0.2369269 0.4786287 0.4786287 0.5688889 0; 0.1713245 0.1713245 0.3607616 0.3607616 0.4679139 0.4679139;f=inline(func);x=xTable(n,:);A=ATable(n,:);T=zeros(1,n);T=(a+b)/2+(b-a)/2*x;y=(b-a)/2*sum(A.*feval(f,T);%function y,Ck,Ak=NewtonCotes(func,a,b,n)% y,Ck,Ak=NewtonCotes(fun,a,b,n)% func:积分表达式，这里有两种选择% (1)积分函数句柄，必须能够接受矢量输入，比如func=(x)sin(x).*cos(x)% (2)x,y坐标的离散点，第一列为x，第二列为y，必须等距，且节点的个数小于9，% 比如：func=1:8;sin(1:8)如果func的表采用第二种方式，% 那么只需要输入第一个参数即可，否则还要输入a,b,n三个参数% a:积分下限% b:积分上限% n:牛顿-科特斯数公式的阶数，必须满足1n=8时不能保证公式的稳定性% (1)n=1，即梯形公式% (2)n=2，即辛普森公式% (3)n=4，即科特斯公式% y:数值积分结果% Ck:科特斯系数% Ak:求积系数if nargin=1 mm,nn=size(func); if mm=8 error(为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！) elseif nn=2 error(func应为：第一列为x，第二列为y，并且个数为小于10的等距节点！) end xk=func(:,1); fk=func(:,2); a=min(xk); b=max(xk); n=mm-1;elseif nargin=4 % 计算积分节点xk和节点函数值fx xk=linspace(a,b,n+1); if isa(func,function_handle)|isa(func,inline) fk=func(xk); %当func=sin(x)/x且x=0时出错,fx(1)=NaN if sym(func)=sym(x)(sin(x)./x) for i=1:length(xk) if isnan(fk(i) fk(i)=1;end end end %当func=sin(x)/x且x=0时出错,定义sin(0)/0=1 else error(输入函数应为匿名函数或内联函数！) endelse error(输入参数错误，请参考函数帮助！)end% 计算科特斯系数Ck=cotescoeff(n);%函数见cotescoeff.m% 计算求积系数Ak=(b-a)*Ck;% 求和算积分y=Ak*fk;%function x,err,k,y=NewtonDownhill(func,x0,delta,step)%x,err,k,y=newton(func,x0,delta,step)%func是非线性函数表达式,例如1-exp(-x2)%x0是初值%delta是给定的允许误差%step是给定的允许迭代最大次数%x是牛顿法求的的方程近似解%err是误差估计%k是实际迭代次数%y=f(x)f=inline(func);df=inline(diff(sym(func);C=1;%绝对误差或相对误差控制常数,一般取C=1for k=1:step lambda=1;%牛顿下山因子 x(k)=x0-lambda*f(x0)/df(x0); while abs(f(x(k)=abs(f(x0) lambda=lambda/2; x(k)=x0-lambda*f(x0)/df(x0); end if x(k)C err(k)=abs(x(k)-x0); else err(k)=abs(x(k)-x0)/abs(x(k); end x0=x(k); y(k)=f(x(k); if(err(k)delta)|(y(k)=0) break; endendfigure;fplot(f,-2.5 2.5);grid;hold onplot(x(k),y(k),r*);h=legend(f(x)=,func,根x=,num2str(x(k),迭代次数k=,num2str(k);set(h,Location,SouthEast);title(牛顿下山法解非线性方程);%function y=mulNewtonCotes(func,a,b,n,m,delta)% y=mulNewtonCotes(func,a,b,n,m)或% y=mulNewtonCotes(func,a,b,n,Nan,delta)% 复化Newton-Cotes数值积分公式，即在每个子区间上使用Newton-Cotes公式，然后求和% func:被积函数表达式例如exp(x)% a,b:积分下上限% m:将区间a,b等分的子区间数量% delta:给定误差,当有m输入时不输入delta,否则根据delta重新计算m% n:采用的Newton-Cotes公式的阶数，必须满足n8，否则积分没法保证稳定性% (1)n=1，即复化梯形公式% (2)n=2，即复化辛普森公式% (3)n=4，即复化科特斯公式f=inline(func);if nargin=6 switch n case 1 d2f=inline(diff(sym(func),2); m=abs(b-a)3*d2f(max(a,b)/60/delta; m=ceil(sqrt(m); case 2 d4f=inline(diff(sym(func),4); m=abs(b-a)5*d4f(max(a,b)/14400/delta; m=ceil(sqrt(sqrt(m); endendxk=linspace(a,b,m+1);for i=1:m s(i)=NewtonCotes(f,xk(i),xk(i+1),n);%函数见NewtonCotes.mendy=sum(s);%function x,k,flag=Jacobi(A,b,delta,step)%函数格式：x,k,flag=jacobi(A,b,delta,step)%A为方程组的系数矩阵%b为方程组右端项%delta为精度要求%step为最大迭代次数,缺省值为100%x为方程组的解%k为迭代次数%flag为收敛标志if nargin4 step=100;end;if nargin3 step=1e-5;end;D=diag(diag(A);L=tril(A,-1);U=triu(A,1);M=D;N=M-A;B=MN;max(abs(eig(B);f=Mb;flag=1;%1代表迭代格式收敛k=1;x0=0 0 0;x=B*x0+f;while 1 x0=x; x=B*x0+f; if norm(x-x0,inf)delta break; end if k=step flag=0;%0代表迭代格式在规定的迭代次数下不收敛 break; end k=k+1;end%function y,T,err,h=Romberg(func,a,b,n,delta)%y,T,err,h=Romberg(func,a,b,n,delta)%func是被积函数表达式,例如x./(4+x.2)%a,b是积分下限和上限%n+1是T数表的列数%delta是允许的误差,缺省值为1e-5%T是T数表%y是所求积分值%h是子区间长度if nargindelta)&(Jn)|(J=abs(f(x1) lambda=lambda/2; x(k)=x1-lambda*f(x1)*(x1-x0)/(f(x1)-f(x0); end if x(k)C err(k)=abs(x(k)-x1); else err(k)=abs(x(k)-x1)/abs(x(k); end x0=x1; x1=x(k); y(k)=f(x(k); if(err(k)delta)|(y(k)=0) break; endendfigure;fplot(f,-2.5 2.5);grid;hold onplot(x(k),y(k),r*);h=legend(f(x)=,func,根x=,num2str(x(k),迭代次数k=,num2str(k);set(h,Location,SouthEast);title(弦截法(含下山因子)解非线性方程);%function q=rectint(x,y,a,b,type)%q=rectint(x,y,a,b,type)%(x,y)离散积分节点%a,b:积分上下限%type:矩形公式类型,left为左矩形,right为右矩形,middle为中矩形%若不输入type,则缺省值为middle%q:积分结果N=length(x);t=zeros(1,N-1);if narginlength(p2) p2=zeros(1,length(p1)-length(p2),p2; p_out=p1+p2;else p1=zeros(1,length(p2)-length(p1),p1; p_out=p1+p2;end%function f=intfun(t,n,k)% 科特斯系数中的积分表达式f=1;for j=0:k-1,k+1:n f=f.*(t-j);end%function