大连理工大学-矩阵大作业

1.设Sw=4/'其精确值为;C-，1VTT)，工)编制技从大到小的顺序5盘=——+——十…+T—，计算Sw的通用2』—1炉一I 闻士—1程序也)编制按从小到大的顺序S1V= ,+...4 计算Su—1(N—一1 /一1的通用程序口｝按两种顺序分别计算51四5寸，5W,并指出有效位数(编制程序时用单精度)也)通过本上机题，你明白了什么(1)从大到小的顺序的计算程序：functiony=snd(n)formatlongy=0;ifn<2disp('请输入大于1的数！')elses=0;i=2;whilei<=ns=single(s+(1/(/2T)));i=i+1;endy=s;end(2)从小到大的顺序的计算程序：functiony=snx(n)formatlongy=0;ifn<2disp('请输入大于1的数！')elses=0;i=n;while1s=single(s+(1/(/2T)));i=i-1;ifi==1breakendendy=s;end(3)按两种顺序分别计算5” 并指出有效位数(编制程序时用单精度)①S的计算结果：102»y=snd(LOO)****就学M大至I］小*****肘¥=0.7100195»y=sns(L001++++*兴从.小到大率率率率今率y-0.7400495②S的计算结果：104»y=snd(］01)******从大到小ww**ww7=0.7^98521»y=sn3(104J卓卓卓卓卓卓MM、至l［大息章才*..V=0."900。③S的计算结果：106»y=snd(10,61常京算**±加大至同\京*常常才京y=0.7498621»^snx(10"6)A*jc**K小小至I］大军芯常常才享y=0.7499990计算时的有效位数为七位数。L秦丸貂算法.已知h次多项式六幻/上口酒,用秦九韶算法编写通用的程序计算函数i=n在》点的值T并廿算〃£)=7/+3/- 11在23点的值.［提示：编写程序时,输入系数向量和点看,输出结果，多项式的次整可以通过向殳的长度来判断I①秦九昭算法计算程序：functiony=qjz(a,x)j=3；i=size(a,2);switchiy=a(1);y=a(1)*x+a(2);otherwisep=a(1)*x+a(2);whilej<=ip=p*x+a(j);j=j+1;endy=p;end②计算,/(『)I3〃 ％ 11在点23的值。»2=［乙为一5,11］7 3 -5 11»y=q.j2(a.23)富的值为二2386652当x=23时fQ)=86652。3.分别用Gauss消元法和列主元消•去■法编程求解方程组M=l3,其中-31-13 0 0 0 -10 00 0-一1335 -9 0 -11(J00 0 27o-9 -ioofld0 0-230t)-107D-30a00-fl0[A.fr]= 0 0U—3。57 -7t)-50TO.0 0 0D-7 47-300 0L20 0 (9 0 0 -30410 0-70 0U0 -5 027—2 70 0 0-9000-2如L0①Gauss法计算程序和结果：程序：A(1,:)=[31,-13,0,0,0,-10,0,0,0];A(2,:)=[-13,35,-9,0,-11,0,0,0,0];A(3,:)=[0,-9,31,-10,0,0,0,0,0];A(4,:)=[0,0,-10,79,-30,0,0,0,-9];A(5,:)=[0,0,0,-30,57,-7,0,-5,0];A(6,:)=[0,0,0,0,-7,47,-30,0,0];A(7,:)=[0,0,0,0,0,-30,41,0,0];A(8,:)=[0,0,0,0,-5,0,0,27,-2];A(9,:)=[0,0,0,-9,0,0,0,-2,29];B=[-15;27;-23;0;-20;12;-7;7;10];[a,b]=gauss(A,B);j=size(a,2);whilej>=1k=j+1;s=b(j);whilek<=9s=s-x(k)*a(j,k);k=k+1;endx(j)=s/a(j,j);j=j-1;enddisp(x)function[x,y]=gauss(a,b)num_i=size(a,1);j=1;whilej<=(num_i-1)i=j+1;whilei<=num_ir=a(i,j)/a(j,j);a(i,:)=a(i,:)-r*a(j,:);b(i,:)=b(i,:)-r*b(j,:);i=i+1;endj=j+1；endx二a;y=b；运行的结果为：%=(-0.2890.345-0.713-0.221-0.4300.154-0.0580.2010.290万。②列主元消去法计算程序和结果：A(1,:)=[31,-13,0,0,0,-10,0,0,0];A(2,:)=[-13,35,-9,0,-11,0,0,0,0];A(3,:)=[0,-9,31,-10,0,0,0,0,0];A(4,:)=[0,0,-10,79,-30,0,0,0,-9];A(5,:)=[0,0,0,-30,57,-7,0,-5,0];A(6,:)=[0,0,0,0,-7,47,-30,0,0];A(7,:)=[0,0,0,0,0,-30,41,0,0];A(8,:)=[0,0,0,0,-5,0,0,27,-2];A(9,:)=[0,0,0,-9,0,0,0,-2,29];B=[-15;27;-23;0;-20;12;-7;7;10];[a,b]=lzy(A,B);j=size(a,2);whilej>=1k=j+1;s=b(j);whilek<=9s=s-x(k)*a(j,k);k=k+1;endx(j)=s/a(j,j);j=j-1;enddisp(x)function[a1,b1]=lzy(a,b)[num_i,num_j]=size(a);ab=zeros(num_i,num_j+1);fork=1:num_jab(:,k)=a(:,k);endab(:,num_j+1)=b(:,1);j=1;whilej<num_j[max,max_i]=searmax(j,j,ab);i_nu=ab(j,:);ab(j,:)=ab(max_i,:);ab(max_i,:)=i_nu;m=j+1;whilem<=num_iforn=j:num_j+1ab(m,n)=ab(m,n)-(ab(m,j)/max)*ab(j,n);endm=m+1;endj=j+1;enda1=zeros(num_i,num_j);fork=1:num_ia1(:,k)=ab(:,k);endb1=ab(:,num_i+1);function[b,c]=searmax(i,j,a)num_i=size(a,1);k=i;m=abs(a(k,j));c=k;whilek<num_ik=k+1;ifm>=abs(a(k,j))continueelsem=abs(a(k,j));c=k;endendb=a(c,j);运行的结果为：%=J0.2370.477—0.767—0.078—0.3080.146—0.1710.2850.345万14.编程求解题3中矩阵A的LU分解及列主兀的LU分解(求出LrU和P),并用LU分解的方法求A的逆矩阵及A的行列式(1)LU分解的计算程序及结果：function[l,u]=lufz(a)[num_i,num_j]=size(a);l=eye(num_i);u=eye(num_i);forj=l:num_ju(l,j)=a(l,j);1(j,l)=a(j,l)/u(l,1);endi=2;whilei<=num_ij=i；whilej<num_is=0;fork=l:i-1s=s+l(i,k)*u(k,j);endu(i,j)=a(i,j)-s;s=0;fork=l:i-1s=s+l(j+1,k)*u(k,i);end(j+1,i)=(a(j+l,i)-s)/u(i,i);j=j+l；ends=0;fork=l:i-1s=s+l(i,k)*u(k,num_i);endu(i,num_i)=a(i,num_i)-s;i=i+l;end输入要求解的矩阵后运行的结果为：(1-0.419 10-0.305100-0.3541L=000-0.39810000-0.157100000-0.65510000-0.112-0.018-0.0251<000-0.119-0.083-0.014-0.020-0.092b

