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

1、大连理工大学矩阵与数值分析上机作业18478矩阵与数值分析上机作业 学校： 大连理工大学 学院： 班级： 姓名： 学号： 授课老师： 注：编程语言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的范数T

2、est1.mclear 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

3、(n=100时);disp(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

4、);disp(y的1-范数:);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.282

5、2； x的无穷-范数: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(

6、1,n+1);for k=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

7、s=s*x+a(i);endy=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)

8、;ylabel(p(x);运行结果：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分解，LUDe.mfunction L,U=LUDe.(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=

9、1:i-1 z=z+L(i,k)*U(k,j); 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分解，PLUDe.mfunction P,L,U =PLUDe.(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=

10、abs(t(i); for j=i+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=1000

11、0;D=diag(diag(A); %求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

12、(精确解);x=Abdisp(迭代初始值);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.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