大连理工大学矩阵与数值分析上机作业

1、矩阵与数值分析上机作业 学校： 大连理工大学 学院： 班级： 姓名： 学号： 授课老师： 注：编程语言Matlab程序：Norm.m函数function s=Norm(x,m)%求向量x的范数%m取1,2,inf分别 表示1,2，无穷范数n=length(x);s=0;switch m case 1 %1-范数 for i=1:n s=s+abs(x(i); end case 2 %2-范数 for i=1:n s=s+x(i)2; end s=sqrt(s); case inf %无穷-范数 s=max(abs(x);end 计算向量x，y的范数Test1.mclear all;clc;n1

2、=10;n2=100;n3=1000;x1=1./1:n1;x2=1./1:n2;x3=1./1:n3;y1=1:n1;y2=1:n2;y3=1:n3;disp(n=10时);disp(x的1-范数:);disp(Norm(x1,1);disp(x的2-范数:);disp(Norm(x1,2);disp(x的无穷-范数:);disp(Norm(x1,inf);disp(y的1-范数:);disp(Norm(y1,1);disp(y的2-范数:);disp(Norm(y1,2);disp(y的无穷-范数:);disp(Norm(y1,inf);disp(n=100时);disp(x的1-范数:)

3、;disp(Norm(x2,1);disp(x的2-范数:);disp(Norm(x2,2);disp(x的无穷-范数:);disp(Norm(x2,inf);disp(y的1-范数:);disp(Norm(y2,1);disp(y的2-范数:);disp(Norm(y2,2);disp(y的无穷-范数:);disp(Norm(y2,inf);disp(n=1000时);disp(x的1-范数:);disp(Norm(x3,1);disp(x的2-范数:);disp(Norm(x3,2);disp(x的无穷-范数:);disp(Norm(x3,inf);disp(y的1-范数:);disp(N

4、orm(y3,1);disp(y的2-范数:);disp(Norm(y3,2);disp(y的无穷-范数:);disp(Norm(y3,inf);运行结果：n=10时x的1-范数:2.9290；x的2-范数:1.2449； x的无穷-范数:1y的1-范数:55； y的2-范数:19.6214； y的无穷-范数:10n=100时x的1-范数:5.1874；x的2-范数: 1.2787； x的无穷-范数:1y的1-范数:5050； y的2-范数:581.6786； y的无穷-范数:100n=1000时x的1-范数:7.4855； x的2-范数:1.2822； x的无穷-范数:1y的1-范数: ；

5、y的2-范数:1.8271e+004；y的无穷-范数:1000程序Test2.mclear all;clc;n=100;%区间h=2*10(-15)/n;%步长x=-10(-15):h:10(-15);%第一种原函数f1=zeros(1,n+1);for k=1:n+1 if x(k)=0 f1(k)=log(1+x(k)/x(k); else f1(k)=1; endendsubplot(2,1,1);plot(x,f1,-r);axis(-10(-15),10(-15),-1,2);legend(原图);%第二种算法f2=zeros(1,n+1);for k=1:n+1 d=1+x(k);

6、 if(d=1) f2(k)=log(d)/(d-1); else f2(k)=1; endendsubplot(2,1,2);plot(x,f2,-r);axis(-10(-15),10(-15),-1,2);legend(第二种算法);运行结果：显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当时，计算机进行舍入造成恒等于1，结果函数值恒为1。程序：秦九韶算法:QinJS.mfunction y=QinJS(a,x)%y输出函数值%a多项式系数，由高次到零次%x给定点n=length(a);s=a(1);for i=2:n s=s*x+a(i);endy=s;计算p(x):tes

7、t3.mclear all;clc;x=1.6:0.2:2.4;%x=2的邻域disp(x=2的邻域：);xa=1 -18 144 -672 2016 -4032 5376 -4608 2304 -512;p=zeros(1,5);for i=1:5 p(i)=QinJS(a,x(i);enddisp(相应多项式p值：);pxk=1.95:0.01:20.5;nk=length(xk);pk=zeros(1,nk);k=1;for k=1:nk pk(k)=QinJS(a,xk(k);endplot(xk,pk,-r);xlabel(x);ylabel(p(x);运行结果：x=2的邻域：x =

8、1.6000 1.8000 2.0000 2.2000 2.4000相应多项式p值：p = 1.0e-003 * -0.2621 -0.0005 0 0.0005 0.2621p(x)在1.95,20.5上的图像程序：LU分解，LUDecom.mfunction L,U=LUDecom(A)%不带列主元的LU分解N = size(A);n = N(1);L=eye(n);U=zeros(n);for i=1:n U(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1);endfor i=2:n for j=i:n z=0; for k=1:i-1 z=z+L(i,k)*U(k,j)

9、; end U(i,j)=A(i,j)-z; end for j=i+1:n z=0; for k=1:i-1 z=z+L(j,k)*U(k,i); end L(j,i)=(A(j,i)-z)/U(i,i); endendPLU分解，PLUDecom.mfunction P,L,U =PLUDecom(A)%带列主元的LU分解m,m=size(A);U=A;P=eye(m);L=eye(m);for i=1:m for j=i:m t(j)=U(j,i); for k=1:i-1 t(j)=t(j)-U(j,k)*U(k,i); end end a=i;b=abs(t(i); for j=i+

10、1:m if b=eps) x0=x; x=J*x0+f n=n+1; err=norm(x-x0,inf) if(n=M) disp(Warning: 迭代次数太多，可能不收敛?); return; endendGauss_Seidel迭代：Gauss_Seidel.mfunction x,n=Gauss_Seidel(A,b,x0)%-Gauss-Seidel迭代法解线性方程组%-方程组系数阵 A%-方程组右端项 b%-初始值 x0%-求解要求的精确度 eps%-迭代步数控制 M%-返回求得的解 x %-返回迭代步数 neps=1.0e-5;M=10000;D=diag(diag(A);

11、%求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵G=(D-L)U;f=(D-L)b;x=G*x0+fn=1; %迭代次数err=norm(x-x0,inf)while(err=eps) x0=x; x=G*x0+f n=n+1; err=norm(x-x0,inf) if(n=M) disp(Warning: 迭代次数太多，可能不收敛！); return; endend解方程组，test7.mclear all;clc;A=5 -1 -3; -1 2 4; -3 4 15;b=-2;1;10;disp(精确解);x=Abdisp(迭代初

12、始值);x0=0;0;0disp(Jacobi迭代过程：);xj,nj=Jaccobi(A,b,x0);disp(Jacobi最终迭代结果：);xjdisp(迭代次数);njdisp(Gauss-Seidel迭代过程：);xg,ng=Gauss_Seidel(A,b,x0);disp(Gauss-Seidel最终迭代结果：);xgdisp(迭代次数);ng运行结果：精确解x = -0.0820 -1.8033 1.1311迭代初始值x0 = 0 0 0Jacobi迭代过程：x = -0.4000 0.5000 0.6667err = 0.6667x = 0.1000 -1.0333 0.453

13、3err =1.5333.x = -0.0820 -1.8033 1.1311err = 9.6603e-006Jacobi最终迭代结果：xj = -0.0820 -1.8033 1.1311迭代次数nj = 281Gauss-Seidel迭代过程：x = -0.4000 0.3000 0.5067err = 0.5067x = -0.0360 -0.5313 0.8012err = 0.8313x = -0.0256 -1.1151 0.9589err =0.5838.x = -0.0820 -1.8033 1.1311err = 9.4021e-006Gauss-Seidel最终迭代结果：

14、xg = -0.0820 -1.8033 1.1311迭代次数ng =20程序：Newton迭代法：Newtoniter.mfunction x,iter,fvalue=Newtoniter(f,df,x0,eps,maxiter)%牛顿法 x得到的近似解%iter迭代次数%fvalue函数在x处的值%f，df被求的非线性方程及导函数%x0初始值%eps 允许误差限%maxiter 最大迭代次数fvalue=subs(f,x0);dfvalue=subs(df,x0);for iter=1:maxiter x=x0-fvalue/dfvalue err=abs(x-x0) x0=x; fval

15、ue=subs(f,x0) dfvalue=subs(df,x0); if(erreps)|(fvalue=0),break,endend弦截法：secant.mfunction x,iter,fvalue=secant(f,x0,x1,eps,maxiter)%弦截法 x得到的近似解%iter迭代次数%fvalue函数在x处的值%f被求的非线性方程%x0，x1初始值%eps 允许误差限%maxiter 最大迭代次数fvalue0=subs(f,x0);fvalue=subs(f,x1);for iter=1:maxiter x=x1-fvalue*(x1-x0)/(fvalue-fvalue

16、0) err=abs(x-x1) x0=x1;x1=x; fvalue0=subs(f,x0);fvalue=subs(f,x1) if(err=eps) mf=subs(f,(a+b)/2); if(mf=0) x=mf;n=n+1;return end if(mf*fa0) b=(a+b)/2; else a=(a+b)/2; end iter=iter+1;endx=(a+b)/2;iter=iter+1;求方程的实根：test9.mclear all;clc;syms xf=exp(x).*cos(x)+2;a=0;a1=pi;a2=2*pi;a3=3*pi;b=4*pi;eps=10

17、e-6;x1,iter1=resecm(f,a,a1,eps)x2,iter2=resecm(f,a1,a2,eps)x3,iter3=resecm(f,a2,a3,eps)x4,iter4=resecm(f,a3,b,eps) 运行结果：0,pi区间的根x1 =1.8807； 迭代次数iter1 = 20pi,2*pi区间的根x2 =4.6941； 迭代次数iter2 =202*pi,3*pi区间的根x3 =7.8548； 迭代次数iter3 =203*pi,4*pi区间的根x4 =10.9955；迭代次数iter4 =20程序：Newton插值：Newtominter.mfunction

18、f=Newtominter(x,y,x0)%牛顿插值 x插值节点%y为对应的函数值%函数返回Newton插值多项式在x_0点的值fsyms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(x和y的维数不相等！); return;end f = y(1);y1 = 0;l = 1;for(i=1:n-1) for(j=i+1:n) y1(j) = (y(j)-y(i)/(x(j)-x(i); end c(i) = y1(i+1); l = l*(t-x(i); f = f + c(i)*l; simplify(f)

19、; y = y1; if(i=n-1) if(nargin = 3) %如果3个参数则给出插值点的插值结果 f = subs(f,t,x0); else %如果2个参数则直接给出插值多项式 f = collect(f); %将插值多项式展开 f = vpa(f, 6); end endend用等距节点做f(x)的Newton插值：test10.mn1=5;n2=10;n3=15;x0=0:0.01:1;y0=sin(pi.*x0);x1=linspace(0,1,n1);%等距节点,节点数5y1=sin(pi.*x1);f01=Newtominter(x1,y1,x0);x2=linspace

20、(0,1,n2);%等距节点,节点数10y2=sin(pi.*x2);f02 = Newtominter(x2,y2,x0);x3=linspace(0,1,n3);%等距节点,节点数15y3=sin(pi.*x3);f03= Newtominter(x3,y3,x0);plot(x0,y0,-r)%原图hold onplot(x0,f01,-g)%5个节点plot(x0,f02,-k)%10个节点plot(x0,f03,-b)%15个节点legend(原图,5个节点Newton插值多项式,10个节点Newton插值多项式,15个节点Newton插值多项式)运行结果：取不同的节点做牛顿插值。得

21、到结果图像如下：可以看出原图与插值多项式的图像近似重合，说明插值效果较好。程序：Lagrange插值：Lagrange.mfunction f,f0 = Lagrange(x,y,x0) %Lagrange插值 x为插值结点，y为对应的函数值，x0为要计算的点。%函数返回L_n(x)表达式f和L_n(x0)的值f0。syms t;if(length(x) = length(y) n = length(x); else disp(x和y的维数不相等！); return;end %检错 f = 0.0;for(i = 1:n) l = y(i); for(j = 1:i-1) l = l*(t-x

22、(j)/(x(i)-x(j); end; for(j = i+1:n) l = l*(t-x(j)/(x(i)-x(j); %计算Lagrange基函数 end; f = f + l; %计算Lagrange插值函数 simplify(f); %化简 if(i=n) if(nargin = 3) f0 = subs(f,t,x0); %如果3个参数则计算插值点的函数值 else f = collect(f); %如果2个参数则将插值多项式展开 f = vpa(f,6); %将插值多项式的系数化成6位精度的小数 end endend用等距节点做Lagrange插值：test11.mclear a

23、ll;clc;n1=5;n2=10;n3=15;x0=-5:0.02:5;y0=1./(1+x0.2);x1=linspace(-5,5,n1);%等距节点,节点数5y1=1./(1+x1.2);f1,f01 = Lagrange(x1,y1,x0);x2=linspace(-5,5,n2);%等距节点,节点数10y2=1./(1+x2.2);f2,f02 = Lagrange(x2,y2,x0);x3=linspace(-5,5,n3);%等距节点,节点数15y3=1./(1+x3.2);f3,f03 = Lagrange(x3,y3,x0);plot(x0,y0,-r)%原图hold on

24、plot(x0,f01,-b)%5个节点plot(x0,f02,-g)%10个节点plot(x0,f03,-k)%15个节点xlabel(x);ylabel(f(x),L(x);legend(原图,5个节点Lagrange插值多项式,10个节点Lagrange插值多项式,15个节点Lagrange插值多项式)运行结果：选取了5,10,15个节点做Lagrange插值，得到原图与插值多项式的图像如下：当选取节点数较多时，选取15个节点时出现Runge现象。程序：复合梯形求积：Comtrap.mfunction s=Comtrap(f,a,b,n)%复合梯形求积 s积分结果%f积分函数%a，b积分区间%n 区间个数h=(b-a)/n;index=(a+h):h:(b-h);s1=sum(subs(f,index);s=(h/2)*(subs(f,a)+2*s1+subs(f,b);复合Simpson求积：function s=Simpson(f,a,b,n)%复合Simpson求积 s积分结果%f积分函数%a，b积分区间%n 区间个数h=(b-a)/(2*n);index1=(a+h):(2*h):(b-h);index2=(a+2*h):(2*h):(b-2*h);s1=sum(subs