王向东数学实验课本(可编辑)2-13

实验十三 线性方程组

【实验目的】

1． 了解线性方程组的基本概念。

2． 了解线性方程组解的结构。

3． 学习掌握MATLAB 软件有关的命令。

【实验内容】

解方程组⎩⎨⎧=-=+2

34232y x y x 【实验准备】

1．线性方程组的基本概念

含有n 个未知数n x x x ,,,21 的线性方程组

⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++n

n nn n n n n n n b x a ax x a b x a x a x a b x a x a x a

2112222212111212111 可以写成如下矩阵的形式

b AX =，

式中

⎥⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n nn n n n n b b b b a a a a a a a a a A 21212222111211, 当0=b 时称方程组为齐次的，0≠b 时则称方程组为非齐次的。

2．线性方程组解的结构

线性方程组的解可能出现三种情形：无解、有唯一解和有无穷多组解。这主要取决于系数矩阵A 的秩)(A R 与增广矩阵)(b A 的秩)(b A R 是否相等、秩与变两个数是否相等，具体地： 若)()(b A R A R ≠，则线性方程组无解； 若n b A R A R ==)()(（n 为变量个数），则线性方程组有唯一解； 若n b A R A R <=)()(，则线性方程组有无穷多组解。

3．线性方程组的解法

线性方程组的解法大概可以划分为两类：直接消去法（这里主要指高斯消去法）和迭代

数值解法。

直接消去法：理论上，经过有限次算术运算能够求出方程组的精确解，实际上由于计算存在舍入误差，因此可能得到的只是近似解。高斯消去法可以用LU 分解来表示。如果方程组的系数矩阵A 满足一定的条件，则A 可以分解为一个下三角矩阵L 与一个上三角矩阵U 的乘积LU A =，于是由b AX =得到

.)(1111b L U b LU b A X ----===

由于上三角矩阵和下三角矩阵的逆很容易求得，因此，这是一种很好的思想。

迭代法：常用的迭代法有雅可比迭代法和高斯─赛德尔迭代法等。将b AX =等价变形为

f BX X +=

由此构造出迭代公式

f X B X n n +=+)()1(

可以证明，在B 满足谱半径1)(

雅可比迭代：设U D L A ++=（D 对角矩阵,L 为上三角矩阵,U 为上三角矩阵）,由b AX =得到

,])([,)(1f BX b X U L D X b X U D L +=++-==++-

其中11),(--=+-=D f U L D B .

赛德尔迭代:将上述的迭代公式写成另外一种形式

f BX b UX L D X b UX X D L +=+-+=+-=+-)()(,)(1

其中b L D F U L D B 11)(,)(--+=+-=.

4．线性方程组求解的MATLAB 命令

MATLAB 中主要用inv,null 分别求矩阵的逆和齐次方程组的基础解系。

可以用help inv,help null 查阅有关这些命令的详细信息

【实验方法与步骤】

练习1 解方程组⎩⎨

⎧=-=+2

34232y x y x . 相应的MA TLAB 代码为： >>clear;

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

>>X=A\b %左除法,解方程组AX=b

算得(x,y)=(6.6364,8.1818).如果用inv 命令, 相应的MATLAB 代码为：

>>clear;

>>A=[1 2;4 -3]; b=[23;2]; x=inv(A)*b

仍算得(x,y)=(6.6364,8.1818). 练习２求不定方程组⎩⎨

⎧=+-=++2

2322z y x z y x 的一个特解.相应的MATLAB 代码为： >>clear; >>A=[1 2 1;3 -2 1]; b=[2;2]; x=A\b

求得一个特解)0,5.0,1(),,(=z y x .

练习3求超定方程组⎪⎩

⎪⎨⎧=-=-=+242312y x y x y x 的最小二乘解.

显然,此方程组无解.所谓方程组b AX =的最小二乘解指使得向量b AX -的长度达到最小的解.相应的MA TLAB 代码为：

>>clear;

>>A=[1 2;3 -2;1 -1]; b=[1;4;2]; x=A\b

求得一最小二乘解)1757.0,2838

.1(),(-=y x . 练习4求奇异方程组⎩

⎨⎧-=--=+24212y x y x 的一个特解. 如果用前面的方法,相应的MATLAB 代码为：

>>clear;

>>A=[1 2;-2 -4]; b=[1;-2]; x=A\b

结果为

x=Inf %说明不能直接求解

Inf

如果用

>>clear;

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

>>x=A\b

仍能求出一特解(x,y)=(0,0.5000).

练习5 求方程组⎪⎩⎪⎨⎧-=+--=-++-=-+-122114321

43214321x x x x x x x x x x x x 的通解.

我们首先看方程组解的结构,计算系数矩阵和增广矩阵的秩,相应的MATLAB 代码为： >>clear;

>>a=[1 -1 1 -1;-1 1 1 -1;2 -2 -1 1]; b=[1;1;-1];

>>rank(a) %系数矩阵的秩

>>rank([a,b]) %增广矩阵的秩

计算表明,系数矩阵和增广矩阵的秩都为2,小于变量的个数4,说明原方程组有无穷组解.有几种方法求原方程组的通解.一种是用rref 命令化为行最简形式求解. 相应的MA TLAB 代码为：

>>clear;

>>a=[1 -1 1 -1;-1 1 1 -1;2 -2 -1 1]; b=[1;1;-1];

>>rref([a,b])

得最简形式为

ans =

1 -1 0 0 0

0 0 1 -1 1

0 0 0 0 0

从而通解为,1,4321+==x x x x 42,x x 为自由变量.

另外一种方法用null 命令求齐次方程组的一个基础解系.由于非齐次方程的通解等于齐次方程的通解加非齐次方程的一个特解,故可用如下的MA TLAB 代码:

>>clear;

>>a=[1 -1 1 -1;-1 1 1 -1;2 -2 -1 1]; b=[1;1;-1];

>>x0=a\b %齐次方程的一个特解

>>x1=null(a) %非齐次方程的通解

结果为

x0 =

1

x1 =

-0.7071 0

-0.7071 -0.0000

-0.0000 0.7071

0 0.7071

原方程组的通解为

)7071.0,7071.0,0,0()0,0,7071.0,7071.0()0,1,0,0(),,,(214321c c x x x x +--+= 式中21,c c 为任意常数.