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