2016矩阵与数值分析上机作业满分(共25页).docx

精选优质文档-倾情为你奉上第一题考虑计算给定向量的范数：输入向量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=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.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范数是.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)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;end%画图，左边是公式算的，右边是算法算的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 problem3()%秦九韶算法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 dex<=n dex=dex+1; temp=temp*x+c(dex); end y(ydex)=temp;endplot(1.95:m:2.05,y);结果：第四题编制计算给定矩阵A的LU分解和PLU分解的通用程序，然后用你编制的程序完成下面两个计算任务：(1) 考虑A=10-1 00 11-1-1 11-1-1 -11Rn×n,自己取定xRn，并计算b = Ax。然后用你编制的不选主元和列主元的Gauss消去法求解该方程组，记你计算出的解为x。对n从5到30估计计算解的精度。(2) 对n从5到30计算其逆矩阵。代码：%LU分解function problem4()clc%生成矩阵Afor n=5:30A=-ones(n);A(:,n)=ones(n,1);for ia=1:n A(ia,ia)=1; for ja=(ia+1):(n-1) A(ia,ja)=0; endendx=rand(n,1);b=A*x;% LU分解n,n=size(A);L=eye(n);M=eye(n);for i=1:n-1 if A(i,i) = 0 fprintf('Error: A(%d,%d)=0!n', i,i); return; end for j=i+1:n L(j,i)=A(j,i)/A(i,i); end for k=i+1:n for l=i:n A(k,l)=A(k,l)-L(k,i)*A(i,l); end endendU=A;M=inv(U)*inv(L)*My1=inv(L)*b;x1=inv(U)*y1;error=norm(x-x1)endreturn;%PLU分解function problem4()clc% 生成矩阵An=5;A=-ones(n);A(:,n)=ones(n,1);for ia=1:n A(ia,ia)=1; for ja=(ia+1):(n-1) A(ia,ja)=0; endendx=rand(n,1);b=A*x;%PLU分解n,n=size(A);Q=eye(n);R=eye(n);S=eye(n);L=zeros(n,n,n-1);for il=1:(n-1) L(:,:,il)=eye(n);endP=zeros(n,n,n-1);for il=1:(n-1) P(:,:,il)=eye(n);endfor i=1:n-1 m=i; a=abs(A(i,i); for k=i+1:n if abs(A(k,i)>a 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(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.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 = 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.8516第五题编制计算对称正定阵的Cholesky分解的通用程序，并用你编制的程序计算Ax = b，其中A=aijRn×n,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.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 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-2uuTRn×n乘以ARn×m的程序HA，注意，你的程序并不显式的计算出H。(3) 考虑矩阵A=123-132-22e 43-102-3027 75/2,用你编制的程序计算H使得HA的第一列为e1的形式，并将HA的结果显示。代码：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结果：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()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 ia<N for iu=(ia+1):N U(ia,iu)=A(ia,iu); end end if ia>1 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; %求误差 for 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.9670wuca1 = 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=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(inS)=xk; if abs(xka-xk)<10(-5)%精度10(-5) break; end xkr=xk; xk=xka;endsubplot(1,2,2);plot(1:inS,f);hold oninSxka结果：inNew = 7xx = 3.4606inS = 10xka = 3.4606结果：第九题用二分法求方程excosx+2=0在区间0,4上的所有根。代码：function problem9()%二分法求解方程clc%查看根的个数N=50;yk=zeros(1,N+1);h=4*pi/N;e=;c=;i=1;for ti=0:h:4*pi e(i)=-2/(exp(ti); c(i)=cos(ti); i=i+1;endplot(0:h:4*pi,e);hold onplot(0:h:4*pi,c);%根的取值范围min=1 4 7 10;max=2 5 8 11;accr=zeros(1,4);acc=0.0001;%精度for r=1:4 while(max(r)-min(r)>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*ymax>0 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=simplify(f);end所以考虑函数f(x)= sin(x)，x 0,1，取二次、三次、四次插值函数。二次：>> x=linspace(0,1,3)>> y=sin(pi*x);>> f2=Newton(x,y)f2 =-4*t*(t - 1)三次：>> x=linspace(0,1,4);>> y=sin(pi*x);>> f3=Newton(x,y)f3 =-(9*3(1/2)*t*(t - 1)/4四次：>>x=linspace(0,1,5);>> y=sin(pi*x);>> f4=Newton(x,y)f4 =(t*(*t3 - *t2 + 32978*t + *2(1/2) + 21487)/40992结果：第十一题对函数fx=11+x2，x -5, 5，取不同的节点数n，用等距节点lagrange插值，观察Runge现象。代码：function p=Lagrange(n,t) x=linspace(-5,5,n);y=1./(1+x.2);m=length(t);for i=1:m s=0; for k=1:n l(k)=1; for j=1:n if j=k l(k)=l(k)*(t(i)-x(j)/(x(k)-x(j); end end s=l(k)*y(k)+s; end p(i)=s;end end>> t=-5:0.1:5;>> p=1./(1+x.2);>> plot(x,p)>> hold on>> p5=Lagrange(5,x);>> plot(x,p5,'r')>> p10=Lagrange(10,x);>> plot(x,p10,'g')>> p15=Lagrange(15,x);>> plot(x,p15,'m')>> p20=Lagrange(20,x);>> plot(x,p20,'y')结果：第十二题令fx=e3x cos(x)，考虑积分02f(x)dx。区间分为50，100,200, 500，1000等，分别用复合梯形以及复合Simpson积分公式计算积分值，将数值积分的结果与精确值比较，列表说明误差的收敛性。代码：function problen12()clcN=input('请输入N=');%改变区间份数h=2*pi/N;%复化梯形公式x=0:h:2*pi;sum=0;for in=2:N sum=sum+exp(3*x(in)*cos(pi*x(in);endTn=2*pi*(1+2*sum+exp(6*pi)*cos(pi*2*pi)/(2*N)%Simapsonsum1=0;sum2=sum;for si=1:N xx=(h/2)+(si-1)*h; sum1=sum1+exp(3*xx)*cos(pi*xx);endSn=2*pi*(1+4*sum1+2*sum2+exp(6*pi)*cos(pi*2*pi)/(6*N)结果：请输入N=50Tn = 3.5125e+07Sn = 3.5231e+07请输入N=100Tn = 3.5205e+07Sn = 3.5232e+07请输入N=200Tn = 3.5226e+07Sn = 3.5232e+07请输入N=500Tn = 3.5231e+07Sn = 3.5232e+07请输入N=1000Tn = 3.5232e+07Sn = 3.5232e+07第十三题分别用2点，3点以及5点的Gauss型积分公式计算如下定积分：(1)-11x21-x2dx, (2) 02sinxxdx代码：function y=Chebyshev(f,n) %使用 Chebyshev 多项式进行 Gauss 型数值积分 %输入： f为表达式， n 为 Gauss 点数 %输出： y 为积分值x=zeros(1,n);A=zeros(1,n);for i=1:n x(i)=cos(2*n-1)*pi/(2*n); A(i)=pi/n;endy=sum(A(1:n).*subs(f,x(1:n);y=double(y);endfunction y=Lengder(f,a,b,x,w) %使用 Lengder 多项式进行 Gauss 积分 %输入： f 为积分表达式， a 和 b 为积分上下限， x 为高斯零点向量， w 为求积系数向量 %输出： y 为积分值T=(a+b)/2+(b-a)/2*x;y=(b-a)/2*sum(w.*subs(f,T);y=double(y);end函数调用：clcformat longsyms xf1=x2;f2=sin(x)/x;a=0;b=pi/2;x1=-0.5574,0.5574; %两节点高斯根向量w1=1,1;x2=-0.7746,0,0.7746; %三节点高斯根向量w2=0.5556,0.8889,0.5556;x3=-0.9062,-0.5385,0,0.5385,0.9062; %五节点高斯根向量w3=0.2369,0.4786,0.5689,0.4786,0.2369;y1=Chebyshev(f1,2),Chebyshev(f1,3),Chebyshev(f1,5)y2=Lengder(f2,a,b,x1,w1),Lengder(f2,a,b,x2,w2),Lengder(f2,a,b,x3,w3)real1=double(int(f1/sqrt(1 -x2),-1,1)real2=double(int(f1/sqrt(1 -x2),-1,1)结果：y1 = 1.4897 1.4897 1.4897y2 = 1.2214 1.6270 1.0828real1= 1.4897real2= 1.4488第十四题考虑微分方程初值问题：dxdt=1(t+1)2ty-y2x0=2分别用Euler法，改进的Euler法，Runge-Kutta法求解该方程。分别取步长为0.1，0.01，0.001，计算到x(1)，画图说明结果。根据题意，首先计算精确积分的结果。得x=ylnt+1+y+y2t+1。将x0=2代入得到y=1或y=-2代码：function problem14()clch=0.1; %步长，改变步长改变此即可N=1/h;x=0:h:1;%Euler法ye=zeros(1,N+1);ye(1)=2;for ie=2:(N+1) ye(ie)=ye(ie-1)+h*(x(ie-1)*ye(ie-1)-ye(ie-1)*ye(ie-1)/(x(ie-1)+1)2);end%画图，红色plot(x,ye,'r');hold on;%改进的Euler法yce=zeros(1,N+1);yce(1)=2;for ice=2:(N+1) f1=(x(ice-1)*yce(ice-1)-yce(ice-1)*yce(ice-1)/(x(ice-1)+1)2; u=yce(ice-1)+h*f1; f2=(x(ice)*u-u*u)/(x(ice)+1)2; sum=f1+f2; yce(ice)=yce(ice-1)+sum*h/2;end%画图,黄色plot(x,yce,'y');hold on;%Runge-Kutta法yr=zeros(1,N+1);yr(1)=2;for ir=2:(N+1) k1=( x(ir-1)*yr(ir-1)- yr(ir-1)*yr(ir-1) )/( x(ir-1) +1 )2; k2=(x(ir-1)+h/2)*(yr(ir-1)+k1/2)-(yr(ir-1)+k1/2)*(yr(ir-1)+k1/2)/(x(ir-1)+h/2)+1)2; k3=(x(ir-1)+h/2)*(yr(ir-1)+k2/2)-(yr(ir-1)+k2/2)*(yr(ir-1)+k2/2)/(x(ir-1)+h/2) +1 )2; k4=( (x(ir-1)+h)*(yr(ir-1)+k3)-(yr(ir-1)+k3)*(yr(ir-1)+k3) )/( (x(ir-1)+h) +1 )2; yr(ir)=yr(ir-1)+h*(k1+2*k2+2*k3+k4)/6;end%画图plot(x,yr);legend('Euler','改进Euler','Runge-Kutta')结果：步长h=0.1时步长h=0.01时步长h=0.001时专心-专注-专业