搜索
数值分析第五章第一题：LU分解法：建立m文件functionh1=zhijieLU(A,b)%h1各阶主子式的行列式值[nn]=size(A);RA=rank(A);ifRA~=ndisp('请注意：因为A的n阶行列式h1等于零，所以A不能进行LU分解。A的秩RA如下：')RA,h1=det(A);returnendifRA==nforp=1:nh(p)=det(A(1:p,1:p));endh1=h(1:n);fori=1:nifh(1,i)==0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：')h1;RAreturnendendifh(1,i)~=0disp('请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：')forj=1:nU(1,j)=A(1,j);endfork=2:nfori=2:nforj=2:nL(1,1)=1;L(i,i)=1;ifijL(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);endendendendh1;RA,U,L,X=inv(U)*inv(L)*bendend输入：A=[10-701;-32.09999962;5-15-1;2102];b=[8;5.900001;5;1];h1=zhijieLU(A,b)输出：请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：RA=4U=10.0000-7.000001.000002.10006.00002.300000-2.1429-4.23810-0.0000012.7333L=1.0000000-0.30001.0000000.50001.19051.0000-0.00000.20001.14293.20001.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=[Ab];n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;ifzhicha0disp('请注意：因为RA~=RB，所以方程组无解')returnwarningoffMATLAB:return_outside_of_loopendifRA==RBifRA==ndisp('请注意：因为RA=RB，所以方程组有唯一解')X=zeros(n,1);C=zeros(1,n+1);forp=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;fork=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);forq=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RBn，所以方程组有无穷多解')endend输入：A=[10-701;-32.09999962;5-15-1;2102];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）输入：formatcompactA=[3.016.031.99;1.274.16-1.23;0.987-4.819.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.006.031.99;1.274.16-1.23;0.990-4.819.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第三题：输入：clearA=[10787;7565;86109;75910];b=[32233331]’;dA=det(A),lamda=eig(A),Ac2=cond(A,2)输出：dA=1.0000lamda=0.01020.84313.858130.2887Ac2=2.9841e+03下面分析误差性态：建立m文件：functionAcp=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=A\b;deltaA=jA-A;pndA=norm(deltaA,p);deltab=jb-b;pndb=norm(deltab,p);ifpndb0jX=A\jb;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;endifpndA0jX=jA\b;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)&(dA0.1)disp('请注意：AX=b是病态的，A的p条件数Acp，A的行列式值dA，解X，近似解jX，解的相对误差xX，解的相对误差估计Xgxx，b或A的相对误差xAb依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'elsedisp('请注意：AX=b是良态的，A的p条件数Acp，A的行列式值dA，解X，近似解jX，解的相对误差xX，解的相对误差估计Xgxx，b或A的相对误差xAb依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'end输入：jA=[1078.17.2;7.085.0465;85.989.899;6.99599.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.00001.00001.00001.0000ans=-9.586318.3741-3.22583.5240xX=10.4661jxX=0.9842Xgxx=22.7396xAb=0.0076xAbj=0.0076Acp=2.9841e+03第四题：（1）输入：建立m文件:forn=2:6a=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=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,deltax=X-x输出：X=1.00001.00001.0000r=1.0e-15*0.222000deltax=1.0e-13*-0.02000.1221-0.1255n=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.3886n=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.00001.00001.00001.00001.0000r=1.0e-15*00.2220000.1110deltax=1.0e-11*-0.00350.0524-0.19370.2591-0.1148n=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.1110deltax=1.0e-11*-0.00350.0524-0.19370.2591-0.1148n=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.1677n=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.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.0978n=9;H=hilb(n);x=(linspace(1,1,n))';b=H*x;[RA,RB,n,X]=gauss(H,b)r=b-H*X,del