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

1、共享知识分享快乐大连理工大学矩阵与数值分析上机作业课程名称：矩阵与数值分析研究生姓名：交作业日时间：2016年12月20日卑微如螃蚁、坚强似大象1*考虑计算给定向量的范数：愉入向量r = (ii,输出|工111, |I应|H|QQ 请编制一个通用程序，并用你编制的程序计算如下向量的范数：=1<)> KX), 1IKK甚至更丸的泛计算其范数、你会发现什么结果？你能否修故你的程序使得计算结栗 相对精确呢？1.1 程序：Clear all;n=input('请输入向量的长度n :') for i=1: n;v(i)=1/i；end丫仁门 orm(v,1)Y2=n orm(

2、v,2)Y3=no rm(v,i nf)1.2 结果n=10 Y1 =2.9290Y2 =1.2449Y3=1n=100 Y1 =5.1874Y2 =1.2787Y3=1n=1000 Y1 =7.4855Y2 =1.2822Y3=1N=10000 Y1 =9.7876Y2 =1.2825Y3=11.3 分析一范数逐渐递增，随着n的增加，范数的增加速度减小；二范数随着n的增加，逐渐趋于定值，无群范数都是1.2 .考虑站=f(H)=吟空 其中定义=匕此时川眄是连续函瓠用此公式计算当工“J旷巴1旷诗时 的函数值，画出图像.另一方面考虑下面算法：d1 + 比if d = 1 thenk1y = lnr

3、f/(d-l)end if用此算法计茸m mo-叫呵时的函数色 画出图像“比较一下发生了什么？2.1 程序clear all;x(1)=-10a-15;dx=10A-18;L=2*10A3;for i=1:Ly1(i)=log(1+x(i)/x(i);d=1+x(i);if d = 1y2(i)=1 ；elsey2 =log(d)/(d-1);endx(i+1)=x(i)+dx;endx=x(1 :le ngth(x)-1);plot(x,y1,'r')；hold onplot(x,y2);2.2 结果2.3 分析红色的曲线代表未考虑题中算法时的情况，如果考虑题中的算法则数值大

4、小始终为 这主要是由于大数加小数的原因。3 .首先编剖一个利用秦九匏算法计算一个多项式在给宣点的禹数值的通用程序”你的程序包括输入多 顼式的系数以及给定点，输出画数值.利用你编制的程序计算-4608X2+ 230 牡-512p(j)= Q-2)B =x9_1>>zb + i 汝7 . 672JA + 201fijs-4032d +在丁二2邨域附ii的值.画出0神厩1.95,20.5上的图像.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:le ng

