数值分析(第五版)计算实习题第五章作业.doc

1、数值分析第五章第一题：LU 分解法：建立 m 文件function h1=zhijieLU(A,b) %h1各阶主子式的行列式值n n=size(A);RA=rank(A);if RA=ndisp(请注意：因为A的n 阶行列式h1等于零，所以A不能进行LU分解。A 的秩RA如下：)RA,h1=det(A);returnendif RA=nfor p=1:nh(p)=det(A(1:p,1:p);endh1=h(1:n);for i=1:nif h(1,i)=0disp(请注意：因为A 的r 阶主子式等于零，所以A 不能进行 LU分解。A的秩RA和各阶顺序主子式h1依次如下：)h1;RAretu

2、rnendendif h(1,i)=0disp(请注意：因为A的r阶主子式都不等于零，所以A 能进行 LU分解。A的秩RA和各阶顺序主子式h1依次如下：)for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:nfor j=2:nL(1,1)=1;L(i,i)=1;if ijL(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k)/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendend

3、h1;RA,U,L, X=inv(U)*inv(L)*bendend输入： A=10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2; b=8;5.900001;5;1; h1=zhijieLU(A,b)输出：请注意：因为 A 的 r 阶主子式都不等于零，所以 A 能进行 LU 分解。A 的秩 RA 和各阶顺序主子式 h1 依次如下：RA =4U =10.0000 -7.0000 0 1.00000 2.1000 6.0000 2.30000 0 -2.1429 -4.23810 -0.0000 0 12.7333L =1.0000 0 0 0-0.3000

4、1.0000 0 00.5000 1.1905 1.0000 -0.00000.2000 1.1429 3.2000 1.0000X =-0.2749-1.32981.29691.4398h1 =10.0000 -0.0000 -150.0001 -762.0001列主元高斯消去法：建立 m 文件function RA,RB,n,X=liezhu(A,b)B=A b;n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;if zhicha0disp(请注意：因为RA=RB ，所以方程组无解 )returnwarning off MATLAB:return

5、_outside_of_loopendif RA=RBif RA=ndisp(请注意：因为RA=RB，所以方程组有唯一解)X=zeros(n,1);C=zeros(1,n+1);for p=1:n-1Y,j=max(abs(B(p:n,p);C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum

6、(A(q,q+1:n)*X(q+1:n)/A(q,q);endelsedisp(请注意：因为RA=RB A=10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2; b=8;5.900001;5;1; RA,RB,n,X=liezhu(A,b),H=det(A)输出：请注意：因为 RA=RB，所以方程组有唯一解RA =4RB =4n =4X =0.0000-1.00001.00001.0000H =-762.0001第二题：建立列主元高斯消去法 m 文件（题一中已有）（1 ）输入： format compact A=3.01 6.03 1.99;1.27 4.1

7、6 -1.23;0.987 -4.81 9.34; b=1;1;1; RA,RB,n,X=liezhu(A,b),h=det(A),C=cond(A)输出：请注意：因为 RA=RB，所以方程组有唯一解RA =3RB =3n =3X =1.0e+03 *1.5926-0.6319-0.4936h =-0.0305C =3.0697e+04(2)输入： A=3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34; b=1;1;1; RA,RB,n,X=liezhu(A,b),h=det(A)输出：请注意：因为 RA=RB，所以方程组有唯一解RA =3RB =

8、3n =3X =119.5273-47.1426-36.8403h =-0.4070 第三题：输入： clear A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10; b=32 23 33 31; dA=det(A),lamda=eig(A),Ac2=cond(A,2)输出：dA =1.0000lamda =0.01020.84313.858130.2887Ac2 =2.9841e+03下面分析误差性态：建立 m 文件：function Acp=pjwc(A,jA,b,jb,p)%Acp矩阵A的p条件数 cond%pjwc： p范数解的误差性态分析%jA是A的近似矩阵jA

9、=A+A，jb=b+bAcp=cond(A,p);dA=det(A);X=Ab;deltaA=jA-A;pndA=norm(deltaA,p);deltab=jb-b;pndb=norm(deltab,p);if pndb0jX=Ajb;Pnb=norm(b,p);pnjx=norm(jX,p);deltaX=jX-X;pnjdX=norm(deltaX,p);jxX=pnjdX/pnjX;pnX=norm(X,p);xX=pnjdX/pnX;pndb=norm(deltab,p);xAb=pndb/pnb;pnbj=norm(jb,p);xAbj=pndb/pnbj;Xgxx=Acp*xAb

10、;endif pndA0jX=jAb;deltaX=jX-X;pnX=norm(X,p);pnjdX=norm(deltaX,p);pnjX=norm(jX,p);jxX=pnjdX/pnjX;xX=pnjdX/pnX;pnjA=norm(jA,p);pnA=norm(A,p);pndA=norm(deltaA,p);xAbj=pndA/pnjA;xAb=pndA/pnA;Xgxx=Acp*xAb;endif (Acp50)7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98; jb=b;p=2; Acp=pjwc(A,jA,b,jb,p)输出：请注意：AX=b是

11、良态的，A的p条件数 Acp，A的行列式值dA，解 X，近似解jX，解的相对误差xX，解的相对误差估计Xgxx，b或A的相对误差xAb依次如下：Acp =2.9841e+03dA =1.0000ans =1.0000 1.0000 1.0000 1.0000ans =-9.5863 18.3741 -3.2258 3.5240xX =10.4661jxX =0.9842Xgxx =22.7396xAb =0.0076xAbj =0.0076Acp =2.9841e+03第四题：（1）输入：建立 m 文件:for n=2:6a=hilb(n);pnH(n-1)=cond(a,inf);endpn

12、Hn=2:6;plot(n,pnH);可见条件数随着 n 的增大而急剧增大（2 ）输入： n=2;H=hilb(n); x=(linspace(1,1,n); b=H*x; RA,RB,n,X=gauss(H,b)输出：请注意：因为 RA=RB，所以方程组有唯一解RA =2RB =2n =2X =1.00001.0000输入： r=b-H*X,deltax=X-x输出：r =00deltax =1.0e-15 *0.4441-0.6661输入： n=3;H=hilb(n); x=(linspace(1,1,n); b=H*x; RA,RB,n,X=gauss(H,b) r=b-H*X,delt

13、ax=X-x输出：X =1.00001.00001.0000r =1.0e-15 *0.222000deltax =1.0e-13 *-0.02000.1221-0.1255 n=4;H=hilb(n); x=(linspace(1,1,n); b=H*x; RA,RB,n,X=gauss(H,b) r=b-H*X,deltax=X-xX =1.00001.00001.00001.0000r =1.0e-15 *-0.44410-0.11100deltax =1.0e-12 *-0.02220.2485-0.59800.3886 n=5;H=hilb(n); x=(linspace(1,1,n

14、); b=H*x; RA,RB,n,X=gauss(H,b) r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.0000r =1.0e-15 *00.2220000.1110deltax =1.0e-11 *-0.00350.0524-0.19370.2591-0.1148 n=6;H=hilb(n); x=(linspace(1,1,n); b=H*x; RA,RB,n,X=gauss(H,b) r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.0000r =1.0e-15 *00.2220000.111

15、0deltax =1.0e-11 *-0.00350.0524-0.19370.2591-0.1148 n=7;H=hilb(n); x=(linspace(1,1,n); b=H*x; RA,RB,n,X=gauss(H,b) r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.00001.0000r =1.0e-15 *000.2220000.1110deltax =1.0e-09 *-0.00080.0219-0.14820.3854-0.42540.1677 n=8;H=hilb(n); x=(linspace(1,1,n); b=H*x; R

16、A,RB,n,X=gauss(H,b) r=b-H*X,deltax=X-xX =1.00001.00001.00001.00001.00001.00001.00001.0000r =1.0e-15 *00-0.222000-0.1110-0.11100deltax =1.0e-06 *-0.00000.0018-0.02360.1279-0.34420.4870-0.34660.0978 n=9;H=hilb(n); x=(linspace(1,1,n); b=H*x; RA,RB,n,X=gauss(H,b) r=b-H*X,deltax=X-xX =1.00001.00001.00001

17、.00001.00001.00001.00001.00001.0000r =1.0e-15 *0.44410.2220-0.222000.22200.22200-0.11100deltax =1.0e-04 *-0.00000.0002-0.00280.0197-0.07220.1471-0.16870.1017-0.0251 n=10;H=hilb(n); x=(linspace(1,1,n); b=H*x; RA,RB,n,X=gauss(H,b) r=b-H*X,deltax=X-xX =1.00001.00001.00001.00000.99991.00030.99961.00040.99981.0000r =1.0e-15 *0.444100-0.22200.222000-0.11100.11100.1110deltax =1.0e-03 *-0.00000.0001-0.00230.0205-0.09740.2669-0.43690.4214-0.22090.0485