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

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-范

3、数:');disp(Norm(y1,2);disp('y的无穷-范数:');disp(Norm(y1,inf);disp('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(&

4、#39;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(Norm(y3,1);disp('y的2-范数:');disp(Norm(y3,2);disp('y的无穷-范数:');disp(Norm(y

5、3,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-范数: 500500； y的2-范数:1.8271e+004；y的无穷-范数:1000程序Test2.mclear all;clc;n=100;

6、%区间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); if(d=1) f2(k)=log(d)/(d-1); else f2(k)=

7、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):test3.mclear all;clc;x=

8、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(

9、'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分解，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

10、; for k=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分解，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); e

11、nd end a=i;b=abs(t(i); for j=i+1:m if b<abs(t(j) b=abs(t(j); a=j; end end if a=i for j=1:m c=U(i,j); U(i,j)=U(a,j); U(a,j)=c; end for j=1:m c=P(i,j); P(i,j)=P(a,j); P(a,j)=c; end c=t(a); t(a)=t(i); t(i)=c; end U(i,i)=t(i); for j=i+1:m U(j,i)=t(j)/t(i); end for j=i+1:m for k=1:i-1 U(i,j)=U(i,j)-U(

12、i,k)*U(k,j); end endendL=tril(U,-1)+eye(m);U=triu(U,0);（1） （2）程序：Test4.mclear all;clc;for n=5:30 x=zeros(n,1); A=-ones(n); A(:,n)=ones(n,1); for i=1:n A(i,i)=1; for j=(i+1):(n-1) A(i,j)=0; end x(i)=1/i; end disp('当n=');disp(n); disp('方程精确解:'); x b=A*x; %系数b disp('利用LU分解方程组的解:'

13、;); L,U=LUDecom(A); %LU分解 xLU=U(Lb) disp('利用PLU分解方程组的解:'); P,L,U =PLUDecom(A); %PLU分解 xPLU=U(L(Pb) %求解A的逆矩阵 disp('A的准确逆矩阵:'); InvA=inv(A) InvAL=zeros(n); %利用LU分解求A的逆矩阵 I=eye(n); for i=1:n InvAL(:,i)=U(LI(:,i); end disp('利用LU分解的A的逆矩阵:'); InvALEnd运行结果：（1） 只列出n=5,6,7的结果当n= 5方程精确

14、解:x =1.0000 0.5000 0.3333 0.2500 0.2000利用LU分解方程组的解:xLU = 1.0000 0.5000 0.3333 0.2500 0.2000利用PLU分解方程组的解:xPLU = 1.0000 0.5000 0.3333 0.2500 0.2000当n=6方程精确解:x = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667利用LU分解方程组的解:xLU = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667利用PLU分解方程组的解:xPLU = 1.0000 0.5000 0.3333

15、0.2500 0.2000 0.1667当n= 7方程精确解:x = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429利用LU分解方程组的解:xLU = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429利用PLU分解方程组的解:xPLU = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429（2） 只列出n=5,6,7时A的逆矩阵的结果当n= 5A的准确逆矩阵:InvA = 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0

16、 0.5000 -0.2500 -0.1250 -0.1250 0 0 0.5000 -0.2500 -0.2500 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0625利用LU分解的A的逆矩阵:InvAL = 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0.5000 -0.2500 -0.2500 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0625当n= 6A的准确逆矩阵

17、:InvA = 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0313利用LU分解的A的逆矩阵:InvAL = 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0.5000

18、 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0313 0.0313当n= 7A的准确逆矩阵:InvA = 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0 0.5000 -0.25

19、00 -0.1250 -0.0625 -0.0625 0 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156利用LU分解的A的逆矩阵:InvAL = 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0 0.5

20、000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156程序：Cholesky分解：Cholesky.mfunction L=Cholesky(A)N = size(A);n = N(1);L=zeros(n);L(1,1)=sqrt(A(1,1);for i=2:n L(i,1)=A(i,1)/L(1

21、,1);endfor j=2:n s1=0; for k=1:j-1 s1=s1+L(j,k)2; end L(j,j)=sqrt(A(j,j)-s1); for i=j+1:n s2=0; for k=1:j-1 s2=s2+L(i,k)*L(j,k); end L(i,j)=(A(i,j)-s2)/L(j,j); endend计算Ax=b；Test5.mclear all;clc;for n=10:20 A=zeros(n,n); b=zeros(n,1); for i=1:n for j=1:n A(i,j)=1/(i+j-1); end b(i,1)=i; end disp('

22、n=');disp(n); disp('方程组原始解');x0=Ab disp('利用Cholesky分解的方程组的解'); L=Cholesky(A) x=L'(Lb)end运行结果：只列出了n=10,11的结果n=10方程组原始解x0 = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -1.2108 3.5947 -6.3233 6.5114 -3.6233 0.8407利用Cholesky分解的方程组的解x = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -1.210

23、5 3.5939 -6.3219 6.5100 -3.6225 0.8405n= 11方程组原始解x0 = 1.0e+009 * 0.0000 -0.0002 0.0046 -0.0567 0.3687 -1.4039 3.2863 -4.7869 4.2260 -2.0685 0.4305利用Cholesky分解的方程组的解x = 1.0e+009 * 0.0000 -0.0002 0.0046 -0.0563 0.3668 -1.3972 3.2716 -4.7669 4.2094 -2.06080.4290程序：（1）House.mfunction u=House(x)n=length(

24、x);e1=eye(n,1);w=x-norm(x,2)*e1;u=w/norm(w,2);(2)Hou_A.mfunction HA=Hou_A(A)a1=A(:,1);n=length(a1);e1=eye(n,1);w=a1-norm(a1,2)*e1;u=w/norm(w,2);H=eye(n)-2*u*u'HA=H*A;(3)test6.mclear all;clc;A=1 2 3 4; -1 2 sqrt(2) sqrt(3); -2 2 exp(1) pi; -sqrt(10) 2 -3 7; 0 2 7 5/2;HA=Hou_A(A)运行结果：H = 0.2500 -0

25、.2500 -0.5000 -0.7906 0 -0.2500 0.9167 -0.1667 -0.2635 0 -0.5000 -0.1667 0.6667 -0.5270 0 -0.7906 -0.2635 -0.5270 0.1667 0 0 0 0 0 1.0000HA = 4.0000 -2.5811 1.4090 -6.5378 0.0000 0.4730 0.8839 -1.7805 0.0000 -1.0541 1.6576 -3.8836 0.0000 -2.8289 -4.6770 -4.1078 0 2.0000 7.0000 2.5000程序：Jacobi迭代：Jacc

26、obi.mfunction x,n=Jaccobi(A,b,x0)%-·方程组系数阵A%-·方程组右端顶b%- 初始值x0%-求解要求精确度eps%-迭代步数控制M%-·返回求得的解x %-·返回迭代步数nM=1000;eps=1.0e-5;D=diag(diag(A); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵J=D(L+U);f=Db;x=J*x0+fn=1; %迭代次数err=norm(x-x0,inf)while(err>=eps) x0=x; x=J*x0+f n=n+1;

27、 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的对角矩阵L=-tril(A,-1)

28、; %求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=Abd

29、isp('迭代初始值');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 =

30、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.9

31、589err =0.5838.x = -0.0820 -1.8033 1.1311err = 9.4021e-006Gauss-Seidel最终迭代结果：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);dfvalu

32、e=subs(df,x0);for iter=1:maxiter x=x0-fvalue/dfvalue err=abs(x-x0) x0=x; fvalue=subs(f,x0) dfvalue=subs(df,x0); if(err<eps)|(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=s

33、ubs(f,x0);fvalue=subs(f,x1);for iter=1:maxiter x=x1-fvalue*(x1-x0)/(fvalue-fvalue0) err=abs(x-x1) x0=x1;x1=x; fvalue0=subs(f,x0);fvalue=subs(f,x1) if(err<eps)|(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迭代初始值&#

34、39;);xn1_0=0disp('Newton迭代结果');xn1,iter_n1,fxn1=Newtoniter(f,df,xn1_0,eps,maxiter)disp('Newton迭代初始值');xn2_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_

35、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)|x(k)f|x(k)-x(k-1)00-10051251101200103.809522.40581.190526.5714335.86013.42863.48371.39060.325834.546280.75692.02523.46070.00660.0

36、23043.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-0142.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.8004

37、2.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-004103.4606-2.8388e-0101.8302e-007

38、程序：二分法：resecm.mfunction x,iter=resecm(f,a,b,eps)%二分法 x 近似解%iter 迭代次数%f 求解的方程%a，b 求解区间%eps 允许误差限fa=subs(f,a);fb=subs(f,b);iter=0;if(fa=0) x=a;returnendif(fb=0) x=b;returnendwhile(abs(a-b)>=eps) mf=subs(f,(a+b)/2); if(mf=0) x=mf;n=n+1;return end if(mf*fa<0) b=(a+b)/2; else a=(a+b)/2; end iter=it

39、er+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=10e-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； 迭

40、代次数iter2 =202*pi,3*pi区间的根x3 =7.8548； 迭代次数iter3 =203*pi,4*pi区间的根x4 =10.9955；迭代次数iter4 =20程序：Newton插值：Newtominter.mfunction 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);

41、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); 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)的N

42、ewton插值：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(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'

43、;)%原图hold onplot(x0,f01,'-g')%5个节点plot(x0,f02,'-k')%10个节点plot(x0,f03,'-b')%15个节点legend('原图','5个节点Newton插值多项式','10个节点Newton插值多项式','15个节点Newton插值多项式')运行结果：取不同的节点做牛顿插值。得到结果图像如下：可以看出原图与插值多项式的图像近似重合，说明插值效果较好。程序：Lagrange插值：Lagrange.mfunction f,f0 = L

44、agrange(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(j)/(x(i)-x(j); end; for(j = i+1:n) l = l*(t-x(j)/(x(i)-x(j); %计算Lag

45、range基函数 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 all;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,'