5、th(x);y1(i)=f(A,x(i);y2(i)=(x(i)-2)7;endfigure(3);plot(x,y1);hold on;plot(x,y2, 'r')；F.m文件fun cti on y=f(A,x)y=A(i)；for i=2:le ngth(A); y=x*y+A(i);end ;3.2 结果第4题4 .编剖计鼻踣定担阵力分祥和向血分解的通用程序T然焉用你编制的程序完 成下面两个什算任 齐： 考虑10+ a 01-1 a * C_1 . -1 -1 1宜己取定上月即并计算b = 找焉f f 1你编制的年逸主元和列主元的C&ts消去法岸解该 方程俎，

6、记你计算出的解初L对八从5到30牯计计算孵的精度.对畀从斤到3(1计算其逆腿阵.4.1 程序 clear all ; n=input('请输入向量的长度 n ： 1) A=2*eye(n)-tril(ones(n,n),0);for i=1:nA(i,n)=1; end n=length(A); L)=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=eO*U; e1=eye(n);for j=i+1 :n e1 (jJ)=-U(j,i)/

7、U(i,i); endU=e1*U;Pi=eO;%把变换矩阵存到P中Li=e1; e=e1 *e0*e; end for k=1 :n-2 Ldotk=Lk;for i=k+1:n-1Ldotk=Pi*Ldotk*Pi;end end Ldotn-1=Ln-1; LL=eye(n); PP=eye(n); for i=1:n-1PP=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:1endX=UA-1*eA-1*eye( n);AN=X'%计算逆矩

8、阵result2 n-4,1=AN;result1 n-4,1=x;fprintf( '%d:n',n)fprintf( 1%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.50.06251.125-0.75-0.5-0.

9、06250.06251.125-0.75-0.50.06250.1251.25-0.5-0.06250.06250.1251.25-0.50.0625-0.06250.125-0.1250.25-0.251.5-0.5-0.06250.06250.0625-0.06250.125-0.1250.25-0.251.5-0.51.06250.0625-0.8751.125-0.75-0.75-0.5-0.5-0.0625-0.06251.06250.0625-0.8751.125-0.75-0.75-0.5-0.50.06250.1251.25-0.5-0.06250.06250.1251.25-

10、0.50.0625-0.06250.125-0.1250.25-0.251.5-0.5-0.06250.06250.0625-0.06250.125-0.1250.25-0.251.5-0.5-0.0625-0.0625-0.0625-0.06250.0625-0.0625-0.0625-0.0625-0.06250.0625同样的方法可以算出n=20,n=30时的结果,这里就不罗列了。第5题5.编制计算对称正定阵的C加伽切分解的通用衽序，并用你编制的程序计昇，二优其中川二 (呵)W即巴呵二士护可以曲你自己取岂对”从10到刘脸证程序的可绑性.5.1 程序clear all ;n=input(&

11、#39;请输入向量的长度n ： 10至20')for i=1:nfor j=1 :na(ij)=1/(i+j-1)；endendfor j=1:nsum=0;for k=1 :j-1sum = sum + l(j,k)A2;endl(jj)=sqrt(a(jj)-sum);for i=j+1 :nsum=0;for k=1 :j-1sum =sum + l(i,k)*l(j,k);endl(ij)=(a(ij)-sum)/l(j,j)；endendb=ones(n,1);y=zeros(n,1);y(n)=b(n)/l(n,n);for i=n-1 >1:1y(i)=(b(i)-

12、l(i,:)*y)/l(i,i); end l=l'x=zeros(n,1);x(n)=y(n)/l(n,n);for i=n-1 >1:1x(i)=(y(i)-l(i,：)*x)/l(i,i);endfprintf('%dt' ,x);fprintf('n');5.2结果n=10n=11n=12n=13n=14n=15n=16n=17n=18n=19n=20-746517.83111493.4-11884558478355909497329749519377549445748685378885969341571897102094396994971

13、03068823.85.6.5.9.7.80.669820595.-35465847163411090-77484610-80547115-82914903-72327488-62111481-60010758-1.76915E-1.68795E089.95111920907690037946+ 11+11-15874441987554340-549934613.06265E+3.18327E+3.23484E+2.86169E+2.48062E+2.55037E+7.40163E+6.93865E+9775921212121212121212152439746-1.17236E7.93546

14、E+-5.20373E-5.40791 E-5.42806E-4.8713E+-4.26882E-4.65382E-1.3123E+-1.21156E43+1111+13+13+ 1313+13+1314+14-76184620 7.35352E+-6.11155E4.7524E+14.93812E+4.89563E+4.46792E+3.97243E+4.57473E+1.22E+151.11E+1504811+12414141414142.18036E+ -2.70378E 2.80302E+-2.94E+15 -6.68E+15 1.72E+151.63E+16 1.04E+16-2.6

15、2E+15 -2.72E+15 -2.66E+15 -2.48E+15 -2.25E+15 -2.72E+15 -6.64E+15 -5.98E+15-1.99261E 7.40507E+-1.83339E3.40E+16 3.53E+16 3.34E+16 3.52E+16 3.56E+16 4.28E+16 4.06E+16 3.77E+16+1112-3.52275E1.3792E+1-3.55E+16 -3.68E+16 -3.47E+16 -3.89E+16 -4.12E+16 -4.38E+16 -4.89E+15 -1.01 E+1630+127.135

16、65E+-5.87483E2.35E+16 2.44E+16 2.37E+16 2.73E+16 2.94E+16 2.40E+16 -1.64E+16 -5.28E+1511+131.08203E+-8.98E+15 -9.29E+15 -1.03E+16 -9.10E+15 -7.78E+15 -4.50E+15 -1.23E+16 -1.58E+16131.50E+15 1.55E+15 3.06E+152.17168E+ -7.94676E5.08E+15 7.29E+15-1.45E+15 3.84E+16 3.85E+1612+141.58892E+-2.50E+15 -2.19E

17、+15 -8.57E+15 -2.64E+16 -2.01 E+16144.80384E+ -3.1296E+1.48E+16 -8.62E+16-7.98E+1614142.304E+14 -9.01 E+15 1.40E+171.19E+171.99E+15 -8.07E+16 -6.35E+161.70E+16 1.09E+167.5453E+111+1213-3.70513E 6.12295E+ -8.11E+13 9.24E+159.59E+15 9.27E+15 8.90E+15 8.30E+15 1.04E+16 2.17E+16 1.94E+16+W 123.69292E+ -

18、8.64269E 1.51789E+-2.17E+16 -2.25E+16 -2.14E+16 -2.14E+16 -2.07E+16 -2.60E+16 -4.18E+16 -3.74E+1611+12146. . 1）编制程序应做/力其作用是对愉人的向量丁，输出单位向量和使f2波力7二iaihh（2）编制Householder变换阵H = / - 2uuT 口心乘以力即“的雷宁次,注意*你的程序并不加戎的计昇 出/.考虑矩阵(1-1A- -2_7101 °23432C7T2凸 727 而 26.1 程序clearall ;A= 1234;-13 sqrt(2) sqrt(3);-

19、22 exp(1) pi;-sqrt(10)2-37;0275/2;U=f61(A(:,2);HA=f62(U,A);f.m文件fun cti on U=f61(x) e1=eye(le ngth(x)51);U=x-sig n(x(1 )*sqrt(dot(x,x)*e1;U=U./sqrt(dot(U,U);fun cti on 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.561

20、0545830.677680956-0.9856702844.44E-16-4.1572272464.5360883032.1766073764.44E-165.8427727540.036088303了,用/ac仍/代求解下面的方程组，输出迭代每一步的课差I血一T,|:5xi* X14- 2X2 + 八3 一12 Jx,3=2一肚+ 4x2+ 15J: J= 107.1程序clear all ;max=1000;x(1,： )=1 2 3;for i=1 :maxx(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+

21、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) < 10A-6) break ;endendfigure ；plot(err);clear err;x(1,：)=1 2 3;for i=1 :maxx(i+1,1)=(-2+x(i,2)+3*x(i,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

22、,:),x(i+1 ,:)-x(i,:);if (err(i)<10A-6) break ;endend hold on plot(err, Y);7.2结果误差越来越小。第8题&取不同的初值用Netonit代以及弦裁法求方程T3-F2八+ 10x- 100 = 0的实根，列 表或者画图说明收敛 速度.8.1程序clear all ;max=100;X =1；for i=1 :maxx(i+1 )=x(i)-(x(i)A3+2*x(i)A2+10*x(i)-100)/(3*x(iF2+4*x(i)+10);if (abs(x(i+1)-x(i)<10A-6)break ;e

23、ndendfigure(8)plot(x);clear x;x(1)=0;x(2)=1;for i=2:maxx(i+1 )=x(i)-(x(i)A3+2*x(i)A2+10*x(i)-100)/(x(i)A3+2*x(i)A2+10*x(i)-(x (i-1)A3+2*x(i-1)A2+10*x(i-1)*(x(i)-x(i-1);if (abs(x(i+1)-x(i)<10A-6)break endendhold onplot(x, 'r');8.2 结果8.3 分析由计算结果可知，弦截法的收敛速度比牛顿法的收敛速度快。第9题9 用二分法求方程严cos x + 2 =

24、 0在区间0, 4TT上的所有根.9.1 程序clear all;f(0,4*pi);f.m文件fun cti on = f( l,r)if (r-l<10A-6)fprintf( '%g,',(r-l)/2+l);returnendif f9(l)*f9(r) < 0if f9(l+r)/2+l)*f9(l)< 0f(l,(l+r)/2+l)；endif f9(l+r)/2+l)*f9(r)< 0 f(l+r)/2+l,r);endendif f9(l)*f9(r) > 0f(l,d+r)/2+l);f(l+r)/2+lj)；end9.2 结果X

25、=4.71239、8.24668 > 17.6715 > 14.1372第10题ID.考虑函=sinfira?),彳 0,1,用等距节点作丁)的帕值.画出插值多埼式以及人可的图像， 观察收线性.10.1程序clear all ;n=3; %节点个数Xj=0:1/n:1;y=si n( pi*Xj);for i=1: n+1f(i,D=y(i);endfor j=2: n+1for i=1: n-j+2;dx=(j-1)/n);endendfor i=1: n+1a(i)=f(1J)；endx=0:0.001:1;for i=1:1/0.001+1;y1(i)=si n(pi*x(i

26、);y2(i)=f10(a,Xj,x(i);endfigure(10);plot(x,y1, 'r');hold on; plot(x,y2J);10.2结果10.3分析有图像可知插值函数的值已经很接近原函数的值了。第11题11 *对函数7E丘55,取不同的节点数趾观Runge现象，用等距节点作lagra “少：桶值t11.1 程序clear all;n=input('请输入n ： ') % n代表节点Xj=-5:1/n:5;Yj=1./(1.+Xj.A2);x=-5:0.01:5;for i=1:10/0.01+1;y1(i)=1/(1+x(i)A2);y2(

27、i)=f(Yj,XjJx(i);endfigure(11);plot(x,y1,'r')；hold on;plot(x,y2 11 )；f.m文件fun cti on y=f(Yj,Xj,x) y=o；for i=1 :le ngth(Yj)1=1；for j=1 :length(Xj) if i=j con ti nue endl=l*(x-Xj(j)/(Xj(i)-Xj(j); endy=y+Yj(i)*l;end11.2 结果从左往右n依次取1、2、3、411.3 分析随着n的不断增加，插值越来越接近真实值。第12题血令£5=£缸COS (7TX )，

28、考虑积分尸/ (忑)血.区闾分为50t 1。八200, 500, 1000等，分 别用更合样形以及赞合引执西皿赧分公式计算积分值*将散值积分的结呆与精确值比较，列表说明误差的收敛性.12.1程序clear all;n=input('请输入 n ： ') %n=50 100 200 500 1000 x=0:2*pi/n:2*pi;for i=1: n+1X=x(i);Y(i)=exp(3*X)*cos(pi*X);endY(1)=1； for i=1: nX=(x(i)+x(i+1)/2;Y2(i)=exp(3*X)*cos(pi*X);endT=(x( n+1 )-x(1 )

29、/(2* n)*(Y(1)+2*sum(Y(2: n) )+Y( n+1);S=(x( n+1 )-x(1 )/(6* n)*(Y(1 )+4*sum(Y2)+2*sum(Y(2: n)+Y( n+1);13.2结果50100T3.512534119493697e+073.520489119998542e+07S3.523140786833490e+073.523241623782247e+0720050010003.522553497836321e+073.523136936592422e+073.523220478232005e+073.523247916807637e+073.52324

30、8325445201e+073.523248335509114e+07第13题13,分别用2点，3点以及5点的伽唤型枳分公式计算如下定枳分:jr/ 2 sin x(2) / dr4 iJo 工13.1程序clear all;G2=f(1/sqrt(3)+f(-1/sqrt(3);G3=0.555555556*f(-0.7745966692)+0.555555556*f(+0.7745966692)+0.888888 88889*f(0);G5=0.2369268851 *f(-0.9061798459)+0.2369268851 *f(+0.9061798459)+.0.4786286705*f(-0.5384693101 )+0.4786286705*f(+0.5384693101)+. 0.568