MATLAB解线性方程组的直接方法_第1页
MATLAB解线性方程组的直接方法_第2页
MATLAB解线性方程组的直接方法_第3页
MATLAB解线性方程组的直接方法_第4页
MATLAB解线性方程组的直接方法_第5页
已阅读5页,还剩8页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、在这章中我们要学习线性方程组的直接法,特别是适合用数学软件在计算机上求解 的方法.3.1 方程组的逆矩阵解法及其MATLA跑序3.1.3线性方程组有解的判定条件及其MATLA醒序判定线性方程组 A m迎X = b是否有解的MATLA典序function RA,RB,n=jiepb(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=RBif RA=ndisp('请注意:因为RA=RB=n ,所

2、以此方程组有唯一解 .') elsedisp('请注意:因为RA=RB<n ,所以此方程组有无穷多解.)endend例3.1.4判断下列线性方程组解的情况.如果有唯一解,则用表 3-2方法求解2x1 + 3x2 - X3 +5x4 =0,3x1 +4x2 - 5x3 +7x4 =0,(1)3X1 十 X2 + 2X3 7X4 =0, (2)J2X1 -3X2 + 3X3 - 2X4 = 0,4x1 + x2 - 3x3 +6x4 =0,4x1 +11x2 -13x3 +16x4 = 0,J1 2x2 +4x3 7X4 =0;7 X1 -2x2 + x3 +3x4 =0;&

3、#39;4X1 +2X2 -X3 =2,2x + y-z+w=1, 3X1 -X2 + 2x3 =10,(4)4x + 2y - 2z + w = 2,11 x1 + 3x2 =8;、2x + y z w=1.解 在MATLABT作窗口输入程序>> A=2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7;b= 0; 0; 0; 0; RA,RB,n=jiepb(A,b)运行后输出结果为请注意:因为 RA=RB=n ,所以此方程组有唯一解RA = 4,RB =4,n =4在MATLAB:作窗口输入>>X=Ab,运行后输出结果为 X =(0 0 0 0

4、).(2)在MATLABT作窗口输入程序>> A=3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3;b= 0; 0; 0; 0;RA,RB,n=jiepb(A,b)24.运行后输出结果请注意:因为RA=RB<n ,所以此方程组有无穷多解.RA =2,RB =2,n =4(3) 在MATLAB:作窗口输入程序>> A=4 2 -1;3 -1 2;113 0; b=2;10;8; RA,RB,n=jiepb(A,B)运行后输出结果请注意:因为RA=RB,所以此方程组无解.RA =2,RB =3,n =3(4)在MATLABT作窗口输入程序

5、>> A=2 1 -1 1;4 2 -2 1;2 1 -1 -1;b=1; 2; 1; RA,RB,n=jiepb(A,b)运行后输出结果请注意:因为RA=RB<n,所以此方程组有无穷多解.RA =2,RB =2,n =33.2三角形方程组的解法及其MATLA跑序3.2.2解三角形方程组的MATLAB?序解上三角形线性方程组AX = b的MATLA叁序function RA,RB,n,X=shangsan(A,b)B=A b; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:

6、因为RA=RB,所以此方程组无解.') return endif RA=RBif RA=ndisp('请注意:因为RA=RB=n ,所以此方程组有唯一解.')X=zeros(n,1); X(n)=b(n)/A(n,n);for k=n-1:-1:1X(k)=(b(k)-sum(A(k,k+1:n)*X(k+1:n)/A(k,k); end else disp('请注意:因为RA=RB<n ,所以此方程组有无穷多解.)end end例3.2.2 用解上三角形线性方程组的 MATLAB程序解方程组P5x1 -x2 +2 x3 +3x4 =20,-2x2 

