3.000000000000000 X=
1.000000000000000 1.000000000000000
例 2.3 用Gauss?Jordan消元法思想求
?1?10?? A??223?????121??的逆阵.
解:
解法一:直接建立求解的M文件:Gauss_Jordan.m,源程序如下: %Gauss-Jordan法求例2.3 clc;
A=[1 -1 0;2 2 3;-1 2 1];
A1=A;%先保存原来的方阵A [n,m]=size(A); if n~=m
error('A必须为方阵'); return; end
A(:,n+1:2*n)=eye(n);%构造增广矩阵 for k=1:n
[l,m]=max(abs(A(k:n,k)));%按列选主元 if A(k+m-1,k)==0
error('找到列最大的元素为零,错误'); return; end
if m~=1 %交换 Temp=A(k,:);
A(k,:)=A(k+m-1,:); A(k+m-1,:)=Temp; end for i=1:n if i~=k
A(i,:)=A(i,:)-A(k,:)*A(i,k)/A(k,k); end end end
for i=n:(-1):1
A(i,:)=A(i,:)/A(i,i); end
26
A(:,1:n)=[]; disp('A='); disp(A1);
disp('用Gauss-Jandan算得矩阵A的逆矩阵为:'); disp('inv(A)='); disp(A);
clear Temp i k l m n;%清除临时变量
在Matlab命令窗口中输入Gauss_Jordan回车后得到结果如下: A=
1 -1 0 2 2 3 -1 2 1
用Gauss-Jandan算得矩阵A的逆矩阵为: inv(A)=
-4 1 -3 -5 1 -3 6 -1 4
解法二:用通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以借助
rref(C)命令求解,非常方便:
输入矩阵:
>> A=[1 -1 0;2 2 3;-1 2 1] A =
1 -1 0 2 2 3 -1 2 1
>> C=[A,eye(length(A))] C =
1 -1 0 1 0 0 2 2 3 0 1 0 -1 2 1 0 0 1 >> invA=rref(C) invA =
1 0 0 -4 1 -3 0 1 0 -5 1 -3 0 0 1 6 -1 4 后三列即为其逆阵,检验其正确性: >> c=invA(:,4:6) c =
-4 1 -3 -5 1 -3 6 -1 4
27
>> d=c*A d =
1 0 0 0 1 0 0 0 1
可见求解正确.
例 2.4 用分解LU的方法求解方程组
?2426??x1??9????x???49615???2???23??26918??x3??22???????6151840???x4??47?
解:
解线性方程组中LU分解的L,U可以实现矩阵A的三角分解,使得A=L*U. L,U应该是下三角和上三角矩阵的,这样才利于回代求根.但是MATLAB中的LU分解与解线性方程组中的L,U不一样. MATLAB的LU分解命令调用格式为:
[L,U]=lu(A)
MATLAB计算出来的L是\准下三角\交换L的行后才能成为真正的下三角阵),U为上三角矩阵,但它们还是满足A=L*U的. 先录入矩阵系数.
>> A=[2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40] A =
2 4 2 6 4 9 6 15 2 6 9 18 6 15 18 40 >> b=[9 23 22 47]' b = 9 23 22 47
将A作LU分解,方法是使用矩阵分解的LU命令即可: >>[L,D]=lu(A) L =
0.3333 1.0000 -0.6667 1.0000 0.6667 1.0000 0 0 0.3333 -1.0000 1.0000 0 1.0000 0 0 0
28
U =
6.0000 15.0000 18.0000 40.0000 0 -1.0000 -6.0000 -11.6667 0 0 -3.0000 -7.0000 0 0 0 -0.3333 再检验其正确性: >>C=L*U C =
2 4 2 6 4 9 6 15 2 6 9 18 6 15 18 40 解方程组Ly?b >>y=L\\b y =
47.0000 -8.3333 -2.0000 0.3333
解方程组Ux?y得到方程组的最终解: >>x=U\\y x =
0.5000 2.0000 3.0000 -1.0000
故方程组的最终解为: x?(0.5,2,3,?1)T
例 2.5 解方程组
??610???141??x1???x???6??2???24? ??0114????x3????322??试用平方根法,改进的平方根法和追赶法分别解之.
解:
先输入矩阵:
>>A=[6 1 0;1 4 1;0 1 14] A =
6 1 0 1 4 1 0 1 14
29
>>b=[6 24 322]' b = 6 24 322 平方根法:
先对A矩阵作Cholesky分解:
>>L=chol(A) L =
2.4495 0.4082 0 0 1.9579 0.5108 0 0 3.7066 检验其正确性 >>L'*L ans =
6.0000 1.0000 0 1.0000 4.0000 1.0000 0 1.0000 14.0000 将L转化为下三角矩阵: >>L=L' L =
2.4495 0 0 0.4082 1.9579 0 0 0.5108 3.7066 解方程组Ly?b: >>y=L\\b y =
2.4495 11.7473 85.2526
再解方程组LTx?y,得到最终解: >>x=L'\\y x =
1.0000 -0.0000 23.0000
故平方根法解上述方程的结果为: x?(1,0,23)T改进的平方根法:
先对矩阵A作LDL分解: >> [L,D]=ldl(A)
30