大连理工大学矩阵与数值分析大作业.pdf

矩阵与数值分析上机实习矩阵与数值分析上机实习1.设,其精确值为.(1)编制按从大到小的顺序程序(2)编制按从小到大的顺序的通用程序(3)按两种顺序分别计算(4)通过本上机题，你明白了什么？从小到大,代码:%1-SN=%N=input(please input a number(N=2)if(N=2)102N=100S:7.4005e-001 clearplease input a number(N=2)104N=10000S:7.4990e-001 clearplease input a number(N=2)106N=1000000S:7.5000e-001从大到小代码:%1-SN=%eps(single)N=input(please input a number(N=2)if(N=2)102N=100S:7.4005e-001 clearans=1.1921e-007please input a number(N=2)104N=10000S:7.4990e-001 clearans=1.1921e-007please input a number(N=2)106N=1000000S:7.5000e-001(4)计算的顺序不同，结果有可能不同。2.秦九韶算法。已知 n 次多项式，用秦九韶算法编写通用的程序计算函数在点的值，并计算在点的值.(提示：编写程序时，输入系数向量和点，输出结果，多项式的次数可以通过向量的长度来判断)代码代码:A=input(请输入系数,由高次幂开始);n=input(请输入计算变量的值);len=length(A);val=zeros(len);val(1)=A(1);%printf(len=%c,len)for i=2:1:len%disp(val(i-1)%disp(n)val(i)=val(i-1)*n+val(i);end%printf(%f,val(len)val(len)结果结果:请输入系数，由高次幂开始7 3-5 11请输入计算变量的值23ans=851693.分别用 Gauss 消元法和列主元消去法编程求解方程组Ax=b，其中，.高斯消去法代码:A=31-13 0 0 0-10 0 0 0 -13 35-9 0-11 0 0 0 0 0-9 31-10 0 0 0 0 0 0 0-10 79-30 0 0 0-9 0 0 0-30 57-7 0-5 0 0 0 0 0-7 47-30 0 0 0 0 0 0 0-30 41 0 0 0 0 0 0-5 0 0 27-2 0 0 0-9 0 0 0-2 29 ;B=-15,27,23,0,-20,12,-7,7,10;B=B.;%disp(A)%disp(B)C=A,B;n=size(A,1);ra=rank(A);rc=rank(C);x=zeros(1,n);if ra=rc disp(no solution)elsefor i=1:1:(n-1)for j=i+1:1:n temp=(-C(j,i)/C(i,i);for k=i:1:(n+1)C(j,k)=C(i,k)*temp+C(j,k);endendendendfor i=n:-1:1if i=n disp(C(i,(n+1)disp(C(i,(n+1)-sum(C(i,(i+1):n).*x(i+1):n)disp(C(i,i)end x(i)=(C(i,(n+1)-sum(C(i,(i+1):n).*x(i+1):n)/C(i,i);enddisp(solution)disp(x)结果:-2.2051e-0029.4726e-0011.0410e+0007.4481e-002-2.6617e-0012.0020e-001-2.4242e-0022.3844e-0013.8439e-001列主元消去法代码:%t4-max%A=31-13 0 0 0-10 0 0 0 -13 35-9 0-11 0 0 0 0 0-9 31-10 0 0 0 0 0 0 0-10 79-30 0 0 0-9 0 0 0-30 57-7 0-5 0 0 0 0 0-7 47-30 0 0 0 0 0 0 0-30 41 0 0 0 0 0 0-5 0 0 27-2 0 0 0-9 0 0 0-2 29 ;B=-15,27,23,0,-20,12,-7,7,10;B=B.;%disp(A)%disp(B)C=A,B;n=size(A,1);ra=rank(A);rc=rank(C);x=zeros(1,n);if ra=rc disp(no solution)elsefor i=1:1:(n-1)%make max value in col in the row%max=0;inx=i;for j=i:nif(abs(C(j,i)max)max=abs(C(j,i);inx=j;endendif inx ifor j=i:(n+1)swap=C(i,j);C(i,j)=C(inx,j);C(inx,j)=swap;endendfor j=i+1:1:n temp=(-C(j,i)/C(i,i);for k=i:1:(n+1)C(j,k)=C(i,k)*temp+C(j,k);endendendendfor i=n:-1:1if i=n disp(C(i,(n+1)disp(C(i,(n+1)-sum(C(i,(i+1):n).*x(i+1):n)disp(C(i,i)end x(i)=(C(i,(n+1)-sum(C(i,(i+1):n).*x(i+1):n)/C(i,i);enddisp(solution)disp(x)结果:solution-2.2051e-0029.4726e-0011.0410e+0007.4481e-002-2.6617e-0012.0020e-001-2.4242e-0022.3844e-0013.8439e-0014.编程求解题 3 中矩阵的 LU 分解及列主元的 LU 分解（求出 L,U 和 P），并用 LU 分解的方法求 A 的逆矩阵及 A 的行列式LU 分解:代码:clearA=31-13 0 0 0-10 0 0 0 -13 35-9 0-11 0 0 0 0 0-9 31-10 0 0 0 0 0 0 0-10 79-30 0 0 0-9 0 0 0-30 57-7 0-5 0 0 0 0 0-7 47-30 0 0 0 0 0 0 0-30 41 0 0 0 0 0 0-5 0 0 27-2 0 0 0-9 0 0 0-2 29 ;B=-15,27,23,0,-20,12,-7,7,10;B=B.;%disp(A)%disp(B)C=A,B;n=size(A,1);ra=rank(A);rc=rank(C);x=zeros(1,n);L=eye(n,n);if ra=rc disp(no solution)elsefor i=1:1:(n-1)for j=i+1:1:n temp=(-C(j,i)/C(i,i);for k=i:1:(n+1)C(j,k)=C(i,k)*temp+C(j,k);end L(j,i)=-(1*temp);endendend for i=n:-1:1if i=n disp(C(i,(n+1)disp(C(i,(n+1)-sum(C(i,(i+1):n).*x(i+1):n)disp(C(i,i)end x(i)=(C(i,(n+1)-sum(C(i,(i+1):n).*x(i+1):n)/C(i,i);end%disp(solution)%disp(x)disp(C)disp(C)disp(L)disp(L)结果:L1.0000e+0000000000-4.1935e-0011.0000e+000000000000-3.0459e-0011.0000e+00000000000-3.5387e-0011.0000e+00000000000-3.9755e-0011.0000e+00000000000-1.5694e-0011.0000e+00000000000-6.5398e-0011.0000e+000000000-1.1210e-001-1.7545e-002-2.4619e-0021.0000e+0000000-1.1927e-001-8.3391e-002-1.4227e-002-1.9962e-002-9.2316e-0021.0000e+000U3.1000e+001-1.3000e+001000-1.0000e+00100002.9548e+001-9.0000e+0000-1.1000e+001-4.1935e+000000002.8259e+001-1.0000e+001-3.3504e+000-1.2773e+0000000007.5461e+001-3.1186e+001-4.5200e-00100-9.0000e+00000004.4602e+001-7.1797e+0000-5.0000e+000-3.5780e+0000000-8.8818e-0164.5873e+001-3.0000e+001-7.8472e-001-5.6154e-00100000-3.5527e-0152.1381e+001-5.1319e-001-3.6724e-00100000002.6413e+001-2.4200e+000000000002.7390e+001LUP 分解代码：结果L=1.0000e+00000000000-4.1935e-0011.0000e+00000000000-3.0459e-0011.0000e+00000000000-3.5387e-0011.0000e+00000000000-3.9755e-0011.0000e+00000000000-1.5694e-0011.0000e+00000000000-6.5398e-0011.0000e+000000000-1.1210e-001-1.7545e-002-2.4619e-0021.0000e+0000000-1.1927e-001-8.3391e-002-1.4227e-002-1.9962e-002-9.2316e-0021.0000e+000U=3.1000e+001-1.3000e+001000-1.0000e+00100002.9548e+001-9.0000e+0000-1.1000e+001-4.1935e+000000002.8259e+001-1.0000e+001-3.3504e+000-1.2773e+0000000007.5461e+001-3.1186e+001-4.5200e-00100-9.0000e+00000004.4602e+001-7.1797e+0000-5.0000e+000-3.5780e+000000004.5873e+001-3.0000e+001-7.8472e-001-5.6154e-0010000002.1381e+001-5.1319e-001-3.6724e-00100000002.6413e+001-2.4200e+000000000002.7390e+001P=100000000010000000001000000000100000000010000000001000000000100000000010000000001求 A 的逆矩阵：通过 LU 分解获得 LU，LUX=I，LY=I，UX=Y 用高斯下去法解得 Y 和 Xfor i=1:n RA,RB,n,X=gaus(L,I(1:n,i);Y(1:n,i)=X;endfor i=1:n RA,RB,n,X=gaus(U,Y(1:n,i);A_(1:n,i)=X;endA_3.8952e-0021.5961e-0025.8083e-0033.6407e-0037.2823e-0031.2867e-0021.4396e-0031.2292e-0031.5863e-0023.7828e-0021.3083e-0026.5129e-0031.2133e-0027.1149e-0032.4090e-0032.1874e-0034.8698e-0031.1613e-0023.8126e-0027.7386e-0036.9188e-0032.8373e-0031.4666e-0032.5028e-0038.1938e-0041.9539e-0036.4150e-0031.8128e-0021.0528e-0022.3922e-0032.3786e-0035.7899e-0034.5601e-0041.0874e-0033.5701e-0031.0089e-0022.4339e-0025.1100e-0034.7635e-0033.4595e-0031.2743e-0043.0388e-0049.9769e-0042.8193e-0036.8017e-0033.0639e-0021.3312e-0039.6677e-0049.3244e-0052.2235e-0047.3002e-0042.0629e-0034.9768e-0034.6809e-0029.7403e-0047.0739e-0041.0381e-0042.4755e-0048.1276e-0042.2967e-0034.7736e-0031.0064e-0033.8169e-0023.3451e-0032.6145e-0046.2346e-0042.0469e-0035.7843e-0033.5966e-0038.1180e-0043.3705e-0033.6510e-002|A|=|U|=6.1817e+0135.编制程序求解矩阵 A 的 cholesky 分解，并用程序求解方程组Ax=b，其中,%5-choleky%A=2 1-5 1 1-5 0 71.7585e-0029.7237e-0033.8776e-0033.2693e-0036.9837e-0034.1874e-0023.0639e-0021.3755e-0031.1095e-003 0 2 1-1 1 6-1-4 ;B=13,-9,6,0;%disp(A)dim=size(A,1);L=zeros(dim,dim);y=x=zeros(dim,1);for i=1:1:dimfor j=1:1:dimif i=(j+1)L(i,j)=(A(i,j)-sum(L(i,1:(j-1).*L(j,1:(j-1)/L(j,j);else L(j,j)=sqrt(A(j,j)-sum(L(j,1:j-1).*L(j,1:j-1);endendendL_T=L.;for i=1:1:dimif i=1 y(1)=B(1)/L(1,1);else y(i)=(B(i)-sum(L(i,1:(i-1).*y(1:(i-1)/L(i,i);endendfor i=dim:-1:1if i=dim x(i)=y(i)/L(i,i);else x(i)=(y(i)-sum(L(i+1):dim,i).*x(i+1):dim)/L(i,i);end%disp()%disp(x(i)enddisp(x)disp(x)结果:x52.2500-38.750030.7500-52.75006.用追赶法编制程序求解方程组Ax=b，其中，代码:%t6-%A=4 2 0 0 3-2 1 0 0 2 5 3 0 0-1 6 ;B=6 2 10 5;B=B.;dim=size(A,1);c=zeros(1,(dim-1);a=zeros(1,dim);b=zeros(1,dim);f=zeros(dim,1);y=zeros(dim,1);%computation for matrix L and Uc(1)=A(1,2);b(1)=A(1,1);for i=2:(dim-1)c(i)=A(i,(i+1);a(i)=A(i,(i-1)/b(i-1);b(i)=A(i,i)-a(i)*A(i-1),i);enda(dim)=A(dim,(dim-1)/b(dim-1);b(dim)=A(dim,dim)-a(dim)*A(dim-1),dim);%end of computation L and U%Ax=f A=LU LUx=f%Ly=fy(1)=B(1);for i=2:dim y(i)=B(i)-a(i)*y(i-1);enddisp(y)disp(y)%end of Ly=ff(dim)=y(dim)/b(dim);for i=(dim-1):-1:1 f(i)=(y(i)-A(i,(i+1)*f(i+1)/b(i);enddisp(f)disp(f)%Yx=f%end of Ux=y%disp(c)%disp(c)%disp(a)%disp(a)%disp(b)%disp(b)%disp(f)%disp(f)结果:f11117.已知.编程求解矩阵 A 的 QR 分解；代码function Q,R=qrhouseholder(A)dim=size(A,1);R=A;Q=eye(dim);for i=1:(dim-1)x=R(i:dim,i);y=1;zeros(dim-i,1);Ht=householder(x,y);H=blkdiag(eye(i-1),Ht);Q=Q*H;R=H*R;endendfunction H,rho=householder(x,y)x=x(:);y=y(:);if length(x)=length(y)error(The Column Vectors X and Y Must Have The Same Length!);endrho=-sign(x(1)*norm(x)/norm(y);y=rho*y;v=(x-y)/norm(x-y);I=eye(length(x);H=I-2*v*v;end结果:C=1.0000 1.0000 0 0 -1.0000 3.0000 -0.5000 0.5000 -2.0000 2.0000 1.5000 0.5000 -2.0000 2.0000 -0.5000 2.5000 Q,R=qrhouseholder(C)Q=-0.3162 -0.7071 -0.2582 0.5774 0.3162 -0.7071 0.2582 -0.5774 0.6325 -0.0000 -0.7746 -0.0000 0.6325 -0.0000 0.5164 0.5774R=-3.1623 3.1623 0.4743 2.0555 0.0000 -2.8284 0.3536 -0.3536 0.0000 0.0000 -1.5492 1.0328 -0.0000 -0.0000 -0.0000 1.15478.分别应用 Jacobi 迭代法和 Gauss-Seidel 迭代法求解如下方程组8:2x1+x2+5x3=15Jacobi 代码A=4 1 1 4 8 1 2 1 5B=7 21 15;B=B.;dim=size(A,1);x=zeros(1,dim);px=zeros(1,dim)disp(x)disp(x)%disp(D)%disp(D)n=input(请输入迭代次数);if n=1for i=1:1:nfor j=1:1:dim temp=1/A(j,j)px(j)=x(j)+temp*(B(j)-sum(A(j,1:dim).*x(1:dim);end x=px;endelse disp(迭代次数不符合规范)enddisp(x)disp(x)请输入迭代次数:10 x0.66231.99492.3284Gauss-seidel 代码A=4 1 1 4 8 1 2 1 5B=7 21 15;B=B.;dim=size(A,1);x=zeros(1,dim);disp(x)disp(x)%disp(D)%disp(D)n=input(please input number of iterations:);for i=1:1:nfor j=1:1:dim temp=1/A(j,j)x(j)=x(j)+temp*(B(j)-sum(A(j,1:dim).*x(1:dim);endenddisp(x)disp(x)结果x0.66672.00002.33339.分别应用Newton迭代法和割线法计算（1）非线性方程2x3 5x 1=0在1,2上的一个根；（2）exsi n x=0在-4,-3上的一个根。Newton 迭代法（1）%t9-newton-iteration%f=2*(x3)-5*x-1;%f_=6*(x2)-5;x0=0;x1=1.1;while 1 x0=x1;if(6*(x02)-5=0)x1=x0-(2*(x03)-5*x0-1)/(6*(x02)-5);else disp(can not use this method)break;endif(abs(x1-x0)0.0005)breakendenddisp(x1)disp(x1)结果x11.6730（2）%t9-newton-iteration%f=2*(x3)-5*x-1;%f_=6*(x2)-5;%f=exp(x)*sin(x)%f_=exp(x)*(sin(x)+cos(x)x0=0;x1=-3.4;while 1 x0=x1;%x1=x0-(2*(x03)-5*x0-1)/(6*(x02)-5);x1=x0-exp(x0)*sin(x0)/(exp(x0)*(sin(x0)+cos(x0);if(abs(x1-x0)0.00005)breakendenddisp(x1)disp(x1)结果：x1-3.1416Newton 截断法（1）代码%t9-newton-hypotenuse%f=exp(x)*sin(x)%f_=exp(x)*(sin(x)+cos(x)%f=2*(x3)-5*x-1;%f_=6*(x2)-5;%x1=-3.4;%x0=-3.5;x0=1.8;x1=1.7;while 1%x2=x1-exp(x1)*sin(x1)*(x1-x0)/(exp(x1)*sin(x1)-exp(x0)*sin(x0);x2=x1-(2*(x13)-5*x1-1)*(x1-x0)/(2*(x13)-5*x1-1)-(2*(x03)-5*x0-1)if abs(x2-x1)0.00005break;else x0=x1;x1=x2;endenddisp(x2);disp(x2)结果x21.6730（2）%t9-newton-hypotenuse%f=exp(x)*sin(x)%f_=exp(x)*(sin(x)+cos(x)x1=-3.4;x0=-3.5;while 1 x2=x1-exp(x1)*sin(x1)*(x1-x0)/(exp(x1)*sin(x1)-exp(x0)*sin(x0);if abs(x2-x1)=b disp(wrong number)flag=0;else fa=a*cos(a)-2;fb=b*cos(b)-2;if fa*fb 0 disp(there is no adequate conditions to judge)flag=0else m=(a+b)/2;if(b-a)accuracy flag=0;else fm=m*cos(m)-2;if fa*fm 14.用 2 点，3 点和 5 点 Gauss 积分法分别计算定积分，并与真值作比较。function ql,Ak,xk=guasslegendre(fun,a,b,n,tol)syms xp=sym2poly(diff(x2-1)(n+1),n+1)/(2n*factorial(n);tk=roots(p);Ak=zeros(n+1,1);for i=1:n+1 xkt=tk;xkt(i)=;pn=poly(xkt);fp=(x)polyval(pn,x)/polyval(pn,tk(i);Ak(i)=quadl(fp,-1,1,tol);endxk=(b-a)/2*tk+(b+a)/2;fun=fcnchk(fun,vectorize);fx=fun(xk)*(b-a)/2;ql=sum(Ak.*fx);disp(ql)disp(ql)end结果:真实值:quadl(fun,a,b)ans=4.6740e-0012点:fun=(x)(x.2.*cos(x);clear fun=(x)(x.2.*cos(x);a=0;b=pi/2;n=1;tol=0.00005;ql,Ak,xk=guasslegendre(fun,a,b,n,tol)ql 4.7464e-001ql=4.7464e-0013点:fun=(x)(x.2.*cos(x);clearfun=(x)(x.2.*cos(x);a=0;b=pi/2;n=2;tol=0.00005;ql,Ak,xk=guasslegendre(fun,a,b,n,tol)ql 4.6724e-001ql=4.6724e-0015点:fun=(x)(x.2.*cos(x);clearfun=(x)(x.2.*cos(x);a=0;b=pi/2;n=4;tol=0.00005;ql,Ak,xk=guasslegendre(fun,a,b,n,tol)ql 4.6740e-001ql=4.6740e-00115.已知常微分方程，分别用 Euler 法，改进的 Euler 法，Runge-Kutta法去求解该方程，步长选为0.1,0.05,0.01.画图观察求解效果。Euler 法代码：%15-euler%du/dx=(2/x)*u+x2*exp(x)%u(1)=0%x in 1,2%0.1 0.05 0.01 a=1;b=2;u(1)=0;len=input(length of step)n=input(please a number)N=(b-a)/len;if n N disp(wrong number)elsefor j=1:1:n x(j)=a+(j-1)*len;u(j+1)=u(j)+(2/x(j)*u(j)+x(j)2*exp(x(j);y(j+1)=u(j+1);endend disp(len(i)disp(len)disp(u)disp(u(n+1)步长：0.1length of step0.1len=1.0000e-001please a number10n=10len(i)1.0000e-001u 1.2260e+004 clearlength of step0.05len=5.0000e-002please a number20n=20len(i)5.0000e-002u 5.8700e+007length of step0.01len=1.0000e-002please a number90n=90len(i)1.0000e-002u 4.3690e+034步长依次为：0.1,0.05,0.01改进的 euler 法代码%15-improve_euler%du/dx=(2/x)*u+x2*exp(x)%u(1)=0%x in 1,2%0.1 0.05 0.01 a=1;b=2;u(1)=0;len=input(length of step)n=input(please a number)N=(b-a)/len;if n N disp(wrong number)elsefor j=1:1:n x(j)=a+(j-1)*len;u_=u(j)+(2/x(j)*u(j)+x(j)2*exp(x(j);u(j+1)=u(j)+(len/2)*(2/x(j)*u(j)+x(j)2*exp(x(j)+(2/(a+j*len)*u_+(a+j*len)2*exp(a+j*len);y(j+1)=u(j+1);endend disp(len(i)disp(len)disp(u)disp(u(n+1)length of step0.1len=1.0000e-001please a number10n=10len(i)1.0000e-001u3.4619e+001 clearlength of step0.05len=5.0000e-002please a number20n=20len(i)5.0000e-002u3.8950e+001 clearlength of step0.01len=1.0000e-002please a number90n=90len(i)1.0000e-002u3.3518e+001 clear 步长：0.01length of step0.01len=1.0000e-002please a number90n=90len(i)1.0000e-002u3.3943e+001步长依次为：0.1,0.05,0.01Runge-kutta 4阶代码:%t15-runge-kutta%du/dx=(2/x)*u+x2*exp(x)%step=input(please input lenght of step)n=input(please input number of iterations)a=1;b=2;N=(b-a)/step;x=zeros(1,n);x(1)=a;u(1)=0;if n N disp(wrong number)elsefor i=1:1:n x(i)=a+(i-1)*step;%tn,un k1=(2/x(i)*u(i)+x(i)2*exp(x(i);%tn+(1/2)*h,un+(1/2)*h*k1 k2=(2/(x(i)+1/2*step)*(u(i)+(1/2)*step*k1)+(x(i)+1/2*step)2*exp(x(i)+1/2*step);%tn+(1/2)*h,un+(1/2)*h*k2 k3=(2/(x(i)+(1/2)*step)*(u(i)+(1/2)*step*k2)+(x(i)+(1/2)*step)2*exp(x(i)+(1/2)*step);%tn+h,un+h*k3 k4=(2/(x(i)+step)*(u(i)+step*k3)+(x(i)+step)2*exp(x(i)+step);u(i+1)=u(i)+(step/6)*(k1+2*k2+2*k3+k4);endenddisp(u)disp(u(n+1)plot(x,u(1:n)结果:please input lenght of step0.1step=1.0000e-001please input number of iterations10n=10u1.8683e+001 clearplease input lenght of step0.05step=5.0000e-002please input number of iterations20n=20u1.8683e+001please input lenght of step0.01step=1.0000e-002please input number of iterations90n=90u1.4323e+001步长依次为:0.1,0.05,0.01