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

《大连理工大学矩阵与数值分析上机作业 .doc》由会员分享，可在线阅读，更多相关《大连理工大学矩阵与数值分析上机作业 .doc（55页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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-范数: 50

5、0500； 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

6、+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=s;计算p(

7、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);运行结果：x=2的

8、邻域：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=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分解，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=abs(t(i); for j=i+1:m

10、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); %求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(迭代初始值);

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.4533err

13、 =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最终迭代结果：xg =

14、 -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; fvalue=s

15、ubs(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-fvalue0) e

16、rr=abs(x-x1) x0=x1;x1=x; fvalue0=subs(f,x0);fvalue=subs(f,x1) if(erreps)|(fvalue=0),break,endend求方程的实根：test8.mclear all；clc;syms xf=x.3+2*x.2+10*x-100;df=diff(f,x,1);eps=10e-6;maxiter=100;disp(Newton迭代初始值);xn1_0=0disp(Newton迭代结果);xn1,iter_n1,fxn1=Newtoniter(f,df,xn1_0,eps,maxiter)disp(Newton迭代初始值);x

17、n2_0=5disp(Newton迭代结果);xn2,iter_n2,fxn2=Newtoniter(f,df,xn2_0,eps,maxiter)disp(弦截法初始值);xk1_0=0xk1_1=1disp(弦截法迭代结果);xk1,iter_k1,fxk1=secant(f,xk1_0,xk1_1,eps,maxiter)disp(弦截法初始值);xk2_0=5 xk2_1=6disp(弦截法迭代结果);xk2,iter_k2,fxk2=secant(f,xk2_0,xk2_1,eps,maxiter) 运行结果：Newton法结果:取两个不同初值0，5kx(k)f|x(k)-x(k-1

18、)|x(k)f|x(k)-x(k-1)00-10051251101200103.809522.40581.190526.5714335.86013.42863.48371.39060.325834.546280.75692.02523.46070.00660.023043.650811.82090.89543.46061.1043e-0041.5098e-00753.46770.42770.18303.4606-2.8422e-0142.5261e-00963.46066.3111e-0040.007173.46061.3805e-0091.0559e-00583.4606-2.8422e-0

19、142.3098e-011弦截法迭代结果：取两种不同初值0,1；5,6kx(k)f|x(k)-x(k-1)|x(k)f|x(k)-x(k-1)00-100512511-87162481.190527.6923550.43246.69233.983734.80042.016331.9134-66.5387 5.77893.654612.07110.329142.5366-45.44240.62323.47981.15540.174853.879127.25841.34253.46130.04500.018563.3758-4.98050.50343.46061.7917e-0047.5046e-00473.4535-0.42060.07783.46062.7963e-0082.9972e-00683.45350.00750.007293.4606-1.0939e-0051.2544e-00410