2 线性方程组的直接解法

2.1  例题解答

例 2.1 用Gauss消元法解方程组:

解:

直接建立求解该方程组的M文件Gauss.m如下:

 

%求解例题2.1

%高斯法求解线性方程组Ax=b

%A为输入矩阵系数,b为方程组右端系数

%方程组的解保存在x变量中

%先输入方程系数

A=[1 2 3;2 7 5;1 49];

b=[1 6 -3]';

[m,n]=size(A);

%检查系数正确性

if m~=n

    error('矩阵A的行数和列数必须相同');

    return;

end

if m~=size(b)

    error('b的大小必须和A的行数或A的列数相同');

    return;

end

%再检查方程是否存在唯一解

if rank(A)~=rank([A,b])

    error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');

    return;

end

%这里采用增广矩阵行变换的方式求解

c=n+1;

A(:,c)=b;

%%消元过程

for k=1:n-1

A(k+1:n,k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c);

End

%%回代结果

x=zeros(length(b),1); 

x(n)=A(n,c)/A(n,n);

for k=n-1:-1:1

x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end

%显示计算结果

disp('x=');

disp(x);

 

直接运行上面的M文件或在Matlab命令窗口中直接输入Gauss即可得出结果.

在Matlab命令窗口中输入Gauss得出结果如下:

>> Gauss

x=

    2.0000

    1.0000

   -1.0000

 

扩展:

Matlab求解线性方程的几种命令如下(方程组的一般形式可用矩阵和向量表示成: ,但运用下列方法的前提必须保证所求解的方程为恰定方程,即方程组存在唯一的一组解):

运用求逆思想 :  或 ;

左除法,原理上是运用高斯消元法求解,但Matlab在实际执行过程中是通过 分解法进行的(即先将矩阵A作 分解,再回代计算): ;

用符号矩阵法计算,这种计算方法最接近精确值,但计算速度最慢:

通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以这样实现之: ; .

上面四种常用的办法示例如下:

>> A=[1 23;2 7 5;1 4 9]  %上面示例方程组系数

A =

     1    2     3

     2    7     5

     1    4     9

>> b=[1 6-3]'   %方程组右端的系数

b =

     1

     6

-3

>>x1_1=inv(A)*b,x1_2=A^(-1)*b  %方法一,求逆思想

x1_1 =

    2.0000

    1.0000

   -1.0000

x1_2 =

    2.0000

    1.0000

   -1.0000

>> x2=A\b  %方法二,左除思想

x2 =

     2

     1

-1

>>x3=sym(A)\sym(b) %方法三,符号法

 x3 =

   2

   1

  -1

>>C=[A,b],rref(C)  %方法四, 行简化阶梯形思想,最后输出结果的一列为解

C =

 

     1    2     3     1

     2    7     5     6

     1    4     9    -3

ans =

 

     1    0     0     2

     0    1     0     1

     0    0     1    -1

 

例 2.2 用选列主元的Gauss消元法解如下方程:

解:

直接建立求解该方程的M文件Gauss_line.m,求解程序编制如下:

 

%求解例题2.2

%高斯列主元消元法求解线性方程组Ax=b

%A为输入矩阵系数,b为方程组右端系数

%方程组的解保存在x变量中

format long;%设置为长格式显示,显示15位小数

A=[0.00001,1;2,1]

b=[1.00001,3]'

[m,n]=size(A);

%先检查系数正确性

if m~=n

    error('矩阵A的行数和列数必须相同');

    return;

end

if m~=size(b)

    error('b的大小必须和A的行数或A的列数相同');

    return;

end

%再检查方程是否存在唯一解

ifrank(A)~=rank([A,b])

    error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');

    return;

end

c=n+1;

A(:,c)=b; %(增广)

for k=1:n-1

[r,m]=max(abs(A(k:n,k)));   %选主元

 m=m+k-1; %修正操作行的值

   if(A(m,k)~=0)

         if(m~=k)

             A([k m],:)=A([m k],:); %换行

         end

         A(k+1:n, k:c)=A(k+1:n,k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c);  %消去

   end

