大连理工大学矩阵上机.docx

1、大连理工大学矩阵上机大连理工大学DALIAN UNIVERSITY OF TECHNOLOGY 研究生作业课 程 名 称： 矩阵与数值分析 研究生姓名 ： 学 号： 作 业 成 绩： 任课教师签名 能源与动力学院 上交作业时间：2015年12月25日1. 设，其精确值为. （1） 编制按从大到小顺序，计算的通用程序；（2） 编制按从小到大顺序，计算的通用程序；（3） 按两种顺序分别计算, ,，并指出有效位数（编制时的单精度）；（4） 通过本次上机题，你明白了什么？MATLAB程序function s1,s2=s(N)format long gs1=0;s2=0;for j=2:Ns1=1.0/

2、(j*j-1)+s1;endfor j=N:-1:2s2=1.0/(j*j-1)+s2;endend计算结果 s1,s2=s(1.0e2)s1 =0.740049504950495s2 =0.740049504950495 s1,s2=s(1.0e4)s1 =0.749900004999506s2 =0.7499000049995 s1,s2=s(1.0e6)s1 =0.749999000000522s2 =0.7499990000005结果分析计算机做加减法时，运算次序会影响结果，计算和时应先安排绝对值小的数参与运算，这样能取得较高的精度。2. 秦九韶算法。已知n 次多项式，用秦九韶算法编写

