数值分析第二章上机.doc

资源ID：11182686 资源大小：69.50KB 全文页数：9页
资源格式： DOC 下载积分：4元

快捷下载

会员登录下载

微信登录下载

三方登录下载：

微信扫一扫登录

下载资源需要4元

邮箱/手机：
温馨提示：	用户名和密码都是您填写的邮箱或者手机号，方便查询和重复下载（系统自动生成）
支付方式：
验证码：	换一换

加入VIP免费专享

账号：
密码：
验证码：	换一换
当日自动登录忘记密码？

友情提示

1、下载资料失败解决办法

2、PDF文件下载后，可能会被浏览器默认打开，此种情况可以点击浏览器菜单，保存网页到桌面，就可以正常下载了。

3、本站不支持迅雷下载，请使用电脑自带的IE浏览器，或者360浏览器、谷歌浏览器下载即可。

4、本站资源下载后的文档和图纸-无水印,预览文档经过压缩，下载后原文更清晰。

5、试题试卷类文档，如果标题没有明确说明有答案则都视为没有答案，请知晓。

网站客服

侵权投诉

数值分析第二章上机.doc

用列主元消元法解线性方程组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 zhica>0,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,:)=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<n，所以此方程组有无穷多解.) endend习题3.3 3.在MATLAB工作窗口输入程序>> 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 = 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 returnendendif 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 i>j 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 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;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=norm(dertaA,p);dertab=b-jb; Pndb=norm(dertab,p); if Pndb>0 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 PndA>0 jX=jAb;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 (Acp>50)&&(dA<0.1) disp(请注意：AX=b是病态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：) Acp,dA,X,jX,xX,jxX,Xgxx,xAb,xAbjelse disp(请注意： AX=b是良态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：) Acp,dA,X,jX,xX,jxX,Xgxx,xAb,xAbjend习题3.61.在MATLAB工作窗口输入程序>>A=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,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主程序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); end X;djwcX=norm(X-X0,P); xdwcX=djwcX/(norm(X,P)+eps); X0=X;X1=Ab; if (djwcX<wucha)&&(xdwcX<wucha) disp(请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：) break endendif (djwcX>wucha)&&(xdwcX>wucha) 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是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛请注意：雅可比迭代次数已经超过最大迭代次数max1 a = -21 -8 -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,inf,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取最大迭代次数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=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); 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 (djwcX<wucha)|(xdwcX<wucha) break else k;X1;k=k+1;X=X1; end end if (djwcX<wucha)|(xdwcX<wucha) disp(请注意：高斯-塞德尔迭代收敛,此A的分解矩阵D,U,L和方程组的精确解jX和近似解X如下： ) else disp(请注意：高斯-塞德尔迭代的结果没有达到给定的精度，并且迭代次数已经超过最大迭代次数max1,方程组的精确解jX和迭代向量X如下： ) X=X;jX=jX endendX=X;D,U,L,jX=jX习题4.33.（1）在MATLAB工作窗口输入程序>> 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 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 = 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);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

注意事项

本文（数值分析第二章上机.doc）为本站会员（scccc）主动上传，三一文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若此文所含内容侵犯了您的版权或隐私，请立即通知三一文库（点击联系客服），我们立即给予删除！

温馨提示：如果因为网速或其他原因下载失败请重新下载，重复下载不扣分。