end

x=zeros(length(b),1);%回代求解

x(n)=A(n,c)/A(n,n);

for k=n-1:-1:1

    x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end

disp('X=');

disp(x);

format short;%设置为默认格式显示,显示5位

 

运行得到输出结果如下:

A =

   0.000010000000000   1.000000000000000

   2.000000000000000   1.000000000000000

b =

   1.000010000000000

   3.000000000000000

X=

   1.000000000000000

   1.000000000000000

 

例 2.3 用 消元法思想求

的逆阵.

解:

解法一:直接建立求解的M文件:Gauss_Jordan.m,源程序如下:

%Gauss-Jordan法求例2.3

clc;

A=[1 -1 0;2 2 3;-12 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

A(:,1:n)=[];

disp('A=');

disp(A1);

disp('用Gauss-Jandan算得矩阵A的逆矩阵为:');

disp('inv(A)=');

disp(A); 

clear Temp i k l mn;%清除临时变量

 

在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

 

解法二:用通过将矩阵施行初等行变换化成行简化阶梯形的办法,可以借助 命令求解,非常方便:

输入矩阵:

>> A=[1 -10;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

>> d=c*A

d =

 

     1    0     0

     0    1     0

     0    0     1

 

可见求解正确.

 

2.4 用分解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 26;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 2322 47]'

b =

     9

    23

    22

    47

将A作 分解,方法是使用矩阵分解的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

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

解方程组

>>y=L\b  

y =

  47.0000

  -8.3333

  -2.0000

   0.3333  

解方程组 得到方程组的最终解:

>>x=U\y  

x =

   0.5000

   2.0000

   3.0000

  -1.0000  

故方程组的最终解为:

 

例 2.5 解方程组

试用平方根法,改进的平方根法和追赶法分别解之.

解:

先输入矩阵:

>>A=[6 1 0;1 4 1;0 1 14]  

A =

    6     1     0

    1     4     1

    0     1    14  

>>b=[6 24 322]'  

b =

    6

   24

  322

平方根法:

先对A矩阵作 分解:

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

解方程组 :

>>y=L\b  

y =

   2.4495

  11.7473

  85.2526  

再解方程组 ,得到最终解:

>>x=L'\y  

x =

    1.0000

  -0.0000

  23.0000  

故平方根法解上述方程的结果为:

改进的平方根法:

先对矩阵A作LDL分解:

>> [L,D]=ldl(A)

L =

   1.0000         0         0

   0.1667    1.0000         0

        0    0.2609    1.0000

D =

   6.0000         0         0

        0    3.8333        0

        0         0   13.7391  

检验其分解正确性:

>>L*D*L'  

ans =

    6     1     0

    1     4     1

    0     1    14  

解方程组 :

>> y=L\b  

y =

    6

   23

  316  

解方程组 :

>> x=(D*L')\y 

x =

    1

    0

   23  

故改进的平方根法解上述方程的结果亦为:

追赶法:

编制追赶法求解该方程的程序如下:

%pursue.m

%三对角线性方程组的追赶法解方程组例2.5

%输入矩阵

clc;

A=[6 1 0;1 4 1;0 114]

f=[6 24 322]

[n,m]=size(A);

%分别取对角元素

a=zeros(1,n);

a(2:n)=diag(A,-1);

c=diag(A,1);

%此处用变量d存储A主对角线上的元素,因已用变量b存储方程右边的系数

b=diag(A);

if b(1)==0

    error('主对角元素不能为0');

    return;

end

%初始计算,式(2.31)

alpha(1)=b(1);

beta(1)=c(1)/b(1);

%按照公式(2.31)计算

for i=2:n-1

    alpha(i)=b(i)-a(i)*beta(i-1);

    if alpha(i)==0

        error('错误:在解方程过程中α为0');

        return;

    end

    beta(i)=c(i)/alpha(i);

end

%对最后一行作计算

alpha(n)=b(n)-a(n)*beta(n-1);

if alpha(n)==0

        error('错误:在解方程过程中最后一个α为0');

        return;

end

%以下按照公式(2.32)计算,解Ly=f

y(1)=f(1)/b(1);

for i=2:n

    y(i)=(f(i)-a(i)*y(i-1))/alpha(i);

end

%以下按照公式(2.33)计算,解Ux=y

X(n)=y(n);

for i=n-1:-1:1

    X(i)=y(i)-beta(i)*X(i+1);

end

disp('X=');

disp(X);

 

在Matlab命令窗口输入pursue计算结果如下:

>>pursue  

A =

    6     1     0

    1     4     1

    0     1    14

f =

    6    24   322

X=

    1     0    23  

其中A为系数矩阵,f为矩阵右端的系数,最后计算结果为X.

由以上计算可知追赶法解该方程的结果亦为:

 

例 2.6 ,求 .

解:

输入矩阵:

>>x=[1 0.5 -0.3]  

x =

   1.0000    0.5000   -0.3000  

计算其1-范数

>>norm_1=norm(x,1)  

norm_1 =

   1.8000  

计算其2-范数

>>norm_2=norm(x)  

norm_2 =

   1.1576  

计算其无穷大范数

>>norm_inf=norm(x,inf)  

norm_inf =

    1  

还可以计算其无穷小范数(即各分量绝对值中的最小值)

>>norm_minusInf=norm(x,-inf)  

norm_minusInf =

   0.3000  

 

例 2.7

求 .

解:

先输入矩阵:

>>A=[-2 1 0;1 -2 1;0 1 -2]  

A =

   -2     1     0

    1    -2     1

    0     1    -2  

A的1-范数(列和范数):

>>norm_1=norm(A,1)  

norm_1 =

    4  

求解A的无穷大范数(行和范数):

>>norm_inf=norm(A,inf)  

norm_inf =

    4  

A的2-范数( 最大特征值):

>>norm_2=norm(A,2)  

norm_2 =

   3.4142  

还可以求解A的F-范数:

>>norm_F=norm(A,'fro')  

norm_F =

    4  

谱半径可以通过按其概念进行计算:对其特征值的绝对值取最大值即可:

>>R_A=max(abs(eig(A)))  

 

R_A =

   3.4142  

 

例 2.8  n阶Hilbert矩阵

考查其条件数.

解:

上述矩阵为抽象矩阵,而计算机只能进行有限次迭代,我们从 考查其条件数,为此编制如下程序Hilb_cond10.m:

 

% Hilb_cond10.m

%考查从3阶至10的 矩阵2—条件数

for i=3:10

   h=hilb(i);

   condA(i)=cond(h,2);

end

disp('n           cond');

for i=3:10   

   s=sprintf('%d         %f',i,condA(i));

   disp(s);

end

 

运行后得到如下结果:

n           cond

3        524.056778

4        15513.738739

5        476607.250242

6        14951058.641005

7        475367356.370393

8        15257575270.772364

9        493153986466.270940

10       16025391750078.617000

 

结果中左边为 的阶数,右边为对应的条件数 ,从以上结果也可分析可知:当n的阶数稍高时, 矩阵即出现严重的病态

.

例 2.9 求

的条件数

解:

建立2阶 矩阵:

>>H=hilb(2)  

H =

   1.0000    0.5000

   0.5000    0.3333  

求其2-条件数

>>cond_2=cond(H,2)  

cond_2 =

  19.2815  

求其1-条件数

>>cond_1=cond(H,1)  

cond_1 =

  27.0000  

求其无穷大条件数

>>cond_inf=cond(H,inf)  

cond_inf =

  27.0000  

还可求其F-条件数

>>cond_f=cond(H,'fro')  

cond_f =

  19.3333  

阅读更多
想对作者说点什么?

博主推荐

换一批

没有更多推荐了,返回首页