7、9; 7x3 - 4x4 = -7, .6x3 5x4 = 4,3x4 = 6.解 在MATLABT作窗口输入程序>>A=5 -1 2 3;0 -2 7 -4;0 0 6 5;0 0 0 3;b=20; -7; 4;6;RA,RB,n,X=shangsan(A,b)运行后输出结果请注意:因为RA=RB=n,所以此方程组有唯一解.RA = RB =4,4,n =25.4,X =2.4 -4.0 -1.0 2.0'3.3 高斯(Gaus9消元法和列主元消元法及其MATLA能序3.3.1 高斯消元法及其MATLAB?序用高斯消元法解线性方程组AX =b的MATLA翼序functi

8、on RA,RB,n,X=gaus(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=RBif RA=ndisp('请注意:因为RA=RB=n ,所以此方程组有唯一解.')X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1for k=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)

9、;end endb=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); endelsedisp('请注意:因为RA=RB<n ,所以此方程组有无穷多解.)endend例3.3.2用高斯消元法和 MATLA叁序求解下面的非齐次线性方程组,并且用逆矩阵解方程组的方法验证.Xi X2 * X3 - 3x4 1,一 X2 - X3 + X4 = 0, 2x1 -2x2 -4x3 6x4 =-1, I x1 - 2x2 - 4x3x4 - -1.

10、解在MATLABT作窗口输入程序>> A=1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1; b=1;0; -1;-1; RA,RB,n,X =gaus (A,b)运行后输出结果RA =4 RB =4 n =4请注意:因为RA=RB=n ,所以此方程组有唯一解X = 0 -0.5000 0.5000 03.3.2 列主元消元法及其MATLAB?序用列主元消元法解线性方程组AX= b的MATLA翼序function RA,RB,n,X=liezhu(A,b)B=A b; n=length(b); RA=rank(A);26.RB=rank(B);zhi

11、ca=RB-RA;if zhica>0,disp('请注意:因为RA=RB,所以此方程组无解.,) returnendif RA=RBif 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:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);end endb=

12、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); endelsedisp('请注意:因为RA=RB<n ,所以此方程组有无穷多解.,)endend例3.3.3用列主元消元法解线性方程组的MATLAB程序解方程组-X 2 - X3 + X4 = 0,x1一 X2 + X33x4 =1,|2Xi2 x2 -4 X3+ 6x4 = -11.Xi-2x2 -4x3+ x4 = 1解在MATLABT作窗口输入程序>> A=0 -1-

13、1 1;1 -1 1 -3;2 -2 -4 6;1 -2 -4 1;b=0;1;-1;-1; RA,RB,n,X=liezhu(A,b)运行后输出结果请注意:因为RA=RB=n ,所以此方程组有唯一解.RA = 4 , RB = 4 , n = 4 , X =0 -0.5 0.5 0'3.4 LU分解法及其MATLA能序3.4.1 判断矩阵LU分解的充要条件及其MATLAB?#判断矩阵A能否进行LU分解的MATLA翼序function hl=pdLUfj(A)n n =size(A); RA=rank(A);if RA=ndisp('请注意:因为A的n阶行列式hl等于零,所以A

14、不能进行LU分解.A的秩 RA如下:'),RA,hl=det(A);returnendif RA=nfor p=1:n,h(p)=det(A(1:p, 1:p);,endhl=h(1:n);for i=1:nif h(1,i)=0disp('请注意:因为A白r阶主子式等于零,所以 A不能进行LU分 解.A的秩RA和各阶顺序主子式值 hl依次如下:,),hl;RA, return27.end end ifh(1,i)=0disp('请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值 hl依次如下:)hl;RA end end3 ,12

15、3 '1123、7;127;(3)1236d56d56例3.4.1判断下列矩阵能否进行 LU分解,并求矩阵的秩1211245解 (1)在MATLAB作窗口输入程序>> A=1 2 3;1 12 7;4 5 6;hl=pdLUfj(A)运行后输出结果为A能进行LU分解.A的秩RA和请注意:因为A的各阶主子式都不等于零,所以 各阶顺序主子式值hl依次如下:RA = 3 , hl = 1 10 -48(2)在MATLAS:作窗口输入程序>> A=1 2 3;1 2 7;4 5 6;hl=pdLUfj(A)运行后输出结果为请注意:因为A的r阶主子式等于零,所以 A不能进

16、行LU分解.A的秩RA和各阶 顺序主子式值hl依次如下:RA = 3 , hl =1 0 12(3)在MATLAS:作窗口输入程序>> A=1 2 3;1 2 3;4 5 6;hl=pdLUfj(A)运行后输出结果为请注意:因为A的n阶行列式hl等于零,所以A不能进行LU分解.A的秩RA如下RA = 2 , hl = 03.4.2 直接LU分解法及其MATLA醒序将矩阵A进行直接LU分解的MATLA叁序function hl=zhjLU(A)n n =size(A); RA=rank(A);if RA=ndisp('请注意:因为A的n阶行列式hl等于零,所以A不能进行LU分

17、解.A 的秩 RA如下:'),RA,hl=det(A);returnendif RA=nfor p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp('请注意:因为A的r阶主子式等于零,所以A不能进行LU分解.A 的秩RA和各阶顺序主子式值hl依次如下:'),hl;RAreturnend end ifh(1,i)=0disp('请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A28.的秩RA和各阶顺序主子式值hl依次如下:)for j=1:nU(1,j尸A(1,j);endfo

18、r k=2:nfor i=2:n for j=2:nL(1,1)=1;L(i,i)=1; if i>jL(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); end end end end hl;RA,U,L end end例3.4.3用矩阵进行直接 LU分解的MATLAB1011001212040程序分解矩阵0 11 .33;解 在MATLAB:作窗口输入程序>&

19、gt; A=1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3; hl=zhjLU(A)运行后输出结果请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:RA = 4L = 1 0 0 0U = 1 0 2 001000 1 0 112100 0 2 10101hl = 1 1 2 43.4.4 判断正定对称矩阵的方法及其 MATLAB?序 判断矩阵A是否是正定对称矩阵的 MATLA叁序function hl=zddc(A)n n =size(A);for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n)

20、;zA=A'for i=1:nif h(1,i)<=0disp(,请注意:因为A的各阶顺序主子式hl不全大于零,所以林是正定的.A的转置矩阵zA和各阶顺序主子式值 hl依次如下:'),hl;zA, returnend endif h(1,i)>0disp('请注意:因为A的各阶顺序主子式hl都大于零,所以A是正定的.A的 转置矩阵zA和各阶顺序主子式值hl依次如下:,)hl;zAend29.0.1 /2234、 4;(2)1 1201、;(3)至1172 100;(4)C2 11-61 ,01311211341209-6飞2/000<5789;<

21、1与-6190011 了.2A0011例3.4.5判断下列矩阵是否是正定对称矩阵:(1)(1)在MATLAB作窗口输入程序>> A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9;hl=zddc (A)运行后输出结果请注意:A不是对称矩阵请注意:因为A的各阶顺序主子式 hl不全大于零,所以A不是正定的.A的转 置矩阵zA和各阶顺序主子式值 hl依次如下:zA = 1/10-1115222173-313844419hl = 1/1011/5-1601/103696/5因此,A即不是正定矩阵,也不是对称矩阵.(2)在MATLAB:作窗口输入程序>&

22、gt; A=1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19,hl=zddc(A)运行后输出结果A = 1-121-130-3209-61-3-619请注意:A是对称矩阵请注意:因为A的各阶顺序主子式hl都大于零,所以A是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下:zA = 1-121-130-3209-61-3-619hl = 12624(3)在MATLAB:作窗口输入程序>> A=1/sqrt(2) -1/sqrt(2) 0 0; -1/sqrt(2) 1/sqrt(2) 0 0; 00 1/sqrt(2) -1/sqrt(2); 0 0

23、 -1/sqrt(2) 1/sqrt(2), hl=zddc (A)运行点金加结果A= 985/1393 -985/139300-985/13939851393 -985/139300-985/1393 985/1393请注意:A是对称矩阵请注意:因为A的各阶顺序主子式 hl不全大于零,所以A不是正定的.A的转置 矩阵zA和各阶顺序主子式值 hl依次如下:zA = 985/1393 -985/139300-985/13939851393 -985/139300-985/1393985/1393hl = 985/1393000可见,A不是正定矩

24、阵,是半正定矩阵;因为 A = A T因此,A是对称矩阵.(4)在MATLAB:作窗口输入程序>> A=-2 1 1;1 -6 0;1 0 -4;hl=zddc (A)运行后输出结果A = -21130.1 -6010 -4请注意:A是对称矩阵请注意:因为A的各阶顺序主子式 hl不全大于零,所以A不是正定的.A的转置 矩阵zA和各阶顺序主子式值 hl依次如下:zA = -211 hl = -2 11 -381 -6010 -4可见A不是正定矩阵,是负定矩阵;因为 A = A T因此,A是对称矩阵.3.5求解线性方程组的LU方法及其MATLA程序3.5.1解线性方程组的楚列斯基(Ch

25、olesky )分解法及其MATLAB?序例 验证.3.5.1先将矩阵A进行楚列斯基分解,然后解矩阵方程 AX = b,并用其他方法1-121、1、-130-3,b =3209-651-3-6193在工作窗口输入>>A=1 -1 2 1;-1 3 0 -3; 2 0 9 -6;1 -3 -6 19;b1=1:2:7; b=b1' R=chol(A);C=A-R'*R,R1=inv(R);R2=R1'x=R1*R2*b,Rx=Ab-x3.5.2 解线性方程组的直接LU分解法及其MATLAB?序例3.5.2 首先将矩阵A直接进行1020、10101'b

26、=21243-110103<5)分解,然后解矩阵方程 AXLULU分解.在MATLAB作窗口输入程序解(1)首先将矩阵A直接进行运行后输出方程组的解和验证结果 x=Rx :=1.0e-014 *C=1.0e-015 *-8.0000-0.710500000.3333-0.08330-0.4441003.66670.222000002.00000.13320000>> A=10 2 0;01 0 1;12 4 3;01 0 3;b=1;2;-1;5;hl=zhjLU(A),A-L*U运行卡输出LU分解请注意:因为A的各阶主子式都不等于零, 所以A能进行LU分解.A的秩RA 和各

27、阶顺序主子式值 hl依次如下:RA = 4L = 1000U = 10200 1000 1011 21 000210 10 10002hl = 1124A分解升-个单位下三角形矩阵L和一个上三角形矩阵U的积A = LU(2)在工作窗口输入>> U=1 0 2 0;0 1 0 1;0 0 2 1;0 0 0 2; L=1 0 0 0;0 1 0 0;12 1 0;0 1 0 1;b=1;2;-1;5;U1=inv(U); L1=inv(L); X=U1*L1*b,x=Ab运行后输出方程组的解X = 8.500000000000000.50000000000000-3.75000000

28、0000001.500000000000003.5.3 解线性方程组的选主元的 LU方法及其MATLA勰序例 3.5.3先将矩阵A进彳T LU分解,然后解矩阵方程 AX =b其中解方法10.1-111I5根据(3.28 )式编写222173-3138441912-15MATLAB序,然后在工作窗口输入>> A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9;b=1;2;-1;5;L U P=LU(A),L1*P*b运行后输出结果L = 1.0000000-0.0909 1.0000000.0091 0.4628 1.000000.4545 -0.65

29、12 0.2436 1.0000U =11.0000 21.0000 13.0000 41.000003.9091-1.8182 7.7273003.7233 0.0512U1=inv(U); L1=inv(L);P = 0010010010000001X =-1.20133.3677-1.4440'X=U1*0.0536000 -4.6171方法2根据(3.29)式编写MATLA程序,然后在工作窗口输入>> A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9;b=1;2;-1;5; F U=LU(A), U1=inv(U); F1=inv(

30、F); X=U1*F1*b运行后输出结果F=0.0091 0.4628 1.00000-0.0909 1.0000001.00000000.4545 -0.6512 0.2436 1.0000X =-1.2013 3.3677 0.0536 -U=11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.7273 00 3.7233 0.0512000 -4.61711.4440'用LU分解法解线T方程组 A m冲X = b的MATLA察序function RA,RB,n,X,Y=LUjfcz(A,b)n n =size(A);B=A b;