3、通用的程序计算函数在点的值，并计算在23点的值。(提示：编写程序时，输入系数向量和点，输出结果，多项式的次数可以通过向量的长度来判断).MATLAB程序A = input(,);n = input();len = length(A);val = zeros(len);val(1) = A(1);%printf(len = %c, len)for i = 2: 1: len %disp(val(i-1) %disp(n) val(i) = val(i-1)* n + val(i);end%printf(?%f, val(len)val(len)计算结果请输入系数，由高次幂开始7 3 -5 11请

4、输入计算变量的值23ans = 851693. 分别用Gauss 消元法和列主元消去法编程求解方程组Ax=b，其中。Gauss消元法function x=gauss(a,b)n=length(b);for k=1:n-1 if a(k,k)=0 fprintf(Error:the pivot element equal to zero!n,k); return; end index=k+1:n; m=-a(index,k)/a(k,k); a(index,index)=a(index,index)+m*a(k,index); b(index)=b(index)+m*b(k);endx = ze

5、ros(n,1);x(n) = b(n)/a(n,n); for i = n-1:-1:1 x(i) = ( b(i) - a(i,i+1:n)*x(i+1:n) )/a(i,i); end计算结果 a=31,-13,0,0,0,-10,0,0,0;-13,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,

6、29 b=-15;27;-23;0;-20;12;-7;7;10x=gau(a,b)x = -0.289233816015754 0.345435715779115 -0.712811731086879 -0.220608510570529 -0.430400432704022 0.154308739838311.0578* 0.201053894823681 0.290228661879745列主元消去法functionx=gauss(a,b)n=length(a);x=zeros(n,1);a=a b;for k=1:n-1max=k;for i=k+1:nif a(i,k)a(max,k

7、)max=i;endendtemp=a(k,k:n+1);a(k,k:n+1)=a(max,k:n+1);a(max,k:n+1)=temp;for i=k+1:na(i,k)=-a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)+a(i,k)*a(k,k+1:n+1);endendx(n,1)=a(n,n+1)/a(n,n);for i=n-1:-1:1sum=0;for j=i+1:nsum=sum+x(j,1)*a(i,j);endx(i,1)=(a(i,n+1)-sum)/a(i,i);end计算结果 a=31,-13,0,0,0,-10,0,0,0;-13

8、,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,29 b=-15;27;-23;0;-20;12;-7;7;10 x=gauss(a,b)x =-0.2892338160157540.345435715779115-0.712811731086879-0.220608510570529-0.430400

9、4327040220.154308739838311.0578*0.2010538948236810.2902286618797454. 编程求解题3中矩阵A的LU分解及列主元的LU分解（求出L,U和P），并用LU 分解的方法求A 的逆矩阵及A 的行列式.LU分解function L,U=mylu(a) n,n=size(a);for i=1:n L(i,i)=1;endU(1,1)=a(1,1)/L(1,1);if L(1,1)*U(1,1)=0 fprintf(Factorization impossible) else for j=2:n U(1,j)=a(1,j); L(j,1)=a(

10、j,1)/U(1,1); endend for i=2:n-1 sum=0; for k=1:i-1 sum=sum+L(i,k)*U(k,i); end U(i,i)=(a(i,i)-sum)/L(i,i); if L(i,i)*U(i,i)=0 fprintf(Factorization impossible) else for j=i+1:n h=0; s=0; for k=1:i-1 h=h+L(i,k)*U(k,j); s=s+L(j,k)*U(k,i); end U(i,j)=1/L(i,i)*(a(i,j)-h); L(j,i)=1/U(i,i)*(a(j,i)-s); end

11、endend sum=0; for k=1:n-1 sum=sum+L(n,k)*U(k,n); U(n,n)=(a(n,n)-sum)/L(n,n); endend计算结果 L,U=mylu(a)L = Columns 1 through 3 1 0 0 -0.419354838709677 1 0 0 -0.304585152838428 1 0 0 -0.353872899362565 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 4 through 6 0 0 0 0 0 0 0 0 0 1 0 0 -0.397554925856813 1 0 0 -0.

12、156943635949482 1 0 0 -0.653976718762685 0 -0.112102597106773 -0.0175453757582425 -0.119266477757044 -0.0833908821051642 -0.0142268147798013 Columns 7 through 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -0.0246185256433648 1 0 -0.0199621375629645 -0.0923164702590968 1U = Columns 1 through 3 31 -13 0

13、 0 29.5483870967742 -9 0 0 28.2587336244541 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 4 through 6 0 0 -10 0 -11 -4.19354838709677 -10 -3.35043668122271 -1.27729257641921 75.4612710063743 -31.185*5 -0.451999227351748 0 44.6019996774714 -7.17969451931716 0 0 45.8731926371318 0 0 0 0 0 0 0 0 0 Column

14、s 7 through 9 0 0 0 0 0 0 0 0 0 0 0 -9 0 -5 -3.57799433271131 -30 -0.784718179747412 -0.561543439982356 21.3806984371195 -0.513187420344639 -0.367236336322372 0 26.4130849214705 -2.41999576495225 0 0 27.3895043327769LU分解计算a的行列式和逆矩阵 det(U)ans = 6.1817e+13 inv(U)*inv(L)ans = 0.0390 0.0160 0.0058 0.003

15、6 0.0073 0.0176 0.0129 0.0014 0.0012 0.0159 0.0378 0.0131 0.0065 0.0121 0.0097 0.0071 0.0024 0.0022 0.0049 0.0116 0.0381 0.0077 0.0069 0.0039 0.0028 0.0015 0.0025 0.0008 0.0020 0.0064 0.0181 0.0105 0.0033 0.0024 0.0024 0.0058 0.0005 0.0011 0.0036 0.0101 0.0243 0.0070 0.0051 0.0048 0.0035 0.0001 0.00

16、03 0.0010 0.0028 0.0068 0.0419 0.0306 0.0013 0.0010 0.0001 0.0002 0.0007 0.0021 0.0050 0.0306 0.0468 0.0010 0.0007 0.0001 0.0002 0.0008 0.0023 0.0048 0.0014 0.0010 0.0382 0.0033 0.0003 0.0006 0.0020 0.0058 0.0036 0.0011 0.0008 0.0034 0.0365列主元LU分解function l,u,p=lielu(a) m,n=size(a);if m=n error() re

17、turnendif det(a)=0 error()endu=a;p=eye(m);l=eye(m);for i=1:m for j=i:m t(j)=u(j,i); for k=1:i-1 t(j)=t(j)-u(j,k)*u(k,i); end end a=i;b=abs(t(i); for j=i+1:m if b a=31,-13,0,0,0,-10,0,0,0;-13,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,

18、0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,29; l,u,p=lielu(a)l = 1.0000 0 0 0 0 0 0 0 0 -0.4194 1.0000 0 0 0 0 0 0 0 0 -0.3046 1.0000 0 0 0 0 0 0 0 0 -0.3539 1.0000 0 0 0 0 0 0 0 0 -0.3976 1.0000 0 0 0 0 0 0 0 0 -0.1569 1.0000 0 0 0 0 0 0 0 0 -0.6540 1.0000 0 0 0 0 0 0 -0.1121 -0.

19、0175 -0.0246 1.0000 0 0 0 0 -0.1193 -0.0834 -0.0142 -0.0200 -0.0923 1.0000u = 31.0000 -13.0000 0 0 0 -10.0000 0 0 0 0 29.5484 -9.0000 0 -11.0000 -4.1935 0 0 0 0 0 28.2587 -10.0000 -3.3504 -1.2773 0 0 0 0 0 0 75.4613 -31.1856 -0.4520 0 0 -9.0000 0 0 0 0 44.6020 -7.1797 0 -5.0000 -3.5780 0 0 0 0 0 45.

20、8732 -30.0000 -0.7847 -0.5615 0 0 0 0 0 0 21.3807 -0.5132 -0.3672 0 0 0 0 0 0 0 26.4131 -2.4200 0 0 0 0 0 0 0 0 27.3895p = 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 15. 编制程序求解矩阵A 的ch

21、olesky 分解，并用程序求解方程组Ax=b，其中， .MATLAB程序A= 2 1 -5 1 1 -5 0 7 0 2 1 -1 1 6 -1 -4 ;B = 13,-9, 6, 0;%disp(A) dim = size(A,1); L = zeros(dim, dim);y= x = zeros(dim, 1);for i = 1:1:dim for j= 1: 1: dim if i = (j + 1) L(i,j) = (A(i,j) - sum(L(i,1:(j-1).*L(j,1:(j-1)/L(j,j); else L(j,j) = sqrt(A(j,j) - sum(L(j

22、,1:j-1).*L(j,1:j-1); end endend L_T = L.; for i = 1: 1:dim if i = 1 y(1) = B(1) / L(1,1); else y(i) = (B(i) -sum(L(i,1:(i-1).*y(1:(i-1)/L(i,i); endend for i = dim:-1: 1 if i = dim x(i) = y(i)/L(i,i); else x(i) = (y(i) - sum(L(i+1):dim,i).*x(i+1):dim)/L(i,i); end % disp() %disp(x(i)enddisp(x)disp(x)计

23、算结果x 52.2500 -38.7500 30.7500 -52.75006. 用追赶法编制程序求解方程组Ax=b，其中， .MATLAB程序A=4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6;a=0 3 2 -1;c=2 1 3;b=4 -2 5 6;d=6 2 10 5;n=length(b);u0=0;y0=0;a(1)=0;L(1)=b(1)-a(1)*u0;y(1)=(d(1)-y0*a(1)/L(1);u(1)=c(1)/L(1);for i=2:(n-1)L(i)=b(i)-a(i)*u(i-1);y(i)=(d(i)-y(i-1)*a(i)/L(i);u(

24、i)=c(i)/L(i);endL(n)=b(n)-a(n)*u(n-1);y(n)=(d(n)-y(n-1)*a(n)/L(n);x(n)=y(n);for i=(n-1):-1:1x(i)=y(i)-u(i)*x(i+1);endx计算结果ans = 1 1 1 17. 已知编程求解矩阵A的QR分解.MATLAB程序function Q,R=qrhouseholder(A)dim=size(A,1);R=A;Q=eye(dim);for i=1:(dim-1) x=R(i:dim,i); y=1;zeros(dim-i,1); Ht=householder(x,y); H=blkdiag(

25、eye(i-1),Ht); Q=Q*H; R=H*R;end end function H,rho=householder(x,y)x=x(:);y=y(:);if length(x)=length(y) error(The Column Vectors X and Y Must Have The Same Length!);endrho=-sign(x(1)*norm(x)/norm(y);y=rho*y;v=(x-y)/norm(x-y);I=eye(length(x);H=I-2*v*v; End 计算结果C = 1.0000 1.0000 0 0 -1.0000 3.0000 -0.5

26、000 0.5000 -2.0000 2.0000 1.5000 0.5000 -2.0000 2.0000 -0.5000 2.5000 Q,R= qrhouseholder(C)Q = -0.3162 -0.7071 -0.2582 0.5774 0.3162 -0.7071 0.2582 -0.5774 0.6325 -0.0000 -0.7746 -0.0000 0.6325 -0.0000 0.5164 0.5774R = -3.1623 3.1623 0.4743 2.0555 0.0000 -2.8284 0.3536 -0.3536 0.0000 0.0000 -1.5492

27、1.0328 -0.0000 -0.0000 -0.0000 1.15478. 分别应用Jacobi 迭代法和Gauss-Seidel 迭代法求解如下方程组Jacobi法主程序functionx,k,succ=Jacobi(A,b,eps,it_max)n=length(A);k=0; x=zeros(n,1);y=zeros(n,1);succ=1; while 1 for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-A(i,j)*x(j); end end if k=it_max succ=0; return; end y(i)=y(i)/A(i,i); end if norm(y-x,2) A=4,-1,1;4,8,1;-2,1,5; b=7,-21,15; x,k,succ=Jacobi(A,b,1e-6,1000)x = 0.0606 -3.1111 3.6465k = 21succ = 1G-