(31-13000-10000)29.548-90-11-4.19400028.259-10-3.350-1.27700075.461-31.186-0.45200-9U=44.602-7.1800-5-3.57845.8732-30-0.785-0.56221.381-0.513-0.36726.413-2.420V27,390)(2)带列主元的LU分解计算程序及结果由于Matlab中自带带列主元的LU分解函数，故计算程序如下:A(1,:)=[31,-13,0,0,0,-10,0,0,0];A(2,:)=[-13,35,-9,0,-11,0,0,0,0];A(3,:)=[0,-9,31,-10,0,0,0,0,0];A(4,:)=[0,0,-10,79,-30,0,0,0,-9];A(5,:)=[0,0,0,-30,57,-7,0,-5,0];A(6,:)=[0,0,0,0,-7,47,-30,0,0];A(7,:)=[0,0,0,0,0,-30,41,0,0];A(8,:)=[0,0,0,0,-5,0,0,27,-2];A(9,:)=[0,0,0,-9,0,0,0,-2,29];[L,U,P]=lu(A);运行结果如下：’1)-0.41910-0.305100-0.3541L=000-0.39810000-0.157100000-0.65510000-0.112-0.018-0.0251V000-0.119-0.083-0.014-0.020-0.0921)f31 -13000-100029.548-90-11-4.1940028.259-10-3.350-1.2770075.461-31.186-0.45200U=44.602-7.1800-545.8732-30-0.78521.381-0.51326.4130I00-9-3.578-0.562-0.367-2.42027,390)0000P=00000000100000000100010001000000000000000000000010010000001000000为单位阵。00001001)(3)求A的逆矩阵的计算程序及结果:A(1,:)=[31,-13,0,0,0,-10,0,0,0];A(2,:)=[-13,35,-9,0,-11,0,0,0,0];A(3,:)=[0,-9,31,-10,0,0,0,0,0];A(4,:)=[0,0,-10,79,-30,0,0,0,-9];A(5,:)=[0,0,0,-30,57,-7,0,-5,0];A(6,:)=[0,0,0,0,-7,47,-30,0,0];A(7,:)=[0,0,0,0,0,-30,41,0,0];A(8,:)=[0,0,0,0,-5,0,0,27,-2];A(9,:)=[0,0,0,-9,0,0,0,-2,29];[num_i,num_j]=size(A);A_n=eye(num_i);E=eye(num_i);[L,U]=lufz(A);fori=1:num_iy=hd(1,L,E(:,i));A_n(:,i)=hd(0,U,y);enddisp(A_n)function[l,u]=lufz(a)[num_i,num_j]=size(a);l=eye(num_i);u=eye(num_i);forj=1:num_ju(1,j)=a(1,j);l(j,1)=a(j,1)/u(1,1);endi=2;whilei<=num_ij=i;whilej<num_is=0;fork=1:i-1s=s+l(i,k)*u(k,j);endu(i,j)=a(i,j)-s;s=0;fork=1:i-1s=s+l(j+1,k)*u(k,i);endl(j+1,i)=(a(j+1,i)-s)/u(i,i);j=j+1;ends=0;fork=1:i-1s=s+l(i,k)*u(k,num_i);endu(i,num_i)=a(i,num_i)-s;i=i+1;endfunctionx=hd(f,a,b)[num_i,num_j]=size(a);x=zeros(num_i,1);switchfcase0i=num_i-1;x(num_i)=b(num_i)/a(num_i,num_j);whilei>=1s=0;fork=i+1:num_is=s+a(i,k)*x(k);endx(i)=(b(i)-s)/a(i,i);i=i-1;endcase1x(1)=b(1)/a(1,1);j=2；whilej<=num_is=0;fork=1:j-1s=s+a(j,k)*x(k);endx(j)=(b(j)-s)/a(j,j);j=j+1;endotherwisedisp('error！请输入正确的参数！')end运行结果：/0.03900.01600.00580.00360.00730.01760.01290.00140.0012、0.01590.03780.01310.00650.01210.00970.00710.00240.00220.00490.01160.03810.00770.00690.00390.00280.00150.00250.00080.00200.00640.01810.01050.00330.00240.00240.0058A-i=0.00050.00110.00360.01010.02430.00700.00510.00480.00350.00010.00030.00100.00280.00680.04190.03060.00130.001000.00020.00070.00210.00500.03060.04680.00100.00070.00010.00020.00080.00230.00480.00140.00100.03820.0033[0.00030.00060.00200.00580.00360.00110.00080.00340.0365)(4)求A的行列式的计算程序及结果:A(1,:)=[31,-13,0,0,0,-10,0,0,0];A(2,:)=[-13,35,-9,0,-11,0,0,0,0];A(3,:)=[0,-9,31,-10,0,0,0,0,0];A(4,:)=[0,0,-10,79,-30,0,0,0,-9];A(5,:)=[0,0,0,-30,57,-7,0,-5,0];A(6,:)=[0,0,0,0,-7,47,-30,0,0];A(7,:)=[0,0,0,0,0,-30,41,0,0];A(8,:)=[0,0,0,0,-5,0,0,27,-2];A(9,:)=[0,0,0,-9,0,0,0,-2,29];[num_i,num_j]=size(A);[L,U]=lufz(A);s=1;fori=1:num_is=s*U(i,i);enddisp(['矩阵的行列式值为：'，num2str(s)])运行程序后结果与调用matlab中det()函数结果如下:矩阵的行列式值为：€1517405101e50.01»d&t⑷an吕=61817405101650(1)求cholesky分解程序及结果:functionl=chole(a)[num_i,num_j]=size(a);l=zeros(num_i);l(1,1)=sqrt(a(1,1));fork=2:num_il(k,1)=a(k,1)/l(1,1);endj=2;whilej<=num_js1=0;fork=1:j-1s1=s1+l(j,k)"2;endl(j,j)=sqrt(a(j,j)-s1);i=j+1;whilei<=numis2=0;fork=1:j-1s2=s2+l(i,k)*l(j,k);endl(i,j)=(a(i,j)-s2)/l(j,j);i=i+1;endj=j+1;end运行程序后的结果为:T.4142000)0.70712.121300L=—0.70711.17852.84800[0.70713.0641—0.74130.7494,(2)求解方程组程序及结果：a=[2,1,-1,1;1,5,2,7;-1,2,10,1;1,7,1,11];b=[13;-9;6;0];L=chole(a);y=hd(1,L,b);x=hd(0,L',y);disp(x)运行后的结果为：x=(26.4146-65.390212.512238.0732)氏用追超法编制程序求解方程组加寸，其中~12I)if3-2102A=02 531000-16程序和结果如下：A=[4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6];B=[6;2;10;5];[L,U]=sjfj(A);%y=hd(1,L,B);x=hd(0,U,y);disp('方程组的解为：')disp(x)function[l,u]=sjfj(a)num_i=size(a,1);i=2;whilei<=num_ia(i,i-1)=a(i,i-1)/a(i-1,i-1);a(i,i)=a(i,i)-a(i,i-1)*a(i-1,i);i=i+1;endl=tril(a);u=triu(a);forj=1:num_i;l(j,j)=1;end运行程序后结果为：X=G11LT7,已知「II。U'—13-1/21/2-223/2L/2-22-1/25/2编程求解矩阵R的QR分解；编程求解矩阵A的QR分解：(1)QR分解计算程序：function[Q,R]=hqr(A)[n,n]=size(A);Q=eye(n);fori=1:(n-1)E=eye(n-i+1);e1=E(:,1);a=zeros(1,n-i+1)';forj=1:(n-i+1)a(j)=A(j+i-1,i);endav=sqrt(a，*a);w=a-av*e1;h=E-(2/(w'*w))*(w*w');q=eye(n);fork=i:nforj=i:nq(k,j)=h(k-i+1,j-i+1);endendA=q*A;Q=q*Q;endR=A;Q=Q';(2)输入矩阵A后的计算结果:'0.31620.70710.2582-0.5774、Q=-0.31620.7071—0.25820.5774—0.632500.77460[-0,63250—0.5164-0.5774,'3.1623-3.1623—0.4743—2.0555、R=02.8284—0.35360.3536001.5492—1.0328、000-1.1547)8.分别应用Jacobi迭代法和Gaus>Seidel迭代法求解如下方程组Ilij-工2+6=7<ir：i+租值+均=—211-2工1+g+HEJacobi迭代法求解程序与结果：a=0;b=0;c=0;while1a0=0.25*(7+b-c);b0=(1/8)*(-21-4*a-c);c0=(1/5)*(15+2*a-b);A=[abs(a-a0),abs(b-b0),abs(c-c0)];f=max(A);iff<=0.001x=[a0,b0,c0];breakelsea=a0;b=b0;c=c0;endenddisp('方程组的解为x')disp(x)运行后的结果为：％=6.060—3.1113.647万(2)Gauss-Seidel迭代法求解的程序与结果：a=0;b=0;c=0;while1a0=a;b0=b;c0=c;a=(1/4)*(7+b-c);b=(1/8)*(-21-4*a-c);c=(1/5)*(15+2*a-b);A=[abs(a-a0),abs(b-b0),abs(c-c0)];f=max(A);iff<=0.001x=[a,b,c];breakendenddisp('方程组的解为：')disp(x)运行后的结果为：％=6.061—3.1113.646万.分别应用Newt前迭代法和割线法计算(1)非线性方程加一E一1=。在[14上的一个根t(2) 由1£=0在卜4周上的一个根口(1)计算2x3-5x-1=0在1,21上的根的程序：symsxf=2*x"3-5*x-1;df=diff(f,x);g=x-(f/df);x0=1;while1x1=subs(g,x0);dx=abs(x1-x0);ifdx<=0.001disp(['Newton法的解为：x=',num2str(x1)])breakelsex0=x1;endendx0=1;x1=1.1;while1x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)))*(x1-x0);dx=abs(x2-x1);ifdx<=0.001disp(['割线法的解为：x=',num2str(x2)])breakelsex0=x1;x1=x2;endend运行后结果为：Newton法x=1.673,割线法x=1.673。(2)计算ssinx=0在L4,-31上的根的程序：symsxf=exp(x)*sin(x);df=diff(f,x);g=x-(f/df);x0=-2;while1x1=subs(g,x0);dx=abs(x1-x0);ifdx<=0.001disp(['Newton法的解为：x=',num2str(x1)])breakelsex0=x1;endendx0=-2;x1=-2.1;while1x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)))*(x1-x0);dx=abs(x2-x1);ifdx<=0.001disp(['割线法的解为：x=',num2str(x2)])breakelsex0=x1;x1=x2;endend运行后结果为：Newton法x=—3.1416，割线法x=—3.1416。.采用二分法计算非线性方程H侬金一？=Q,叠我区间为symsxf=x*cos(x)-2;x0=-4;x2=-2;while1x1=0.5*(x0+x2);f0=subs(f,x0);f1=subs(f,x1);c=f0*f1;ifc<0x2=x1;elsex0=x1;endl=abs(x2-x0);ifl<=0.001