31、RA=rank(A); RB=rank(B);for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp('请注意:因为A的r阶主子式等于零,所以A不能进行LU分解.A 的秩R解口各阶顺序主子式值hl依次如下:,)hl;RA returnendend if h(1,i)=0disp('请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:)32.X=zeros(n,1); Y=zeros(n,1); C=zeros(1,n);r=1:1;for p=1:n

32、-1max1,j=max(abs(A(p:n,p); C=A(p,:);A(p,:)= A(j+P1,:); C= A(j+P1,:);g=r(p); r(p)= r(j+P1); r(j+P1)=g;for k=p+1:nH= A(k,p)/A(p,p); A(k,p) = H; A(k,p+1:n)=A(k,p+1:n)- H* A(p,p+1:n);endendY(1)=B(r(1);for k=2:nY(k)= B(r(k)- A(k,1:k-1)* Y(1:k-1);endX(n)= Y(n)/ A(n,n);for i=n-1:-1:1X(i)= (Y(i)- A(i, i+1:n

33、) * X (i+1:n)/ A(i,i);end end RA,RB,n,X,Y';3.6 误差分析及其两种MATLA能序3.6.1 用MATLA歆件作误差分析例3.6.2 解一T7U矩阵方程 AX =b,并比较方程(1)和(2)有何区别,它们的解有何变化.其中11 / 21 / 31 /41 / 51 / 61 /7j1、1/21 /31/41 /51 / 61 /71 / 821/ 31 /41/51 /61 / 71 /81 /92(1) A =1/41 /51/61 /71 /81 /91 /10,b 二21/ 51 / 61 / 71 /81 / 91/101/1121/

34、61 / 71 / 81 /91 /101 / 111 /122/ 71 /81 / 91/101/111/121 / 13,01I .0011 / 21 /31 / 41 /51 / 61 / 71 /21 /31 /41/51 / 61 / 71 /821 / 31 / 41 /51 / 61 / 71 / 81 / 92(2)A =1 / 41 / 51 / 61 / 71 / 81 / 91 / 10b = ,2 .1 / 51 / 61 / 71 / 81 / 91 /101 / 1121 / 61 / 71 /81 / 91 /101/111 / 122J/71 / 81 / 91

35、 / 101/111 /121 /13解(1)矩阵方程 AX=b的系数矩阵A为7阶希尔伯特(Hilbert )矩阵,我们可以用下列命令计算n阶希尔伯特矩阵>>h=hilb(n)% 输出 h 为n 阶 Hilbert 矩阵在MATLABL作窗口输入程序>> A=hilb(7);b=1;2;2;2;2;2;2;X=Ab运行后输出 AX=b 的解为 X =(-35 504 -1260-4200 20790 -2772012012 )T .(2)在MATLAB:作窗口输入程序>> B =0.001,zeros(1,6);zeros(6,1),zeros(6,6);A

36、=(B+hilb(7); b=1;2;2;2;2;2;2;X=Ab运行后输出方程的解为X= (-33 465 -966 -5181 22409 -29015 12413)T.在MATLABL作窗口输入程序33.>> X =-33, 465,-966,-5181,22409,-29015,12413'X1 =-35,504,-1260,-4200,20790,-27720,12012' wu=X1'- X'运行后输出方程(1)和(2)的解的误差为X =(-2 39 -294 981 -1619 1295 - 401 )T.、一 ,、一 ,、治.001

37、0徵)万程(1)和(2)的系数矩阵的差为 A=淤,常数向量相同,则Ax = b(。6或06淘/的解的差为 整=(2 39 -294 981 -1619 1295 -401 )T . A的微小变化, 引起X的很大变化,即 X对A的扰动是敏感的.3.6.2 求P条件数和讨论AX =b解的性态的MATLAB序求P条件数和讨论AX =b解的性态的MATLA翼序function Acp=zpjxpb(A)Acw = cond (A, inf);Ac1= cond (A,1);Ac2= cond (A,2);Acf = cond (A,'fro' );dA=det(A);if (Acw&g

38、t;50)&(dA<0.1)disp('请注意:AX=b是病态的,A的松件数,1条件数,2条件数,弗罗贝 尼乌斯条件数和A的行列式的值依次如下:')Acp=Acw Ac1 Ac2 Acf dA'elsedisp( ' AX=b是良态的,A的og条件数,1条件数,2条件数,弗罗贝尼乌斯条件数和A的行列式的值依次如下:')Acp=Acw Ac1 Ac2 Acf dA'end例3.6.3根据定理3.10 ,讨论线性方程组 AX = b解的性态,并且求出 A的4种 条件数.其中5一7 .6-7MATLAB程序保存名为23-1(1) A为7阶

39、希尔伯特矩阵;(2)3 12A =4 1- 311- 24解 (1)首先将求 P条件数和讨论 AX =b解的性态的zpjxpb.m 的 四件,然后在 MATLAB:作窗口输入程序>> Acp =zpjxpb(hilb(7); Acp',det(hilb(7)运行后输出结果请注意:AX=b是病态的,A的8条件数,1条件数,2条件数,弗罗贝尼乌斯条 件数和A的行列式的值依次如下:ans = 1.0e+008 *9.8519 9.8519 4.7537 4.8175 0.0000 ans = 4.8358e-025(2)在 MATLAB:作窗口输入程序>> A=2 3

40、 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7;Acp=zpjxpb(A); Acp'运行后输出结果AX=b是良态的,A的8条件数,1条件数,2条件数,弗罗贝尼乌斯条件数和 A 的行列式的值依次如下:ans =14.1713 19.4954 8.2085 11.4203 327.00003.6.3 用P范数讨论AX = b解和A的性态的MATLAB?序 用P范数讨论 AX=b解和A的性态的MATLA翼序function Acp=zpjwc(A,jA,b,jb,P)Acp = cond (A,P);dA=det(A); X=Ab;dertaA=A-jA;34.PndA

41、=norm(dertaA, P);dertab=b-jb;Pndb=norm(dertab,P);if Pndb>0jX=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;

42、endif PndA>0jX=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,近

43、似解jX ,解的相对误差jxX ,解的相对误差估计值 Xgxx , b或A的相对误差 xAb依次如下:')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'elsedisp('请注意:AX=b是良态的,A的P条件数Acp,A的行列式值dA, 解X,近似解jX ,解的相对误差jxX ,解的相对误差估计值 Xgxx , b或A的相对误 差xAb依次如下:')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'e

44、nd例3.6.4根据定理3.10 ,讨论线性方程组 AX =b解的性态,并利用(3.32 )式讨论当A的每个元都取4位有效数字时,其解的相对误差.其中A为7阶希尔伯特矩阵, b =:11 342222 T.解(1)取畸范数和七条件数,线性方程组 AX= b的b不变时,取空范数和8条 件数,系数矩阵 A为7阶希尔伯特矩阵,A中的每个元素取4位有效数字.用P范数讨论 AX=b解和A的性态的MATLABg序保存名为zpjwc.m的文件,然后 在工作窗口输入 MATLA典序>> jA =1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.5000 0.33330.25000.20000.16670.14290.12500.3333 0.25000.20000.16670.14290.12500.11110.2

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论