《数值分析第二章上机》由会员分享，可在线阅读，更多相关《数值分析第二章上机（9页珍藏版）》请在金锄头文库上搜索。
1、用列主元消元法解线性方程组AX=b的MATLAB程序functionRA,RB,n,X=liezhu(A,b)B=A b;n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica0,disp(请注意：因为RA=RB，所以此方程组无解.)returnendif RA=RB if RA=ndisp(请注意：因为RA=RB=n，所以此方程组有唯一解.)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,:)=
2、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);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(A(q,q+1:n)*X(q+1:n)/A(q,q);end else disp(请注意：因为RA=RB A=1 1 1;1 3 9;1 7 49;b=6;5;2;RA,RB,n,X=liezhu(A,b)运行后输出结果请注意：因为RA=RB=n，所以此方程组有唯一解.RA = 3RB = 3n = 3X =
3、6.3750-0.3333-0.0417将矩阵A进行直接LU分解的MATLAB程序function hl=zhjLU(A)n n=size(A);RA=rank(A);if RA=ndisp(请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:),RA,hl=det(A);returnendif RA=n for p=1:nh(p)=det(A(1:p,1:p); endhl=h(1:n);for i=1:nif h(1,i)=0 disp(请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：),hl;RA return
4、endendif h(1,i)=0disp(请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：)for j=1:nU(1,j)=A(1,j);endfor 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 en
5、d endend hl;RA,U,Lendend习题3.41.(1) 在MATLAB工作窗口输入程序 A=2 4 -6;1 5 3;1 3 2;hl=zhjLU(A)运行后输出结果请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA = 3U =2.0000 4.0000 -6.0000 0 5.0000 6.0000 0 0 3.8000L =1.0000 0 0 0.5000 1.0000 0 0.5000 0.2000 1.0000hl =2 6 18(2) 在MATLAB工作窗口输入程序 A=1 1 6;-1 2 9;1 -2 3;
6、hl=zhjLU(A) 运行后输出结果请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA =3U = 1.0000 1.0000 6.0000 0 2.0000 15.0000 0 0 19.5000L = 1.0000 0 0 -1.0000 1.0000 0 1.0000 -1.5000 1.0000hl = 1 3 36用P范数讨论AX=b解和A的性态的MATLAB程序function Acp=zpjwc(A,jA,b,jb,p) Acp=cond(A,p);dA=det(A);X=Ab; dertaA=A-jA; PndA=nor
7、m(dertaA,p);dertab=b-jb; Pndb=norm(dertab,p); if Pndb0 jX=Ajb;Pnb=norm(b,p);PnjX=norm(jX,p);dertaX=X-jX; PnjdX=norm(dertaX,p);jxX=PnjdX/PnjX;PnjX=norm(jX,p); PnX=norm(X,p);jxX=PnjdX/PnjX;xX=PnjdX/PnX; Pndb=norm(dertab,p); xAb=Pndb/Pnb;Pnbj=norm(jb,p);xAbj=Pndb/Pnbj; Xgxx=Acp*xAb; end if PndA0 jX=jAb
8、;dertaX=X-jX;PnX=norm(X,p); PnjdX=norm(dertaX,p); PnjX=norm(jX,p);jxX=PnjdX/PnjX;xX=PnjdX/PnX; PnjA=norm(jA,p);PnA=norm(A,p); PndA=norm(dertaA,p);xAbj=PndA/PnjA;xAb=PndA/PnA; Xgxx=Acp*xAb; endif (Acp50)&(dAA=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;jA=A;b=32 23 33 31;jb=32.1 22.9 22.2 30.9;Acp=zpjwc(A,jA,
9、b,jb,1)运行后输出结果请注意：AX=b是良态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：Acp = 4.4880e+03dA = 1.0000ans = 1.0000 1.0000 1.0000 1.0000ans = -99.8000 172.7000 -50.0000 31.6000xX = 88.5250jxX = 1.0000Xgxx = 418.6286xAb = 0.0933xAbj = 0.1027Acp = 4.4880e+03用雅可比迭代解线性方程组AX=b的MATLAB主程
10、序function X=jacdd(A,b,X0,P,wucha,max1)n m=size(A); for j=1:m a(j)=sum(abs(A(:,j)-2*(abs(A(j,j); end for i=1:n if a(i)=0 disp(请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛) return end endif a(i)0 disp(请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛)endfor k=1:max1 k; for j=1:m X(j)=(b(j)-A(j,1:j-1,j+1:m)*X0(1: j-1,j+1:m)/A(j,j
11、); end X;djwcX=norm(X-X0,P); xdwcX=djwcX/(norm(X,P)+eps); X0=X;X1=Ab; if (djwcXwucha)&(xdwcXwucha)&(xdwcXwucha) disp(请注意：雅可比迭代次数已经超过最大迭代次数max1)enda,X=X;jX=X1, 习题4.22.（1）取最大迭代次数Max1=100在MATLAB工作窗口输入程序A=23 -1 -2;-1 10 -2;-1 -1 5;b=1.7;8.3;4.2;X0=0 0 0;X=jacdd(A,b,X0,inf,0.001,100)运行后输出结果请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：a = -21 -8 -1jX = 0.2159 1.0711 1.0974X =0.2159 1.0710 1.0973取最大迭代次数Max=5在MATLAB工作窗口输入程序 A=23 -1 -2;-1 10 -2;-1 -1 5;b=1.7;8.3;4.2;X0=0 0 0;X=jacdd(A,b,X0,inf,0.001,5)运行后输出结果请注意：系数矩阵A是严格对角