disp(['二分法的解为x=',num2str(x0)])breakendend运行后的结果为x=-2.499。.已知函数在[-5,5]上分别取2,L支为单位长度的等距节点作为插值节点，用匕即的处方法插值，并把原函数图与插值函数图比较，观察插值效果(1)Lagrange插值程序：functionyv=lagran(xd,yd,xv)symsxnum=length(xd);%计算节点的个数f=0;fork=1:num%创建lagrange插值函数index=setdiff(1:num,k);f=f+yd(k)*prod((x-xd(index))./(xd(k)-xd(index)));endyv=subs(f,xv);%计算xv点插值函数值(2)绘制函数图程序：symsxf=1./(1+x「2);xd2=-5:2:5;yd2=subs(f,xd2);x0=-5:0.1:5;y=subs(f,x0);plot(x0,y)holdony2=lagran(xd2,yd2,x0);plot(x0,y2)(3)(3)运行程序后得到的步长分别为2,1,1 …i，…、E，-的插值函数图与原图形的比较图如下:1Z用三次样条插值上题中的插值节点，并画图比较插值效果,（提示：原函数在两个端点-5』5的导数值可作为边界条件）(1)三次样条插值程序：functionyv=csi(xd,yd,xv,h,mf,ml)num_xd=length(xd);n=num_xd-1;a_m=1:num_xd-2;fori=1:numxd-2a_m(i)=2;enda_np=1:num_xd-3;fori=1:num_xd-3a_np(i)=0.5;endA=diag(a_m)+diag(a_np,-1)+diag(a_np,1);g=1:n-1;fork=1:n-1g(k)=3*((0.5/h)*(yd(k+2)-yd(k+1))+(0.5/h)*(yd(k+1)-yd(k)));endb=1:n-1;b(1)=g(1)-0.5*mf;b(n-1)=g(n-1)-0.5*ml;fork=2:n-2b(k)=g(k);endb=b.';M=A\b;m=1:n+1;m(1)=mf;m(n+1)=ml;fori=2:nm(i)=M(i-1);endnum_xv=length(xv);yv=1:num_xv;fori=1:num_xvforj=1:num_xdifxv(i)<xd(j)k1=j-1;k2=k1+1;breakendendyv1=(h+2*(xv(i)-xd(k1)))*(xv(i)-xd(k2))"2*(yd(k1)/h"3);yv2=(h-2*(xv(i)-xd(k2)))*(xv(i)-xd(k1))"2*(yd(k2)/h"3);yv3=(xv(i)-xd(k1))*(xv(i)-xd(k2))"2*(m(k1)/h"2);yv4=(xv(i)-xd(k2))*(xv(i)-xd(k1))"2*(m(k2)/h"2);yv(i)=yv1+yv2+yv3+yv4;end(2)绘图程序(改变程序中的h值即可改变步长)：clearallclch=0.5;xd=-5:h:5;yd=1./(1+xd「2);xv=-5:0.1:5;y=1./(1+xv「2);mf=0.014793;ml=-0.014793;yv=csi(xd,yd,xv,h,mf,ml);plot(xv,y)holdonplot(xv,yv)(3)步长为2，1，1时原函数图与插值函数图的比较:13.已知=X13.已知=X2sinxi分别用复合梯形法和复合Simp60n公式计算积分Cja”3区间分为20・40,80,加0个小区间，并计算其精确值，比较计算精度情［兄口(1)复化梯形公式计算积分程序：functiont=trape(f,a,b,n)h=(b-a)/n;index=(a+h):h:(b-h);t1=sum(subs(f,index));t=((b-a)/(2*n))*(subs(f,a)+2*t1+subs(f,b));(2)复化Simpson公式计算程序：functions=simpson(f,a,b,n)h=(b-a)/n;index1=(a+0.5*h):h:(b-0.5*h);index2=(a+h):h:(b-h);s1=sum(subs(f,index1));s2=sum(subs(f,index2));s=((b-a)/(6*n))*(subs(f,a)+4*s1+2*s2+subs(f,b));（3）区间数不同时的计算结果：（其中』2x2sinxdx=0）-2**区间个数为**区间个数为2。时科卓**夏化桃形计算转果***-8,8918#-017**复化Simpson计算结果***2.0724e-016**区间个数为40时*****复化梯肥计算结果***1.3323e-016**豆化Simpart计箕结果*.**l.«2a：3e-016区间个料为80时*比***夏优梯形计算结果***S.SSlBe-017**复化5i呻鸟皿计茸结果不****区间个翻为2□。时*****原北稀杷计算鳍果-8.8818e-018**官化SiMpwon计耳结果*”6.2172e-01714.用2点，