2016矩阵与数值分析上机作业满分.docx

1、第一题考虑计算给定向量的范数：输入向量x=x1,x2,xnT,输出x1,x2,x。请编制一个通用程序，并用你编制的程序计算如下向量的范数：x=1,12,13,1nT, y=1,2,nT.对n = 10，100，1000甚至更大的n计算其范数，你会发现什么结果?你能否修改你的程序使得计算结果相对精确呢？代码：function Problem1()clcN=input(input N=);X=zeros(1,N);Y=1:N;for i=1:N X(i)=1/i;endx1=0;x2=0;y1=0;y2=0;for i=1:N x1=x1+abs(X(i); x2=x2+X(i)*X(i); y1

2、=y1+abs(Y(i); y2=y2+Y(i)*Y(i);endx1;x2=sqrt(x2);x3=max(abs(X);y1;y2=sqrt(y2);y3=max(abs(Y);fprintf(x的1范数是%.2fn,x1);fprintf(x的2范数是%.2fn,x2);fprintf(x的无穷范数是%.2fn,x3);fprintf(y的1范数是%.2fn,y1);fprintf(y的2范数是%.2fn,y2);fprintf(y的无穷范数是%.2fn,y3);结果：input N=10x的1范数是2.93x的2范数是1.24x的无穷范数是1.00y的1范数是55.00y的2范数是19

3、.62y的无穷范数是10.00input N=100x的1范数是5.19x的2范数是1.28x的无穷范数是1.00y的1范数是5050.00y的2范数是581.68y的无穷范数是100.00input N=1000x的1范数是7.49x的2范数是1.28x的无穷范数是1.00y的1范数是500500.00y的2范数是18271.11y的无穷范数是1000.00第二题考虑y=fx=ln(1+x)x,其中定义f(0)=1，此时f(x)是连续函数。用此公式计算当x-10-15,10-15时的函数值，画出图像。另一方面，考虑下面算法：d=1+xif d=1 theny=1else y=lnd/(d-1

4、)end if用此算法计算x-10-15,10-15时的函数值，画出图像。比较一下发生了什么？代码：function problem2()clcN=1000;n=2*10(-15)/N;%公式计算f=zeros(1,N+1);t=1;for i=-10(-15):n:10(-15) if(i=0) f(t)=log(1+i)/i; else f(t)=1; end t=t+1;end%算法计算a=zeros(1,N+1);t=1;for i=-10(-15):n:10(-15) d=1+i; if(d=1) a(t)=log(d)/(d-1); else a(t)=1; end t=t+1;e

5、nd%画图，左边是公式算的，右边是算法算的subplot(1,2,1);plot(-10(-15):n:10(-15),f);hold onsubplot(1,2,2);plot(-10(-15):n:10(-15),a);hold on结果：第三题首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用程序，你的程序包括输入多项式的系数以及给定点，输出函数值。利用你编制的程序计算px=x-29=x9-18x8+144x7-672x6+2016x5-4032x4+5376x3-4608x2-512在x=2邻域附近的值。画出px在x1.95,20.5上的图像。代码：function pro

6、blem3()%秦九韶算法clcn=9;c=1,-18,144,-672,2016,-4032,5376,-4608,2304,-512;N=100;m=(2.05-1.95)/N;y=zeros(1,N+1);ydex=0;for x=1.95:m:2.05 ydex=ydex+1; dex=1; temp=c(dex); while dexa a=abs(A(k,i); m=k; end end A(i,m,:)=A(m,i,:); P(i,m,:,i)=P(m,i,:,i); for j=i+1:n L(j,i,i)=-A(j,i)/A(i,i); for l=i:n A(j,l)=A(

7、j,l)+L(j,i,i)*A(i,l); end endendU=A;Ufor x=1:n-1 Q=L(:,:,x); for y=x+1:n-1 Q=P(:,:,y)*Q*P(:,:,y); end S=S*inv(Q);end%L矩阵S for z=1:(n-1) R=P(:,:,z)*R;end%P矩阵Rb1=R*b;b2=inv(S)*b1;x2=inv(U)*b2;M=inv(U)*inv(S)*Rerror=norm(x-x2)结果：由于篇幅有限仅列出n=5时的精度及逆矩阵LU分解M = 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0.500

8、0 -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.0625error = 8.8818e-16PLU分解M = 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.5000 0.5000 0.2500 0.1250 0.0625 0.0625error =7.85

9、16第五题编制计算对称正定阵的Cholesky分解的通用程序，并用你编制的程序计算Ax = b，其中A=aijRnn,aij=1i+j-1。b可以由你自己取定，对n从10到20验证程序的可靠性。代码：function problem5()%cholesky分解clc%生成矩阵AN=10; A=zeros(N);for ia=1:N for ja=1:N A(ia,ja)=1/(ia+ja-1); endendL=chol(A);L 结果：由于篇幅有限仅列出n=10时的L矩阵L = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0

10、.1111 0.1000 0 0.2887 0.2887 0.2598 0.2309 0.2062 0.1856 0.1684 0.1540 0.1417 0 0 0.0745 0.1118 0.1278 0.1331 0.1331 0.1304 0.1265 0.1220 0 0 0 0.0189 0.0378 0.0525 0.0630 0.0702 0.0748 0.0777 0 0 0 0 0.0048 0.0119 0.0195 0.0265 0.0326 0.0378 0 0 0 0 0 0.0012 0.0036 0.0068 0.0103 0.0139 0 0 0 0 0 0

11、0.0003 0.0011 0.0022 0.0038 0 0 0 0 0 0 0 0.0001 0.0003 0.0007 0 0 0 0 0 0 0 0 0.0000 0.0001 0 0 0 0 0 0 0 0 0 0.0000第六题(1) 编制程序House(x)，其作用是对输入的向量x，输出单位向量u使得I-2uuTx=x2e1。(2) 编制Householder变换阵H=I-2uuTRnn乘以ARnm的程序HA，注意，你的程序并不显式的计算出H。(3) 考虑矩阵A=123-132-22e 43-102-3027 75/2,用你编制的程序计算H使得HA的第一列为e1的形式，并将HA的

12、结果显示。代码：function problem6()clcx=1;-1;-2;-10(1/2);0;if x(1)0 sg=1;else sg=-1;endNx=length(x);e1=zeros(Nx,1);e1(1)=1;sum=0;for i=1:Nx sum=sum+x(i)*x(i);endx2=sqrt(sum);clear sum;w=x-sg*x2*e1;u=w/sqrt(w*w);u%生成矩阵AA=1,2,3,4;-1,3,2(1/2),3(1/2);-2,-2,exp(1),pi;-10(1/2),2,-3,7;0,2,7,5/2;H=eye(Nx)-2*u*u;H*A

13、结果：u = -0.6124 -0.2041 -0.4082 -0.6455 0ans = 4.0000 -0.8311 1.4090 -6.5378 0.0000 2.0563 0.8839 -1.7805 0.0000 -3.8874 1.6576 -3.8836 0.0000 -0.9843 -4.6770 -4.1078 0 2.0000 7.0000 2.5000第七题用Jacobi和Gauss-Seidel迭代求解下面的方程组，输出迭代每一步的误差xk-x:5x1-x2-3x3=-2-x1+2x2+4x3=1-3x1+4x2+15x3=10代码：function problem7(

14、)clcinter=10;%迭代次数rx=69*7/61-8, 69*3/(61*2)-7/2, 69/61; %根据题意生成A,bA=5 -1 -3; -1 2 4; -3 4 15;b=-2 1 10;N=length(b);x=ones(N,1);%Jacobi迭代法L=zeros(N,N);D=L;U=L;for ia=1:N D(ia,ia)=A(ia,ia); if ia1 for il=1:(ia-1) L(ia,il)=-A(ia,il); end endendfor intt=1:inter x=inv(D)*(L+U)*x+inv(D)*b; wuca1=0; %求误差 f

15、or ix=1:N wuca1= wuca1+(x(ix)-rx(ix)2; end wuca1end%G-S迭代法for intt=1:inter x=inv(D-L)*U*x+inv(D-L)*b; wuca2=0; %求误差 for ix=1:N wuca2= wuca2+(x(ix)-rx(ix)2; end wuca2end结果:wuca1 = 24.6036wuca1 = 12.3259wuca1 = 1.9276wuca1 = 5.4520wuca1 = 15.9787wuca1 = 10.0748wuca1 = 3.6731wuca1 = 6.2770wuca1 = 11.96

16、70wuca1 = 8.9591wuca2 = 5.5096wuca2 = 9.1425wuca2 = 6.9372wuca2 = 8.0757wuca2 = 7.4675wuca2 = 7.7807wuca2 = 7.6170wuca2 = 7.7018wuca2 = 7.6577wuca2 =7.6806第八题取不同的初值用Newton迭代以及弦截法求方程x3+2x2+10x-100=0的实根，列表或者画图说明收敛速度。代码：function problem8()%inNew记录Newton迭代法的收敛速度%inS记录弦截法的收敛次数clcchuzhi=4;%初值f=;%Newton迭代x

17、=chuzhi;for inNew=1:inter xx=x-(x3+2*x2+10*x-100)/(3*x2+4*x+10); f(inNew)=x; if abs(xx-x)10(-5) break; end x=xx;endsubplot(1,2,1);plot(1:inNew,f);hold onclear f;inNewxx%弦截法xk=chuzhi;xkr=chuzhi-1;for inS=1:inter xka=xk-(xk3+2*xk2+10*xk-100)*(xk-xkr)/(xk3+2*xk2+10*xk-100)-(xkr3+2*xkr2+10*xkr-100); f(i

18、nS)=xk; if abs(xka-xk)acc) ymax=exp(max(r)*cos(max(r)+2; temp=(max(r)+min(r)/2; ytemp=exp(temp)*cos(temp)+2; if ytemp*ymax=0 break; else if ytemp*ymax0 max(r)=temp; else min(r)=temp; end end end accr(r)=temp;endaccr结果：accr = 1.8807 4.6940 7.8548 10.9955第十题考虑函数f(x)= sin(x)，x 0, 1。用等距节点作f(x)的Newton插值，画出插值多项式以及f(x)的图像，观察收敛性。代码：function f=Newton(x,y)syms t;n=length(x);c(1:n)=0;f=y(1);y1=0;p=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); p=p*(t-x(i); f=f+c(i)*p; y=y1;endf=simplif