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

1、大连理工大学矩阵与数值分析报告上机作业大 连 理 工 大 学矩阵与数值分析上机作业课程名称： 矩阵与数值分析 研究生某某： 交作业日时间：2016 年 12 月 20日 第1题1.1程序：Clear all;n=input(请输入向量的长度n：)for i=1:n;v(i)=1/i;endY1=norm(v,1)Y2=norm(v,2)Y3=norm(v,inf)1.2结果n=10 Y1 =2.9290 Y2 =1.2449 Y3 =1n=100 Y1 =5.1874 Y2 =1.2787 Y3 =1n=1000 Y1 =7.4855 Y2 =1.2822 Y3 =1N=10000 Y1 =9

2、.7876 Y2 =1.2825 Y3 =11.3 分析一X数逐渐递增，随着n的增加，X数的增加速度减小；二X数随着n的增加，逐渐趋于定值，无群X数都是1.第2题2.1程序clear all;x(1)=-10-15;dx=10-18;L=2*103;for i=1:L y1(i)=log(1+x(i)/x(i); d=1+x(i);if d = 1 y2(i)=1;else y2(i)=log(d)/(d-1);end x(i+1)=x(i)+dx;endx=x(1:length(x)-1);plot(x,y1,r);hold onplot(x,y2);2.2 结果2.3 分析红色的曲线代表未

3、考虑题中算法时的情况，如果考虑题中的算法如此数值大小始终为1，这主要是由于大数加小数的原因。第3题3.1 程序clear all;A=1 -18 144 -672 2016 -4032 5376 -4608 2304 -512;x=1.95:0.005:2.05;for i=1:length(x); y1(i)=f(A,x(i); y2(i)=(x(i)-2)9;endfigure(3);plot(x,y1);hold on;plot(x,y2,r);F.m文件function y=f(A,x)y=A(1);for i=2:length(A); y=x*y+A(i);end;3.2 结果第4题

4、4.1 程序clear all;n=input(请输入向量的长度n：)A=2*eye(n)-tril(ones(n,n),0);for i=1:n A(i,n)=1;endn=length(A);U=A;e=eye(n);for i=1:n-1 max_data,max_index=max(abs(U(i:n,i); e0=eye(n); max_index=max_index+i-1; U=e0*U; e1=eye(n);for j=i+1:n e1(j,i)=-U(j,i)/U(i,i);end U=e1*U; Pi=e0;%把变换矩阵存到P中 Li=e1; e=e1*e0*e;endfo

5、r k=1:n-2 Ldotk=Lk;for i=k+1:n-1 Ldotk=Pi*Ldotk*Pi;endendLdotn-1=Ln-1;LL=eye(n);PP=eye(n);for i=1:n-1 PP=Pi*PP; LL=Ldoti*LL;endb=ones(n,2);b=e*b; %解方程x=zeros(n,1);x(n)=b(n)/U(n,n);for i=n-1:-1:1 x(i)=(b(i)-U(i,:)*x)/U(i,i);endX=U-1*e-1*eye(n);%计算逆矩阵AN=X;result2n-4,1=AN;result1n-4,1=x;fprintf(%d:n,n)

6、fprintf(%d ,AN);4.2 结果n=51.0625-0.875-0.75-0.5-0.06250.06251.125-0.75-0.5-0.06250.06250.1251.25-0.5-0.06250.06250.1250.251.5-0.0625-0.0625-0.125-0.25-0.50.0625n=101.0625-0.875-0.75-0.5-0.06251.0625-0.875-0.75-0.5-0.06250.06251.125-0.75-0.5-0.06250.06251.125-0.75-0.5-0.06250.06250.1251.25-0.5-0.06250

7、.06250.1251.25-0.5-0.06250.06250.1250.251.5-0.06250.06250.1250.251.5-0.0625-0.0625-0.125-0.25-0.50.0625-0.0625-0.125-0.25-0.50.06251.0625-0.875-0.75-0.5-0.06251.0625-0.875-0.75-0.5-0.06250.06251.125-0.75-0.5-0.06250.06251.125-0.75-0.5-0.06250.06250.1251.25-0.5-0.06250.06250.1251.25-0.5-0.06250.06250

8、.1250.251.5-0.06250.06250.1250.251.5-0.0625-0.0625-0.125-0.25-0.50.0625-0.0625-0.125-0.25-0.50.0625同样的方法可以算出n=20,n=30时的结果，这里就不罗列了。第5题5.1 程序clear all;n=input(请输入向量的长度n：10至20)for i=1:nfor j=1:n a(i,j)=1/(i+j-1);endendfor j=1:n sum=0;for k=1:j-1 sum = sum + l(j,k)2;end l(j,j)=sqrt(a(j,j)-sum);for i=j+1

9、:n sum=0;for k=1:j-1 sum =sum + l(i,k)*l(j,k);end l(i,j)=(a(i,j)-sum)/l(j,j);endendb=ones(n,1);y=zeros(n,1);y(n)=b(n)/l(n,n);for i=n-1:-1:1 y(i)=(b(i)-l(i,:)*y)/l(i,i);endl=l;x=zeros(n,1);x(n)=y(n)/l(n,n);for i=n-1:-1:1 x(i)=(y(i)-l(i,:)*x)/l(i,i);endfprintf(%dt,x);fprintf(n);5.2 结果n=10n=11n=12n=13n

10、=14n=15n=16n=17n=18n=19n=20-746517.86883111493.423-11884558.85478355909.6497329749.5519377549.9445748685378885969.7341571897.81020943960994971030.669820595.08-354658479.91634110905-77484610111-80547115920-82914903907-72327488690-62111481037-60010758946-1.76915E+11-1.68795E+11-15874441979875543407-54

11、9934615923.06265E+123.18327E+123.23484E+122.86169E+122.48062E+122.55037E+127.40163E+126.93865E+12152*-1.17236E+117.93546E+11-5.20373E+13-5.40791E+13-5.42806E+13-4.8713E+13-4.26882E+13-4.65382E+13-1.3123E+14-1.21156E+14-761846200487.35352E+11-6.11155E+124.7524E+144.93812E+144.89563E+144.46792E+143.97

12、243E+144.57473E+141.22E+151.11E+152.18036E+11-2.70378E+122.80302E+13-2.62E+15-2.72E+15-2.66E+15-2.48E+15-2.25E+15-2.72E+15-6.64E+15-5.98E+15-3.70513E+116.12295E+12-8.11E+139.24E+159.59E+159.27E+158.90E+158.30E+151.04E+162.17E+161.94E+163.69292E+11-8.64269E+121.51789E+14-2.17E+16-2.25E+16-2.14E+16-2.

13、14E+16-2.07E+16-2.60E+16-4.18E+16-3.74E+16-1.99261E+117.40507E+12-1.83339E+143.40E+163.53E+163.34E+163.52E+163.56E+164.28E+164.06E+163.77E+1644905979430-3.52275E+121.3792E+14-3.55E+16-3.68E+16-3.47E+16-3.89E+16-4.12E+16-4.38E+16-4.89E+15-1.01E+167.13565E+11-5.87483E+132.35E+162.44E+162.37E+162.73E+1

14、62.94E+162.40E+16-1.64E+16-5.28E+151.08203E+13-8.98E+15-9.29E+15-1.03E+16-9.10E+15-7.78E+15-4.50E+15-1.23E+16-1.58E+161.50E+151.55E+153.06E+15-2.94E+15-6.68E+151.72E+151.63E+161.04E+162.17168E+12-7.94676E+145.08E+157.29E+15-1.45E+153.84E+163.85E+161.58892E+14-2.50E+15-2.19E+15-8.57E+15-2.64E+16-2.01

15、E+164.80384E+14-3.1296E+141.48E+16-8.62E+16-7.98E+162.304E+14-9.01E+151.40E+171.19E+171.99E+15-8.07E+16-6.35E+161.70E+161.09E+167.5453E+14第6题6.1 程序clear all;A= 1 2 3 4; -1 3 sqrt(2) sqrt(3); -2 2 exp(1) pi; -sqrt(10) 2 -3 7; 0 2 7 5/2;U=f61(A(:,2);HA=f62(U,A);f.m文件function U=f61(x)e1=eye(length(x),1

16、);U=x-sign(x(1)*sqrt(dot(x,x)*e1;U=U./sqrt(dot(U,U);function HA=f62(U,A)HA=A-2*U*U*A;6.2 结果-2.26491106454.7358408697.6958675462.2649110644.44E-16-0.321627306-1.9638167380.1766073764.44E-161.5610545830.677680956-0.9856702844.44E-16-4.1572272464.5360883032.1766073764.44E-165.8427727540.036088303第7题7.1

17、 程序clear all;max=1000;x(1,:)=1 2 3;for i=1:max x(i+1,1)=(-2+x(i,2)+3*x(i,3)/5; x(i+1,2)=(1+x(i,1)-4*x(i,3)/2; x(i+1,3)=(10+3*x(i,1)-4*x(i,2)/15; err(i)=sqrt(dot(x(i+1,:)-x(i,:),x(i+1,:)-x(i,:);if(err(i) 10-6)break;endendfigure(7);plot(err);clear err;x(1,:)=1 2 3;for i=1:max x(i+1,1)=(-2+x(i,2)+3*x(i

18、,3)/5; x(i+1,2)=(1+x(i+1,1)-4*x(i,3)/2; x(i+1,3)=(10+3*x(i+1,1)-4*x(i+1,2)/15; err(i)=sqrt(dot(x(i+1,:)-x(i,:),x(i+1,:)-x(i,:);if(err(i) 10-6)break;endendhold onplot(err,r);7.2 结果误差越来越小。第8题8.1 程序clear all;max=100;x(1)=1;for i=1:max x(i+1)=x(i)-(x(i)3+2*x(i)2+10*x(i)-100)/(3*x(i)2+4*x(i)+10);if(abs(x

19、(i+1)-x(i) 10-6)break;endendfigure(8)plot(x);clear x;x(1)=0;x(2)=1;for i=2:max x(i+1)=x(i)-(x(i)3+2*x(i)2+10*x(i)-100)/(x(i)3+2*x(i)2+10*x(i)-(x(i-1)3+2*x(i-1)2+10*x(i-1)*(x(i)-x(i-1);if(abs(x(i+1)-x(i) 10-6)break;endendhold onplot(x,r);8.2结果8.3 分析由计算结果可知，弦截法的收敛速度比牛顿法的收敛速度快。第9题9.1 程序clear all;f(0,4*

20、pi);f.m文件function = f( l,r )if (r-l10-6) fprintf(%g,(r-l)/2+l);returnendif f9(l)*f9(r) 0if f9(l+r)/2+l)*f9(l) 0 f(l,(l+r)/2+l);endif f9(l+r)/2+l)*f9(r) 0 f(l,(l+r)/2+l); f(l+r)/2+l,r);end9.2 结果X=4.71239、8.24668、17.6715、14.1372第10题10.1 程序clear all;n=3;% 节点个数Xj=0:1/n:1;y=sin(pi*Xj);for i=1:n+1 f(i,1)=

21、y(i);endfor j=2:n+1for i=1:n-j+2; dx=(j-1)/n); f(i,j)=(f(i+1,j-1)-f(i,j-1)/dx;endendfor i=1:n+1 a(i)=f(1,i);endx=0:0.001:1;for i=1:1/0.001+1; y1(i)=sin(pi*x(i); y2(i)=f10(a,Xj,x(i);endfigure(10);plot(x,y1,r);hold on;plot(x,y2,);10.2结果10.3分析有图像可知插值函数的值已经很接近原函数的值了。第11题11.1 程序clear all;n=input(请输入n：) %

22、 n代表节点Xj=-5:1/n:5;Yj=1./(1.+Xj.2);x=-5:0.01:5;for i=1:10/0.01+1; y1(i)=1/(1+x(i)2); y2(i)=f(Yj,Xj,x(i);endfigure(11);plot(x,y1,r);hold on;plot(x,y2,);f.m文件function y=f(Yj,Xj,x)y=0;for i=1:length(Yj) l=1;for j=1:length(Xj)if i=j continue;end l=l*(x-Xj(j)/(Xj(i)-Xj(j);end y=y+Yj(i)*l;end11.2 结果从左往右n依次

23、取1、2、3、4 11.3 分析随着n的不断增加，插值越来越接近真实值。第12题12.1 程序clear all;n=input(请输入n：) %n=50 100 200 500 1000x=0:2*pi/n:2*pi;for i=1:n+1 X=x(i); Y(i)=exp(3*X)*cos(pi*X);endY(1)=1;for i=1:n X=(x(i)+x(i+1)/2; Y2(i)=exp(3*X)*cos(pi*X);endT=(x(n+1)-x(1)/(2*n)*(Y(1)+2*sum(Y(2:n)+Y(n+1);S=(x(n+1)-x(1)/(6*n)*(Y(1)+4*sum(

24、Y2)+2*sum(Y(2:n)+Y(n+1);13.2 结果n501002005001000T3.512534119493697e+073.520489119998542e+073.522553497836321e+073.523136936592422e+073.523220478232005e+07S3.523140786833490e+073.523241623782247e+073.523247916807637e+073.523248325445201e+073.523248335509114e+07第13题13.1 程序clear all;G2=f(1/sqrt(3)+f(-1/

25、sqrt(3);G3=0.555555556*f(-0.7745966692)+0.555555556*f(+0.7745966692)+0.88888888889*f(0);G5=0.2369268851*f(-0.9061798459)+0.2369268851*f(+0.9061798459)+. 0.4786286705*f(-0.5384693101)+0.4786286705*f(+0.5384693101)+.0.5688888889*f(0);f.m文件1function y=f(x)y=x2/sqrt(1-x2);f.m文件2function y=f(x)x=pi/2*(x+

26、1)/2;y=sin(x)/x;13.2结果1y=x2/sqrt(1-x2)2y=sin(x)/xG20.816496*G21.74487172393660G31.05409255403515G31.74531021550858G51.05409255403515G51.74530859901330第14题14.1 程序clear all;U0=2;%初值step=0.2;%step表示步长U1(1)=U0;Len=1;for i=1:floor(Len/step) U1(i+1)=U1(i)+step*f14(U1(i),(i-1)*step);end% 改良的Euler法U2(1)=U0;for i=1:floor(Len/step) U2(i+1)=U2(i)+step/2*(f14(U2(i),(i-1)*step)+f14(U1(i+1),(i+1-1)*step);end% Runge_KUttaU3(1)=U0;for i=1:floor(Len/step) t=(i-1)*step; k1=f14(U3(i),t); k2=f14(U3(i)+step/2*k