第二章 解线性方程组的直接方法 matlab用法范文

第二章 解线性方程组的直接方法

在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法.

2.1 方程组的逆矩阵解法及其MATLAB程序

2.1.3 线性方程组有解的判定条件及其MATLAB程序 判定线性方程组Am?nX?b是否有解的MATLAB程序

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，所以此方程

组无解.')

return end

if RA==RB if RA==n

disp('请注意：因为RA=RB=n，所以此方

程组有唯一解.')

else

disp('请注意：因为RA=RB

程组有无穷多解.')

end end

例2.1.4 判断下列线性方程组解的情况.如果有唯一解，则用表 3-2方法求解.

(1)

?2x1?3x2?x3?5x4?0,?3x?x?2x?7x?0,?1234??4x1?x2?3x3?6x4?0,??x1?2x2?4x3?7x4?0; (2)

第二章 解线性方程组的直接方法的MATLAB程序

?3x1?4x2?5x3?7x4?0,?2x?3x?3x?2x?0,?1234 ??4x1?11x2?13x3?16x4?0,??7x1?2x2?x3?3x4?0;?4x1?2x2?x3?2,?2x?y?z?w?1,?(3) ?3x1?x2?2x3?10, (4) ??4x?2y?2z?w?2, ?11x?3x?8;?2x?y?z?w?1.?12?

解 在MATLAB工作窗口输入程序

>> 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=A\\b,

运行后输出结果为 X =(0 0 0 0)’.

(2) 在MATLAB工作窗口输入程序

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

运行后输出结果

请注意：因为RA=RB

多解.

RA =2,RB =2,n =4

(3) 在MATLAB工作窗口输入程序

>> A=[4 2 -1;3 -1 2;11 3 0]; b=[2;10;8];

[RA,RB,n]=jiepb(A,B) 运行后输出结果

请注意：因为RA~=RB，所以此方程组无解. RA =2,RB =3,n =3

25．

(4)在MATLAB工作窗口输入程序

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

多解.

RA =2,RB =2,n =3

2.2 三角形方程组的解法及其MATLAB程序

2.2.2 解三角形方程组的MATLAB程序 解上三角形线性方程组AX?b的MATLAB程序

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('请注意：因为RA~=RB，所以此方程

组无解.')

return end

if RA==RB if RA==n

disp('请注意：因为RA=RB=n，所以此

方程组有唯一解.')

X=zeros(n,1);

X(n)=b(n)/A(n,n);

for k=n-1:-1:1

X(k)=(b(k)-sum(A(k,k+1:n)*X(k+1:n)))/A(k,k);

end else

disp('请注意：因为RA=RB

方程组有无穷多解.')

end end

第二章 解线性方程组的直接方法的MATLAB程序

例2.2.2 用解上三角形线性方程组的MATLAB程序解方程组

?5x1?x2?2x3?3x4?20,??2x?7x?4x??7,?234. ?6x3?5x4?4,??3x4?6.?解 在MATLAB工作窗口输入程序

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

4,

X =[2.4 -4.0 -1.0 2.0]’

2.3 高斯（Gauss）消元法和列主元消元法及其MATLAB程序

2.3.1 高斯消元法及其MATLAB程序

用高斯消元法解线性方程组AX?b的MATLAB程序

function [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，所以此方程组无解.') return end

if RA==RB if RA==n

disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1

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

b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1

X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end

else

27．

disp('请注意：因为RA=RB

例2.3.2 用高斯消元法和MATLAB程序求解下面的非齐次线性方程组，并且用逆矩阵解方程组的方法验证.

?x1?x2?x3?3x4?1,??x2?x3?x4?0, ???2x1?2x2?4x3?6x4??1,??x1?2x2?4x3?x4??1.解 在MATLAB工作窗口输入程序

>> 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=RB=n，所以此方程组有唯一解. RA =

X = 4

0 RB =

-0.5000 4

0.5000 n =

0 4

2.3.2 列主元消元法及其MATLAB程序

用列主元消元法解线性方程组AX?b的MATLAB程序

function [RA,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，所以此方程组无解.') return end

if RA==RB if RA==n

disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1

[Y,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); end end

b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1

X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end

else

disp('请注意：因为RA=RB

例2.3.3 用列主元消元法解线性方程组的MATLAB程序解方程组

第二章 解线性方程组的直接方法的MATLAB程序

?x2?x3?x4?0,??x?x?x?3x?1,?1234. ??2x1?2x2?4x3?6x4??1,??x1?2x2?4x3?x4??1.解 在MATLAB工作窗口输入程序

>> A=[0 -1 -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]’

2.4 LU分解法及其MATLAB程序

2.4.1判断矩阵LU分解的充要条件及其MATLAB程序 判断矩阵A能否进行LU分解的MATLAB程序

function hl=pdLUfj(A)

[n n] =size(A); RA=rank(A); if RA~=n

disp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A

的秩RA如下：'), RA,hl=det(A); return

end

if RA==n

for p=1:n,h(p)=det(A(1:p, 1:p));, end

hl=h(1:n); for i=1:n

if h(1,i)==0

disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分

解.A的秩RA和各阶顺序主子式值hl依次如下：'),hl;RA,return

end end

if h(1,i)~=0

disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分

解.A的秩RA和各阶顺序主子式值hl依次如下：')

hl;RA end end

例2.4.1 判断下列矩阵能否进行LU分解，并求矩阵的秩.

?123??123??123???????（1）?1127?；（2）?127?；（3）?123?.

?456??456??456???????解 （1）在MATLAB工作窗口输入程序

>> A=[1 2 3;1 12 7;4 5 6];hl=pdLUfj(A)

运行后输出结果为

请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和

各阶顺序主子式值hl依次如下：

RA = 3， hl = 1 10 -48

（2）在MATLAB工作窗口输入程序

>> A=[1 2 3;1 2 7;4 5 6];hl=pdLUfj(A)

运行后输出结果为

请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶

顺序主子式值hl依次如下：

29．

RA = 3， hl =1 0 12

（3）在MATLAB工作窗口输入程序

>> A=[1 2 3;1 2 3;4 5 6];hl=pdLUfj(A)

运行后输出结果为

请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下

RA = 2， hl = 0

2.4.2 直接LU分解法及其MATLAB程序 将矩阵A进行直接LU分解的MATLAB程序

function hl=zhjLU(A)

[n n] =size(A); RA=rank(A); if RA~=n

disp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A

的秩RA如下:'), RA,hl=det(A);

return end

if RA==n for p=1:n

h(p)=det(A(1:p, 1:p)); end

hl=h(1:n); for i=1:n

if h(1,i)==0

disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A

的秩RA和各阶顺序主子式值hl依次如下：'), hl;RA

return end end

if h(1,i)~=0

disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A

的秩RA和各阶顺序主子式值hl依次如下：')

for j=1:n

U(1,j)=A(1,j); end

for 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 end end

hl;RA,U,L end end

例2.4.3 用矩阵进行直接LU分解的MATLAB程序分解矩阵

第二章 解线性方程组的直接方法的MATLAB程序

?1020????0101?. A??1243????0103???解 在MATLAB工作窗口输入程序

>> 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依次如下：

L = 1 0 0 0 RA = 4

0 1 0 0 U = 1 0 2 0

1 2 1 0 0 1 0 1

0 1 0 1 0 0 2 1

0 0 0 2 hl = 1 1 2 4

2.4.4 判断正定对称矩阵的方法及其MATLAB程序 判断矩阵A是否是正定对称矩阵的MATLAB程序

function hl=zddc(A) [n n] =size(A); for p=1:n

h(p)=det(A(1:p, 1:p)); end

hl=h(1:n);zA=A'; for i=1:n

if h(1,i)<=0

disp('请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正

定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：'), hl;zA,return

end end

if h(1,i)>0

disp('请注意：因为A的各阶顺序主子式hl都大于零，所以A是正定的.A的

转置矩阵zA和各阶顺序主子式值hl依次如下：')

hl;zA end

例2.4.5 判断下列矩阵是否是正定对称矩阵：

?0.1234??1?121?1?1???211??????00???2?；（1）??12?34?；(2) ??130?3?； (3) ?2（4）?1?60??. 1?1??11211341??209?6?00????10?4?22?????????5789??1?3?619?11??0????0?????0?0?2122?1???2?解 （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 -1 11 5 2 2 21 7 3 -3 13 8 4 4 41 9 hl = 1/10 11/5 -1601/10 3696/5

因此,A即不是正定矩阵，也不是对称矩阵.

（2）在MATLAB工作窗口输入程序

>> A=[1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19],hl=zddc(A)

31．

运行后输出结果

A = 1 -1 2 1 -1 3 0 -3 2 0 9 -6 1 -3 -6 19 请注意： A是对称矩阵

请注意：因为A的各阶顺序主子式hl都大于零，所以A是正定的.A的转置矩阵zA

和各阶顺序主子式值hl依次如下：

zA = 1 -1 2 1 -1 3 0 -3 2 0 9 -6 1 -3 -6 19 hl = 1 2 6 24 （3）在MATLAB工作窗口输入程序

>> A=[1/sqrt(2) -1/sqrt(2) 0 0; -1/sqrt(2) 1/sqrt(2) 0 0; 0

0 1/sqrt(2) -1/sqrt(2); 0 0 -1/sqrt(2) 1/sqrt(2)], hl=zddc(A) 运行后输出结果

A= 985/1393 -985/1393 0 0 -985/1393 985/1393 0 0 0 0 985/1393 -985/1393 0 0 -985/1393 985/1393 请注意： A是对称矩阵

请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置

矩阵zA和各阶顺序主子式值hl依次如下：

zA = 985/1393 -985/1393 0 0 -985/1393 985/1393 0 0 0 0 985/1393 -985/1393 0 0 -985/1393 985/1393 hl = 985/1393 0 0 0

T

可见，A不是正定矩阵，是半正定矩阵；因为A= A 因此，A是对称矩阵.

（4）在MATLAB工作窗口输入程序

>> A=[-2 1 1;1 -6 0;1 0 -4];hl=zddc(A)

运行后输出结果

A = -2 1 1

1 -6 0 1 0 -4 请注意： A是对称矩阵

请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置

矩阵zA和各阶顺序主子式值hl依次如下：

zA = -2 1 1 hl = -2 11 -38 1 -6 0 1 0 -4

T

可见A不是正定矩阵，是负定矩阵；因为A = A 因此，A是对称矩阵.

2.5 求解线性方程组的LU方法及其MATLAB程序

2.5.1 解线性方程组的楚列斯基（Cholesky）分解法及其MATLAB程序

例3.5.1 先将矩阵A进行楚列斯基分解,然后解矩阵方程AX?b，并用其他方法验证.

1??1?12?1??????130?3???3?. A??,b??5?209?6??????1?3?619??7?????第二章 解线性方程组的直接方法的MATLAB程序

解 在工作窗口输入

>>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=A\\b-x

运行后输出方程组的解和验证结果

x = Rx = 1.0e-014 * C = 1.0e-015 * -8.0000 -0.7105 0 0 0 0 0.3333 -0.0833 0 -0.4441 0 0 3.6667 0.2220 0 0 0 0 2.0000 0.1332 0 0 0 0

2.5.2 解线性方程组的直接LU分解法及其MATLAB程序

例3.5.2 首先将矩阵A直接进行LU分解，然后解矩阵方程AX?b

?1??1020?????2?. ?0101?，b??A????1?1243??????5??0103?????解 (1) 首先将矩阵A直接进行LU分解.在MATLAB工作窗口输入程序

>> A=[1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3];b=[1;2;-1;5];

hl=zhjLU(A),A-L*U 运行后输出LU分解

请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA

和各阶顺序主子式值hl依次如下：

L = 1 0 0 0 RA = 4

0 1 0 0 U = 1 0 2 0

1 2 1 0 0 1 0 1

0 1 0 1 0 0 2 1

hl = 1 1 2 4 0 0 0 2

A分解为一个单位下三角形矩阵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;1

2 1 0;0 1 0 1];

b=[1;2;-1;5];U1=inv(U); L1=inv(L); X=U1*L1*b,x=A\\b

运行后输出方程组的解

X = 8.000 0.000 -3.000 1.000

2.5.3 解线性方程组的选主元的LU方法及其MATLAB程序

例3.5.3 先将矩阵A进行LU分解，然后解矩阵方程AX?b 其中

34??0.12?1??????12?34??，?2?. A??b???1?11211341??????5??5?789????解 方法1 根据（3.28）式编写MATLAB程序，然后在工作窗口输入

>> 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), U1=inv(U); L1=inv(L); X=U1*

L1*P*b

运行后输出结果

L = 1.0000 0 0 0 P = 0 0 1 0

0 1 0 0

1 0 0 0 33．

0 0 0 1

X =[-1.2013 3.3677 0.0536

-1.4440]’

-0.0909 1.0000 0 0 0.0091 0.4628 1.0000 0 0.4545 -0.6512 0.2436 1.0000 U =11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.7273 0 0 3.7233 0.0512 0 0 0 -4.6171

方法2 根据（3.29）式编写MATLAB程序，然后在工作窗口输入

>> 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(F); X=U1*F1*b

运行后输出结果

F=0.0091 0.4628 1.0000 0 U=11.0000 21.0000 13.0000 41.0000 -0.0909 1.0000 0 0 0 3.9091 -1.8182 7.7273 1.0000 0 0 0 0 0 3.7233 0.0512 0.4545 -0.6512 0.2436 1.0000 0 0 0 -4.6171 X =[-1.2013 3.3677 0.0536 -1.4440]’ 用LU分解法解线性方程组Am?nX?b的MATLAB程序

function [RA,RB,n,X,Y]=LUjfcz(A,b)

[n n] =size(A);B=[A b]; RA=rank(A); RB=rank(B); for p=1:n

h(p)=det(A(1:p, 1:p)); end

hl=h(1:n); for i=1:n

if h(1,i)==0

disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A

的秩RA和各阶顺序主子式值hl依次如下：')

hl;RA return end end

if h(1,i)~=0

disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分

解.A的秩RA和各阶顺序主子式值hl依次如下：')

X=zeros(n,1); Y=zeros(n,1); C=zeros(1,n);r=1:1; for p=1:n-1

[max1,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:n

H= 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);

end end

Y(1)=B(r(1)); for k=2:n

Y(k)= B(r(k))- A(k,1:k-1)* Y(1:k-1); end

X(n)= Y(n)/ A(n,n); for i=n-1:-1:1

X(i)= (Y(i)- A(i, i+1:n) * X (i+1:n))/ A(i,i); end end

[RA,RB,n,X,Y]’;

2.6 误差分析及其两种MATLAB程序

第二章 解线性方程组的直接方法的MATLAB程序

2.6.1 用MATLAB软件作误差分析

例2.6.2 解下列矩阵方程AX?b，并比较方程（1）和（2）有何区别，它们的解有何变化.其中

1/21/31/41/51/61/7?1??1?????1/31/41/51/61/71/8?2??1/2??2??1/3?1/41/51/61/71/81/9?? ??A??1/41/51/61/71/81/91/10?,b??2?;?2??1/5?1/61/71/81/91/101/11?????2??1/61/71/81/91/101/111/12?????1/71/81/91/101/111/121/13?2????1?1/7?1.0011/21/31/41/51/6?????1/31/41/51/61/71/8??2??1/2?2??1/31/41/51/61/71/81/9??? ??A??1/41/51/61/71/81/91/10?,b??2?.?2??1/51/61/71/81/91/101/11??????2??1/61/71/81/91/101/111/12??????2??1/71/81/91/101/111/121/13?(1)(2)解 (1) 矩阵方程AX?b的系数矩阵A为7阶希尔伯特（Hilbert）矩阵，我们可

以用下列命令计算n阶希尔伯特矩阵

>>h=hilb(n) % 输出h为n阶Hilbert矩阵 在MATLAB工作窗口输入程序

>> A=hilb(7);b=[1;2;2;2;2;2;2];X=A\\b 运行后输出AX?b的解为 X =（-35 504 -1260 -4200 20790 -27720 12012)T.

（2）在MATLAB工作窗口输入程序

>> B =[0.001,zeros(1,6);zeros(6,1),zeros(6,6)]; A=(B+hilb(7)); b=[1;2;2;2;2;2;2];X=A\\b

T运行后输出方程的解为 X=（-33 465 -966 -5181 22409 -29015 12413).

在MATLAB工作窗口输入程序

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

方程（1）和（2）的系数矩阵的差为?A???的解的差为?X?0.001O1?6??，常数向量相同，则Ax?b??O6?1O6?6?T?(?239?294981?16191295?401).A的微小变化，

引起X的很大变化，即X对A的扰动是敏感的.

2.6.2 求P 条件数和讨论AX?b解的性态的MATLAB程序 求P条件数和讨论AX?b解的性态的MATLAB程序

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>50)&(dA<0.1)

disp('请注意：AX=b是病态的，A的∞条件数,1条件数,2条件数, 弗罗贝尼乌斯条件数和A的行列式的值依次如下：') Acp=[Acw Ac1 Ac2 Acf dA]'; else

disp(' AX=b是良态的，A的∞条件数,1条件数,2条件数,弗罗贝尼乌斯条

35．

件数和A的行列式的值依次如下：') Acp=[Acw Ac1 Ac2 Acf dA]'; end

例2.6.3 根据定理3.10，讨论线性方程组AX?b解的性态，并且求出A的4种条件数.其中

?2?（1）A为7阶希尔伯特矩阵；（2）?3A??4??1?311?2?12?345???7?. 6???7??解 （1）首先将求P条件数和讨论AX?b解的性态的MATLAB程序保存名为

zpjxpb.m 的M文件，然后在MATLAB工作窗口输入程序

>> Acp =zpjxpb(hilb(7)); Acp',det(hilb(7))

运行后输出结果

请注意：AX=b是病态的，A的∞条件数,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 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7];Acp=zpjxpb(A); Acp' 运行后输出结果

AX=b是良态的，A的∞条件数,1条件数,2条件数, 弗罗贝尼乌斯条件数和A

的行列式的值依次如下：

ans =

14.1713 19.4954 8.2085 11.4203 327.0000

2.6.3 用P范数讨论AX?b解和A的性态的MATLAB程序 用P范数讨论AX?b解和A的性态的MATLAB程序

function Acp=zpjwc(A,jA,b,jb,P)

Acp = cond (A,P);dA=det(A); X=A\\b; dertaA=A-jA;

PndA=norm(dertaA, P);dertab=b-jb;Pndb=norm(dertab, P); if Pndb>0

jX=A\\jb; 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=jA\\b; 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; end

if (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'

第二章 解线性方程组的直接方法的MATLAB程序

else

disp('请注意： AX=b是良态的，A的P条件数Acp,A的行列式值dA，

解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：')

Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj' end

例2.6.4 根据定理3.10，讨论线性方程组AX?b解的性态，并利用（3.32）式讨论当A的每个元都取4位有效数字时，其解的相对误差.其中A为7阶希尔伯特矩阵，

b?1?1342222?T.

解 （1）取?范数和?条件数，线性方程组AX?b的b不变时，取?范数和?条件数，系数矩阵A为7阶希尔伯特矩阵，A中的每个元素取4位有效数字.

用P范数讨论AX?b解和A的性态的MATLAB程序保存名为zpjwc.m的文件,然后在工作窗口输入MATLAB程序

>> jA =[1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429

0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769]; A=hilb(7);b=[1;1/3;4;2;2;2;2];

jb=[1;0.3333;4;2;2;2;2]; Acp=zpjwc(A,jA,b,jb,inf)

运行后输出结果

请注意：AX=b是病态的，A的P 条件数Acp,A的行列式值dA，解X，近似

解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：

Acp = dA =

9.8519e+008 4.8358e-025 ans =

1.0e+007 *

0.0020 -0.0697 0.6243 -2.3054 4.0677 -3.4123

1.0943

ans =

1.0e+004 *

0.0349 -0.4807 2.1126 -5.1087 7.6557 -6.3239

2.1112

xX = jxX = Xgxx =

0.9981 530.3248 4.9291e+004 xAb = xAbj =

5.0032e-005 5.0031e-005

由此可见，因为?条件数Cond(A)??985 194 889.719 848 3??1，所以此方程组为病态的AX?b的解 X ?(19565, -697256, 6243300,-2.305380,40676790,-34123320,10942932)T， (jA)X?b的解为

jX ?(349, -48079, 21126,-51087,76557,-63239,21112)T 解的相对误差

?XX???0.99812,

?X??X??X?49291.27 ,

?X?X??X?530.32?,

?AA???5.0032?10?5,

即相对误差放大了约985 194 889.72倍.

(2) 如果取2范数和2条件数计算，在MATLAB工作窗口输入程序

37．

>> Acp =zpjwc(A,jA,b,jb,2)

运行后输出结果

请注意：AX=b是病态的，A的P 条件数Acp,A的行列式值dA，解X，近似

解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：

Acp = dA =

4.7537e+008 4.8358e-025 ans =

1.0e+007 *

0.0020 -0.0697 0.6243 -2.3054 4.0677 -3.4123

1.0943

ans =

1.0e+004 *

0.0349 -0.4807 2.1126 -5.1087 7.6557 -6.3239

2.1112

xX = jxX = Xgxx=

0.9981 511.0640 2.9951e+004 xAb = xAbj = 6.3006e-005 6.3005e-005

因为2条件数Cond(A)2? 475 367 356.59??1，所以此方程组为病态的.解的相对误差

?X2?A2?X2?A2? 0.9981,?29950.85. ?6.3005?10?5,?511.06,X2A??A2A2X??X2即相对误差放大了约475 367 356.59倍.

本文来源：https://www.bwwdw.com/article/6gx5.html