1、大连理工大学矩阵与数值分析报告上机作业矩阵与数值分析上机作业 学校： 大连理工大学 学院： 班级： 姓名： 学号： 授课老师： 注：编程语言Matlab程序：Norm.m函数function s=Norm(x,m)%求向量x的范数%m取1,2,inf分别 表示1,2，无穷范数n=length(x);s=0;switch mcase 1 %1-范数for i=1:n s=s+abs(x(i);endcase 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
2、 all;clc;n1=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时);di
3、sp(x的1-范数:);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-
4、范数:);disp(Norm(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的无穷-范数:
5、1y的1-范数: 500500； 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=
6、1:n+1 d=1+x(k);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=
7、s;计算p(x):test3.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);运行结
8、果：x=2的邻域：x =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:nfor j=i:n z=0;for k=1:i-1 z=z+L(i
9、,k)*U(k,j);end U(i,j)=A(i,j)-z;endfor 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:mfor j=i:m t(j)=U(j,i);for k=1:i-1 t(j)=t(j)-U(j,k)*U(k,i);endend a=i;b=abs(t(i);for j=i+1
10、:mif 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); %求A的对
11、角矩阵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(迭代初始值);x0=0
12、;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.4533err =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.583
