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

1、矩阵与数值分析上机作业学校：大连理工大学学院：班级：姓名：学号：授课老师：注：编程语言Matlab1考虑计算给定向量的范数:输人向量力二(必如，切):输出|X|i9|X|29|k|oo”请编制一个通用程序，并用你编制的程序计算如下向童的范数：对n=10,100,1000甚至更大的几计算其范数，你会发现什么结果？你能否修改你的程序使得计算结果相对精确呢？Norm,m函数functions-Noim(x,m)务求向量次的范数务m12mf12n-length(x);s-0;switchmcase1fori-l:ns-s*abs(x);case2磋范数fori-l:nS-s+x(i，2;s-sqrt(

2、s);caseinf蕤无穷范数s-ma:(abs(x);计算向X,y的范数Testi.mclearall:clc;nl-10;n2-100;n3-1000;xl-1-/l:nlf;x2-1/l:n21;:3-l/1:n31;yl二Ty2-l:n2Ty3-l:n3Tdisp(fn-10时)；dispfx的1范数);disp(Norm(xl.1);dispfx2.范数：):disp(Norm(xl.2);disp(x：I4穷.泞”:)：disp(Noim(xl：inf)；dispfy的1范数f);disp(Norm(ylf1);disp*y妁2-范数：1);disp(Norm(yl,2);disp

3、(y:JU:j-、disp(Noixu(ylzinf),dispfu-100)；dispfx：1-2);disp(Norm(x2f1);dispfx:J);disp(Norm(x2.2);dispfx：；1Z；/：;disp(Norm(:21范数：)；disp(Norm(y2,1;dispC丫的2范数：.)；disp(Norm(y22;diq)(y白；：;disp(Norm(y2:inf);disp(111-1000)；disp(fx1);disp(Norm(x3J);dispfic2-H);disp(Norm(x3r2);disp(xtl；Cj-数：)；disp(Noun(:3.inf):d

4、isp(y的2.范数：)；disp(Nonn(y3：1);dispCy的2范数：);disp(Norm(y3r2);disp(y的无方范数：)；disp(Norm(y3.inf);运行结果：n=10时X的1范数29290;X的2-范数：1.2449;x的无穷范数：1y的卜范数：55;丫的2-范数:19.6214;丫的无穷范数：1。X的无穷范数：1y的无穷范数：won=100时X的1,范数51874;X的2.范数：1.2787;1000时y的1-范数：500500;y的2范数48272+004;y的无穷范教：10002.考虑匕二/()二芈也，其中定义/(0)-b此时/(巧是连续函数.用此-公式计

5、算当-10-叫10-15时的函数值，画出图像。另一方面，考虑下面算法：=1+_Lif(1=1theny=1elsev=lnd/(d-l)endif用此算法计算xe-10A10J5时的函数值，画出图像比较一下发生了什么？Test2.mclearall.clc;n-100;%K|Hh-220A(15)也;当步长x-10A(-15)hl0A(-15)：务第-种原函数fl-zeros(ln+1);foik-l:n+lifx(k)-0fl(k尸1:endendsubplot(221,1)；axis(-l0A(-150Alegend勺J;$第二种算法f2-zeros(l.nl);fork-1:n+ld-l

6、+x(k);f2(k)-log(d)(d-1);elsef2(k)-l;subplot(2,l,2);plot(Xzf2-il); TOC o 1-5 h z axis”-10A(-15):10A(-15N0)：legend);运行结果：11111da11第二种算法111uni0.8.0.6.0.4020020.40.60.8X1严显然第二种算法结果不准确，是因为汁算机中的舍入误差造成的，当xeEOP.io,时，d=kx,计算机进行舍人造成d恒等于1,结果函数值恒为1。3.首先编制一个利用秦九昶算法计算一个多项式在给定点的函数值的通用程序，你的程序包括输人多项式的系数以及给定点，输出函数值.利

7、用你编制的程序计算p(x)=可2)9=x9-18r8-f144x7-672护.2016式一4032,4-5376/-4608x2230牡512在工=2邻域附近的值。画出P(t)在1.95,20.5上的图像.酢：秦九韶算法:QinJS.mfunctionyQinJS(a.K)务丫输出函数值%3多项式系数，由高次到零显给定点、n-length(a);s-a(l);fori-2:ns-s*x+a(i);endY-s;计算p(x):test3.mclearaliiclc;6:0-2:2-的邻域dispCx-2:i算域：，);xa-l-18144-6722016-40325376-46082304-51

8、2;p-zeros(l,5);fori-l:5enddisp(应多项式p立：，)；pxk-L95:0.01:20.5;nk-lengtli(xk);pk-zeros(Lnk);k-l;fork-l:nkpk(k)-QinJS(a：xk(k);endplot(xkzpk*-r*);xlabel(H);ylabel(:p(x);运行结果：x=2的邻域:1.60001.80002.00002.20002.4000相应多项式P值：p=1-0e-003*-0.2621-0.000500.00050.26214.编制计算给定矩阵4的么U分解和PEU分解的通用程序，然后用你编制的程序完成下面两个计算任务：考

9、虑自己取定并计算b=Ar-然后用你编制的不选主元和列主元的Gauss消去法求解该方程组，记你计算出的解为筑对72从5到30估计计算解的精度。(2)对几从5到30计算其逆矩阵LU分解,LUDecommfunctionL,U-LUDecom(A)务不带列主元的LU，)解N-size(A);n-N(l);L-eye(n);U-zeros(n);fori-l:nU(lzi)-A(lj);L(iJ)-A(kl)/U(1J);endfori-2:nfoxj-in2-0;fork-l:Mz-z+L(kk)*U(k);endU(i.j)-A(ij)-z;endforz-0;foikTilz-z+L0 xk)U

10、(ki);endL(jzi)-(A(jsi)-z)/U(i,i);endendPLU分解PLUDecom.mfunctionP.LzU-PLUDecom(A)务带列主元的LU分解mzm-size(A)；U-A;P-eye(m);L-eye(m);foii-l:mforj-i:mfork-l:Mendenda-i;b-abs(t(i);forj-i十上mifb-eps)xO-x;x-PxO-fn-n*l;err-norm(x-xO.inf)dispCWarning:迭代次数太多可能不收敛?，)；return;endendGauss_Seidel迭代：Gauss_Seidelmfunctionxz

11、n-Gauss_Seidel(A.b.xO)%-Gauss-Seidel性方程组方程组系数阵A方程组右端项b务一初始值xO求解炎求的精确度eps迭代步数控制M务返回求得的解x返回迭代步数neps-1-。总-5;毛求人的对角矩阵电求人的下三角阵务求离的上三角阵M-10000;D-diag(diag(A);L-tnl(Arl)；U-triu(A.l);G-(D-L)U;f-(D-L)b:G*xO+fwhile(err-eps)n-1 ；err-norm(x-xO.mf)毛迭代次数xO-x;x-G*xO+fn-n-1;err-Qorm(x-xO=inf)M)disp(、Vanung:迭代次数太多.可

12、能不收敛!*)；return:解方程组.test7.mclearall.clc;A-5-l-3;-124;-3415;b-2;l;10;disp(-ini.i：x-AbdispfJ始值)；x0-0;0;0disp(*Jacobi代过fil*);xj,nj-Jaccobi(A=b,xO);disp(-Jacobi二)；XJdisp(*迭代欠Jnjdisp(JGauss-SeidelSfl.:fi:xgzng-Gauss_Seidel(A,b,xO);disp(1Gauss-SeidelA，:；、：*)；dispC迭代次数冷;ng运行结果：精确解x=-0.0820-1.80331.1311迭代初始

13、值xO=000Jacobi迭代过程：X=-0.40000.50000.6667err=0.6667x=0.1000-1.03330.4533err=1.5333-0.0820-1.80331.1311err=Jacobi最终迭代结果：xj=-0.0820-1.80331.1311迭代次数nj=281Gauss-Seidel迭代过程：x=-0.40000.30000.5067err=0.5067x=-0.0360-0.53130.8012err=0.8313-0.0256-1.11510.9589err=05838x=-0.0820-1.80331.1311err=9.4021e-006Gaus

14、s-Seidel最终迭代结果:xg=-0.0820-1.80331.1311迭代次数ng=8.取不同的初值用脂”如送代以及弦截法求方程、3十2好十10乂-100-0的实根，列表或者画图说明收敛速度.悔：Newton迭代法:Newtoniter.mfunctionx.ite【,fvalue-NewK)nitei(fdfxO.eps.maxiter)%牛顿法X得到的近似解%iter迭代次数%fvralue函数在x处的值%f.df彼求的非线性方程及导函数%x0初始值%eps允许误差限%maxicer最人迭代次数fvalue-subs(fjLO);dfvalue-subs(dfxO);foriter-

15、1niaxiterx-KO-fValuedfvahieenabs(x-xO)xO-xifvakiLsubs(f.xO)dfalue-subs(df,xO);if(err：0:1初始值%eps允许误差限%maxit:er最人迭代次数falueO-subs(f,xO);falue-subs(f,xl);foriter-1niaxiterx-xl-fvalue*(xl-xO)(fvahie-fvahieO)enabs(x-xl)xO-xl:xl-x:fvalueO-subs(fjiO);fvralue-subs(f,xl)if(err:n2-Newlomter(f,dfxn2_0,eps,maxit

16、er)dispC1,；：-eps)inf-subs(ft(ab)/2);0)x-inf;n-n+l;retumendif(uiPfa:4,iter4-resecm(f.a32b,eps)运行结果：o.pi区间的根x1=1.8807;迭代次数iterl=20pi.2*pi区间的根x2=4.6941;迭代次数iter2=202*pi.3*pi区间的根x3=7.8548;迭代次数iter3=203*pi.4*pi区间的根x4=10.9955;迭代次数iter4=2010.考虑函数f(x)=sin(?rx)0,1.用等距节点作/，(金)的New如i插值，画出插值多项式以及/(可的图像，观察收敛性.程序

17、：Newton插值：Newtominter,mfunctionf-Newtominter(x,y.xO)务牛顿插值次插值节点务y为对应的函数值务函数返回Newton插值多项式在：_0点的值fsyinst;if(length(x)-length(y)n-longth(x);c(ln)-0.0;elsedisp(5和y为维数不术等!return:endf-yd);yi-o;1-1；forfbiyiU)一(y(j)-y(O);C(1)-yl(i+D;1-f-f+c(i)*1;sin)plify(f);y-yi;if-n-1)if(nargin3)口、如果3If-subsll2xO)；elset -

18、collect (f);f-vpa(f, 6);$如果2个参数则直接给出排值参项式$将插值多项式展开用等距节点做f(x)的Newton插值:testl0.miil-5;n2-10;n3-15;0-(0:001:1;yO-sin(pi-x0);xl-linspace(0,1,nl)；务等距,；,：215yl-sin(pi*xl);fO1-Nexxlominter(xl.yl.xO);:2-lmspace(0,l,n2);勺容汕巧:一10y2-sm(pi-x2);f02-Newtominter(x2：y2、xO);x3-linspace(0,l.n3)；务等蹈巧声，-15y3-sin(pi*:3)

19、;f03-Neuiominter(x3,y3.xO);plot(xO.y0：1-r)原图holdonplot(x0:f01,-g)%5卜plot(x0:fl)2:个节点plot(x0:fl)3:l-b)%15、：legend-以图J5，Nexrton：10kfj,占Newton抬油务项式15个节点NexNon插值多项式r)运行结果：取不同的节点做牛顿插值。得到结果图像如下：120203060708可以看出原图与插值多项式的图像近似重合，说明插值效果较好。1L对函数/(尸匚?7.xE-5.5.取不同的节点ISzb用等距节点作lagrangf插值,观察现象.Lagrange插值：Lagrango.

20、mfunctionf.fD一Lagrange(x,y,xO)%Lagrange插气，y函数值mO为要计算的点。$函数返BL_n(:-:)农达式f和L_n(x0)的值fOosyins二if(length(x)-length(y)n-length(x);elsedisp(巳和y的维数不相等！)return:end$检错f-0.0;for(i-l:n)i-y(0;for(j-1:M)1-P(t-x(j)/(x(i)-x(j);end;务计 /Lagrange 教for(j-i+l:n)1-(X-end;电如果3个参数则计算插值点的函数值毛如果2个参数则将插值多项式展开$将插值多项式的系数化成6位精度

21、的小数f-f+1;simplify；if(i-n)if(nargm-3)fO-subs(fJdzxO)：elsef-collect(f);f-vpa(f,6);务计算Lagrange;甫情函数务化简用等距节点做Lagrange插值：testU.mclearall,clc;nl-5;n2-10;n3-15;x0-5:002:5;yOT(lr0a2)；xl-linspace(-5,5.nl)；务等跖工.肖15/(1+xl-A2);flzfDl-Lagrange(xlzyl：xO);:2-lmspace(-5,5.n2);节点.盯点数10y2-1-/(l+:2人2);f2zf02-Lagrange(

22、x2,y2,:0);x3Tinspace(-5,5,n3);告等距1H15y3-卜/(lr3A2);f3.f03-Lagrange(:3,y3.xO)：plot(xO.y021-r1)务原图holdonplot(x0-f012个节点plot(xO：fD2,:g)o1Oplot(x0：f03-r-kf)%15个节点xlabelx*);ylabel(ff(x),L(x)1);legend(，醺图4ZLagrange插T/T10个节：Lagrange插值多项式,；15|)Lagrange插L务项式I)运行结果：选取了5.10,15个节点做Lagrange插值，得到原图与插值多项式的图像如下:当选取右

23、点数较多时，选取15个点时出现Runge现象，12.令仅=介84?a),考虑积分小区间分为50,100,200,500,1000等，分别用复合梯形以及复合S/如so胪分公式计算积分值，将数值积分的结果与精确值比较，列表说明误差的收敛性.复合梯形求积：Comtrap.mfunctions-Comtiap(f,a.bji)务复合梯形求积S积分结果%只分函数%a,b积分区间%n、二工h-(b-a)/n;index-(a-rh):h:(b-h);sl-sum(subs(f,mdex);s-(h2)*(subs(f,a/2*sl+subs(fzb);复合Simpson求积：functions-Smips

24、on(f,a.bji)$复合Simpson求积s积分结果积分函数%a,b积分区间%11区间个数h-(b-a)(2*n);indexl-(a+h):(2*h):(b-h);mdex2-(a+21i):(2*h):(b-2*h);sl-sum(subs(f.indexl);s2-sum(subs(f,index2);s-(h/3)*(subs(fza)*4*sl+2*s2*subs(f.b);计算f(x)积分：test12.mclearall;clc;svmsxfexp(3:)*cos(pi*x);a-0;b-2*pi;&p(*精确积分值I-vpa(int(f.x,a,b),10)nl-50;n2

25、-100m3-200;n4-500:n5-1000:disp区间为50.100.200.500.1000的复合梯形积分值，);Tl-pa(Comtrap(f,a.b,iil).10)e*tl-abs(I-Tl)T2-x,pa(Comuap(f,a.bji2)J0)et2-abs(I-T2)T3-X-pa(Comtiap(f.a:b.n3),10)et3-abs(I-T3)T4-xrpa(Comtiap(f,a,b,n4),10)ec4-abs(I-T4)T5-x-pa(Comtiap(f,a-b,n5).10)et5-abs(I-T5)dispj区诃为50.100.200,500,1000“U

26、合Smipson:J;sl-pa(Snnpson(f,a,b,nl).10)esl-abs(I-sl)s2-vp“Smipson(fab,n2),10)es2-abs(I-s2)s3-,pa(Sinipson(f.a.b.ii3).10)es3-abs(I-s3)s4-vpa(Simpson(f.a：bji4),10)es4-abs(I-s4)s5-pa(Simpson(f,a,b,n5)40)es5-abs(I-s5)运行结果：精确积分值I=35232483.36复合梯形与复合Simpson求积与误差区间个数n复合梯形求积T误差eT复合Simpson求积误差eS5035125341.1910

27、7142.1735231407.871075.4910035204891.227592.1635232416.2467.1220035225534.986948.3835232479.174.1950035231369.371113.9935232483.250.11100035232204.78278.5835232483.360.013.分别用2点，3点以及5点的du岛型积分公式计算如下定积分:,z sin x f(2) 1dxin 丁GaussLegendre求积：GaussLegendre*mfunctions-GaussLegendre(f.a.b.n)%Gauss-Legendre

28、公式求积分.只给2-5个求积结点%f被积函数%b,a积分上下限%n求积结点个数%s积分结果tarb-a)/2;tb-(a+b)/2;switchnsra*(subs(sym(f).findsym(sym).5773503+1匕)+subs(sym(f).findsym(sym(f)tdO5773503+tb);s-ta*(05555556丹subs(sym(f).findsym(sym(f),tawO-7745967+tb)+05555556*subs(sym(f),findsym(sym(f),-ta*O-7745967*tb)十一08888889*subs(sym(f).findsym(s

29、yni(f),tb);case3,s-ta*(0-3478548subs(sym(f).findsym(sym(f),ki08611363+tb)+03478548subs(sym(f),findsym(sym(f),-ta*08611363+tb)+06521452”subs(sym(f),findsym(sym(f),ta本03399610十tb)十。0-6521452*subs(sym(f),ifindsym(sym(f),ta*O-3399810-tb);case4is-ta*(O-2369269.subs(sym(f)、findsym(sym(f)、ta本09061798+tb)+0

30、2369269字subs(sym(f),findsym(sym(f),-ta*O9061798+tb)+04786287*subs(sym(f).findsym(sym(f),ta*0-5384693-tb)十,04786287*subs(syrm(f).findsym(sym(f).-ta*O-5384693-tb)十,05688889.subs(sym(f)fndsym(sym),tb);endGaussChebyshev求积：GaussChebyshev.mfunctions-GaussChebyshev(f,n)%GaussChebyshev求积s积分值%什只分函数务求积结点个数k1k

31、-0nx-cos(2-率M)*pi/(2*(n+1);a-piz(n+l);s-a*sum(subs(fjc);计算f(x)积分：test13.mcleaiall.clc;syms:fl-xA2;p-sqrt(bx-A2);f2-sin(x)/x:nl-2;n2-3;n3-5;a-O.b-pi/2;disp(,);H-vpa(int(f1p.xL1),8)disp(值，)；I2-vpa(int(f2xa.b).8)disp(；,个235Guass冷;I2-vpa(GaussChebyshev(fl.nl-1).8)13rpMGaussChebyshev(fl,n2-l),8)I5-vpa(Ga

32、ussChebyshev(f1,n3-l),8)disp(第2.3.5Qu*型积分)；GL2-xpa(GaussLegendre(f2,a,bznl-l),8)GL3-vpa(GaussLegendre(f2：a.b.n2-l),8)GL5-vpa(GaussLegendre(f22a：bzn3-l),8)运行结果：(1)实际积分值：11=1.57079632.3.5点GaussChebyshev型积分12=1.570796313=1.570796315=1.5707963(2)实际积分值12=1.37076222.3.5点GaussLegendre型积分GL2=1.370419GL3=1.3

33、707635GL5=1.370762214.考虑微分方程初值问题:(务44p(纫),U(0)=2.分别用奶J、改进的民后法，忌如以俅解该方程.分别取步长为。i.o.ouoooi,计算到父.画图说明结果.酢：Euler法：Euler-mfunctionT-Euler(f.xO,xn.yO：h)%Euler法%x?y微分方程的解%微分方程%xO:yO初始值%xn端点%h步长n-(xn-xO)/h;&区叵x-zeros(lzn+1);y*zeros(Lnl);x(l尸xO:yd)-yo;fork-l:ny(k+lr(k)十h*feval(f：x(k)2y(k);x(kHmk)十h：endT-XzY1

34、;改进Euler法：TranEuler.mfunctionT-TraiiEuler(f,xOz:n-yO,h)务改进Euler法%x：y微分方程的解%微分方程%xdyO初为仔i%xn端点%h步长n-(xii-xO)/h;x-zeios(Ln-t-l);卜zaos(l,n+1);x(l)-xO;yd)-yo;foxk-l:nx(kH尸x(k)+h;zOy(k)h*feval(fx(k),y(k);处+1尸淋)+h/2*(feval(f.x(k)：y(k)+feFKf.x(k+l)：zO);endT-恒卬；经典Runge-Kutta法:R_K4mfunctionR-R_K4(f,xO,xn,yO：h)务经典，4阶RungeKutta法%力微分方程的解%微分方程%:0:5-0初始值%xn端点$h步长n-(xn-xO)/h:；区|1.y-zeros(l：n+1);yd)-yo;x-xO:hxn:fork-l:nkl-feval(frX(k)ty(k);k2-feva