利用matlab解线性方程组

数值计算实验

——解线性方程组

西南交通大学

2012级茅7班

20123257 陈鼎

摘要

本报告主要介绍了基于求解线性方程组的高斯消元法和列主消元法两种数值分析方法的算法原理及实现方法。运用matlab数学软件辅助求解。

实验内容

1．编写用高斯消元法解线性方程组的MATLAB程序，并求解下面的线性方程组，然后用逆矩阵解方程组的方法验证。

2．编写用列主消元法解线性方程组的MATLAB程序，并求解下面的线性方程组，然后用逆矩阵解方程组的方法验证。

给定方程组如下：

①0.325x1+2.564x2+3.888x3+5x4=1.521

②-1.548x1+3.648x2+4.214x3-4.214x4=2.614

③-2.154x1+1.647x2+5.364x3+x4=3.978

④0x1+2.141x2-2.354x3-2x4=4.214

A.高斯消元法

一、算法介绍

高斯消元法是一种规则化的加减消元法。基本思想是通过逐次消元计算把需要求解的线性方程组转化成为上三角方程组，即把现形方程组的系数矩阵转化为上三角矩阵，从而使一般线性方程组的求解转化为等价的上三角方程组的求解。

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

disp(‘因为RA=RB

end

end

三、实验过程与结果

输入的量：系数矩阵A和常系数向量b；

输出的量：系数矩阵A和增广矩阵B的秩RA、RB，方程中未知量的个数n和有关方程组解X及其解的信息。

1.编写高斯消元法的运行MATLAB代码如下：

tic,

A=[0.325 2.564 3.888 5;-1.548 3.648 4.214 -4.214;-2.154 1.647 5.364 1;0 2.141 -2.354 -2];

b=[1.521 ;2.614 ;3.978 ;4.214] ;

[RA,RB,n,X] =gaus (A,b),toc

运行结果为：

因为RA=RB=n，所以此方程组有唯一解.

RA =

4

RB =

4

n =

4

X =

-3.2161

1.3785

-1.0959

0.6585

Elapsed time is 0.004000 seconds.

2.编写逆矩阵检验的运行MATLAB代码如下：

tic,

A=[0.325 2.564 3.888 5;-1.548 3.648 4.214 -4.214;-2.154 1.647 5.364 1;0 2.141 -2.354 -2];

b=[1.521 ;2.614 ;3.978 ;4.214] ;

X=A\b,toc

运行结果为：

X =

-3.2161

1.3785

-1.0959

0.6585

Elapsed time is 0.002000 seconds.

与步骤1中运算结果相同，检验成功。

B.列主消元法

一、算法介绍

列主消元法是为控制舍入误差而提出来的一种算法，在高斯消元的过程中，若出现a=0则消元无法进行，即使其不为0，但很小，把他作为除数，就会导致其他元素量级的巨大增长和舍入误差的扩散，最后使计算结果不可靠，使用列主消元法计算，基本上能控制舍入误差的影响，并且选主元素比较方便。

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

end

end

三、实验过程与结果

输入的量：系数矩阵A和常系数向量b；

输出的量：系数矩阵A和增广矩阵B的秩RA、RB，方程中未知量的个数n和有关方程组解X及其解的信息。

1.编写列主消元法的运行MATLAB代码如下：

tic,

A=[0.325 2.564 3.888 5;-1.548 3.648 4.214 -4.214;-2.154 1.647 5.364 1;0 2.141 -2.354 -2];

b=[1.521 ;2.614 ;3.978 ;4.214] ;

[RA,RB,n,X] =liezhu(A,b),toc

运行结果为：

因为RA=RB=n，所以此方程组有唯一解.

RA =

4

RB =

4

n =

4

X =

-3.2161

1.3785

-1.0959

0.6585

Elapsed time is 0.003000 seconds.

同时，与高斯消元法和逆矩阵检验运算结果一致。

总结与体会：

通过本次实验我掌握了消元法解方程的一些基本算法以及用matlab实现矩阵的几种基本计算，对MATLAB软件有了更深的了解。

同时，在实验中，在输入向量b=[1.521 ;2.614 ;3.978 ;4.214]时错误的输成b=[1.521 ；2.614 ；3.978 ；4.214],致使程序不能运行，无法得到预期的实验结果，后经多番仔细查找，最终发现分号应为英文输入法，以后在编程时，一定更加认真仔细，不犯同样的错误！

在高斯消元法和列主消元法两种方法的选择上，需要仔细斟酌，同时我对两种方法做了如下的分析：

如果给定了一方程组如：a11*x1+a12*x2=b1;a21*x1+a22*x2=b2；如何选择算法？

通常高斯消去法有两个过程:消元过程和回代过程. 如果不选主元,消元过程第1步需要进行下面运算,a22-(a21/a11)*a12。

这是消元后所得新的第2个方程的x2的系数,如果a11较a21小的多,则a21/a11就很大,由于在计算机上编程计算或手算时,舍入误差难以避免,如a12=1/3,计算时需舍入为有限小数,比如保留8位有效数字,1/3用0.33333333代替,误差很小,但是当a21/a11很大,比如a21=10,a11=0.000001,则a21/a11=10000000,此时计算a22-(a21/a11)*a12的值,由于a12有10^(-8)的误差,则(a21/a11)*a12的误差却变为了10^(-1),误差被大大的放大了,算法不稳定,不稳定的算法计算出的结果不可靠了,所以消元时要选主元,列主元消去法是调换方程的次序,使调换后方程的a11较a21大,此时a21/a11的绝对值小于1,这样做的结果,在以后的计算中,误差不但不被放大,反而缩小,这时算法稳定,用稳定的算法计算,计算结果才可能可靠,列主元消去法是稳定算法,但通常的高斯消去法不是稳定的

算法,我认为一般用列主元消去法比较好。