《数值分析第二章上机.doc》由会员分享，可在线阅读，更多相关《数值分析第二章上机.doc（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是严格对角
12、占优的，此方程组有唯一解，且雅可比迭代收敛 请注意：雅可比迭代收敛,此方程组的精确解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是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 请注意：雅可比迭代次数已经超过最大迭代次数max1 a = -21 -8 -
13、1jX = 0.2159 1.0711 1.0974X = 0.2152 1.0697 1.0959（2）取最大迭代次数Max1=100在MATLAB工作窗口输入程序A=15 -1 -2;-1 10 -2;-1 -1 0.5;b=7.2;8.3;4.2;X0=0 0 0;X=jacdd(A,b,X0,inf,0.001,100)运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛取最大迭代次数Max1=5在MATLAB工作窗口输入程序 A=15 -1 -2;-1 10 -2;-1 -1 0.5;b=7.2;8.3;4.2;X0=0 0 0;X=jacdd(A,b,X0,i
14、nf,0.001,5)运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛（3）取最大迭代次数Max1=100在MATLAB工作窗口输入程序 A=10 -1 -2;-1 10 -2;-1 -1 5;b=7.2;8.3;4.2;X0=0 0 0;X=jacdd(A,b,X0,inf,0.001,100)运行后输出结果请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：a = -8 -8 -1jX = 1.1000 1.2000 1.3000X = 1.0998 1.1998 1.2998取最大迭
15、代次数Max1=5在MATLAB工作窗口输入程序 A=10 -1 -2;-1 10 -2;-1 -1 5;b=7.2;8.3;4.2;X0=0 0 0;X=jacdd(A,b,X0,inf,0.001,5)运行后输出结果请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 请注意：雅可比迭代次数已经超过最大迭代次数max1 a = -8 -8 -1jX = 1.1000 1.2000 1.3000X = 1.0951 1.1951 1.2941（4）取最大迭代次数Max1=100在MATLAB工作窗口输入程序A=1 2 3;2 5 2;3 1 5;b=14;18;20;X0=
16、0 0 0;X=jacdd(A,b,X0,inf,0.001,100)运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛取最大迭代次数Max1=5在MATLAB工作窗口输入程序 A=1 2 3;2 5 2;3 1 5;b=14;18;20;X0=0 0 0;X=jacdd(A,b,X0,inf,0.001,5)运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛用高斯-塞德尔迭代定义解线性方程组 的MATLAB主程序function X=gsdddy(A,b,X0,P,wucha,max1) D=diag(diag(A);U=-triu(A,1);
17、 L=-tril(A,-1); dD=det(D);if dD=0 disp(请注意：因为对角矩阵D奇异，所以此方程组无解.)else disp(请注意：因为对角矩阵D非奇异，所以此方程组有解.) iD=inv(D-L); B2=iD*U;f2=iD*b;jX=Ab; X=X0; n m=size(A); for k=1:max1 X1= B2*X+f2; djwcX=norm(X1-X,P); xdwcX=djwcX/(norm(X,P)+eps); if (djwcXwucha)|(xdwcXwucha) break else k;X1;k=k+1;X=X1; end end if (dj
18、wcXwucha)|(xdwcX A=11 -1 -2;-1 10 -2;-1 -1 0.5;b=7.2 ;8.3;4.2;X0=0 0 0;X=gsdddy(A,b,X0,inf,0.00001,100)运行后输出结果请注意：因为对角矩阵D非奇异，所以此方程组有解.请注意：高斯-塞德尔迭代收敛,此A的分解矩阵D,U,L和方程组的精确解jX和近似解X如下： D = 11.0000 0 0 0 10.0000 0 0 0 0.5000U = 0 1 2 0 0 2 0 0 0L = 0 0 0 1 0 0 1 1 0jX = 15.8529 17.3941 74.8941X = 15.8518
19、17.3928 74.8892（2）在MATLAB工作窗口输入程序 A=4 4 -5 7;2 -8 3 -2;4 5 -13 16;7 -2 2 3;b=5;2;-1;21;X0=0 0 0 0;X=gsdddy(A,b,X0,inf,0.00001,100)运行后输出结果请注意：因为对角矩阵D非奇异，所以此方程组有解.请注意：高斯-塞德尔迭代的结果没有达到给定的精度，并且迭代次数已经超过最大迭代次数max1,方程组的精确解jX和迭代向量X如下： jX = 0.3446 0.7555 4.7894 3.5066D = 4 0 0 0 0 -8 0 0 0 0 -13 0 0 0 0 3U =
20、0 -4 5 -7 0 0 -3 2 0 0 0 -16 0 0 0 0L = 0 0 0 0 -2 0 0 0 -4 -5 0 0 -7 2 -2 0jX = 0.3446 0.7555 4.7894 3.5066X = 1.0e+26 * -0.7417 -0.3374 0.0768 1.4545用谱半径判别超松弛迭代法产生的迭代序列的敛散性的MATLAB主程序function H=ddpbj(A,om) D=diag(diag(A);U=-triu(A,1); L=-tril(A,-1); iD=inv(D-om*L); B2=iD*(om*U+(1-om)*D); H=eig(B2);
21、mH=norm(H,inf);if mH=1 disp(请注意：因为谱半径不小于1，所以超松弛迭代序列发散,谱半径mH和B的所有的特征值H如下：)else disp(请注意：因为谱半径小于1，所以超松弛迭代序列收敛,谱半径mH和B的所有的特征值H如下：)endmH习题4.41.（1）当取=1.15时，在MATLAB工作窗口输入程序 A=7 1 -1 -2;2 8 1 3;1 -2 -4 -1;-1 3 2 7;H=ddpbj(A,1.15)运行后输出结果请注意：因为谱半径小于1，所以超松弛迭代序列收敛,谱半径mH和B的所有的特征值H如下：mH = 0.1608H = 0.0715 + 0.1440i 0.0715 - 0.1440i -0.1308 + 0.0498i -0.1308 - 0.0498i（2）当取=1.15时，在MATLAB工作窗口输入程序 A=7 1 -1 -2;2 8 1 3;1 -2 -4 -1;-1 3 2 7;H=ddpbj(A,5)运行后输出结果请注意：因为谱半径不小于1，所以超松弛迭代序列发散,谱半径mH和B的所有的特征值H如下：mH = 13.1892H = -13.1892 + 0.0000i -2.6969 + 0.0000i 0.2460 + 2.6714i 0.2460 - 2.6714i