数值分析第二章上机

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 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,:

2、)=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

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); en

4、dhl=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,

5、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.

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

7、 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

8、); 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=

9、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',xabj'else disp('请注意： ax=b是良态的，a的p条件数acp,a的行列式值da，解x，近似解jx，解的相对误差

10、jxx，解的相对误差估计值xgxx，b或a的相对误差xab依次如下：') acp,da,x',jx',xx',jxx',xgxx',xab',xabj'end习题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，解的

11、相对误差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=

12、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

13、,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

14、工作窗口输入程序>>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=

15、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,

16、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

17、9;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

18、是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 请注意：雅可比迭代次数已经超过最大迭代次数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工作窗口输入程序>>

19、; 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('请注意：因为对

20、角矩阵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的分解矩阵

21、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)运行后输出结果请注

22、意：因为对角矩阵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=

23、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.07