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

数值分析第五章第一题：LU分解法：建立m文件function h1=zhijieLU(A,b) %h1各阶主子式的行列式值 n n=size(A);RA=rank(A);if RA=n disp(请注意：因为A的n阶行列式h1等于零，所以A不能进行LU分解。A的秩RA如下：) RA,h1=det(A); returnendif RA=n for p=1:n h(p)=det(A(1:p,1:p); end h1=h(1:n); for i=1:n if h(1,i)=0 disp(请注意：因为A的r阶主子式等于零，所以A不能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：) h1;RA return end end if h(1,i)=0 disp(请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：) for j=1:n U(1,j)=A(1,j); end for k=2:n for i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if ij L(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); else U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end h1;RA,U,L, X=inv(U)*inv(L)*b endend输入： 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.0000 0 2.1000 6.0000 2.3000 0 0 -2.1429 -4.2381 0 -0.0000 0 12.7333L = 1.0000 0 0 0 -0.3000 1.0000 0 0 0.5000 1.1905 1.0000 -0.0000 0.2000 1.1429 3.2000 1.0000X = -0.2749 -1.3298 1.2969 1.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 zhicha0 disp(请注意：因为RA=RB，所以方程组无解) return warning off MATLAB:return_outside_of_loopendif RA=RB if RA=n disp(请注意：因为RA=RB，所以方程组有唯一解) X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1 Y,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:n m=B(k,p)/B(p,p); B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); end else disp(请注意：因为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.0000 1.0000 1.0000H =-762.0001第二题：建立列主元高斯消去法m文件（题一中已有）（1）输入： format compact A=3.01 6.03 1.99;1.27 4.16 -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 = 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.0102 0.8431 3.8581 30.2887Ac2 = 2.9841e+03下面分析误差性态：建立m文件：function Acp=pjwc(A,jA,b,jb,p)%Acp矩阵A的p条件数cond%pjwc：p范数解的误差性态分析%jA是A的近似矩阵jA=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 pndb0 jX=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;endif pndA0 jX=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)&(dA jA=10 7 8.1 7.2;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是良态的，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:6 a=hilb(n); pnH(n-1)=cond(a,inf);endpnHn=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 = 0 0deltax = 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,deltax=X-x输出：X = 1.0000 1.0000 1.0000r = 1.0e-15 * 0.2220 0 0deltax = 1.0e-13 * -0.0200 0.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.0000 1.0000 1.0000 1.0000r = 1.0e-15 * -0.4441 0 -0.1110 0deltax = 1.0e-12 * -0.0222 0.2485 -0.59800.3886 n=5;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.0000 1.0000 1.0000 1.0000 1.0000r = 1.0e-15 * 0 0.2220 0 0 0.1110deltax = 1.0e-11 * -0.0035 0.0524 -0.1937 0.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.0000 1.0000 1.0000 1.0000 1.0000r = 1.0e-15 * 0 0.2220 0 0 0.1110deltax = 1.0e-11 * -0.0035 0.0524 -0.1937 0.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.0000 1.0000 1.0000 1.0000 1.0000 1.0000r = 1.0e-15 * 0 0 0.2220 0 0 0.1110deltax = 1.0e-09 * -0.0008 0.0219 -0.1482 0.3854 -0.4254 0.1677 n=8;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.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000r = 1.0e-15 * 0 0 -0.2220 0 0 -0.1110 -0.1110 0deltax = 1.0e-06 * -0.0000 0.0018 -0.0236 0.1279 -0.3442 0.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.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000r = 1.0e-15 * 0.4441 0.2220 -0.2220 0 0.2220 0.2220 0 -0.1110 0deltax = 1.0e-04 * -0.0000 0.0002 -0.0028 0.0197 -0.0722 0.1471 -0.1687 0.1017 -0.0251 n=10;H=hilb(n); x=(linspace(1,1,n); b=