解方程组基本思想.docx

解方程组基本思想.docx

《解方程组基本思想.docx》由会员分享，可在线阅读，更多相关《解方程组基本思想.docx（30页珍藏版）》请在冰豆网上搜索。

解方程组基本思想

四：

基本方法

基本思路将在解题的过程中得到体现。

1．（求线性方程组的唯一解或特解），这类问题的求法分为两类：

一类主要用于解低阶稠

密矩阵——直接法；一类是解大型稀疏矩阵——迭代法。

1.1利用矩阵除法求线性方程组的特解（或一个解）

方程：

AX=b，解法：

X=A\b，（注意此处’\’不是’/’）

例1-1   求方程组的解。

解:

   A=;=；b=（1,0,0,0,1）’

由于>>rank（A）=5,rank（）=5   %求秩，此为R（A）=R（）>=n的情形，有唯一解。

>>X=A\b     %求解X=（2.2662,-1.7218,1.0571,-0.5940,0.3188）’或用函数rref

求解，>>sv=rref（A:

b）;所得sv的最后一列即为所要求的解。

1．2利用矩阵的LU、QR和cholesky分解求方程组的解

这三种分解，在求解大型方程组时很有用。

其优点是运算速度快、可以节省磁盘空间、节省内存。

I）LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行交换

）和上三角矩阵的乘积。

即A=LU，L为下三角阵，U为上三角阵。

则：

A*X=b         变成L*U*X=b

所以X=U\（L\b）    这样可以大大提高运算速度。

命令   [L，U]=lu（A）

在matlab中可以编如下通用m文件：

在Matlab中建立M文件如下

%exp1.m

A;b;

[L，U]=lu（A）;

X=U\（L\b）

II）Cholesky分解

若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即：

    其中R为上三角阵。

方程   A*X=b   变成    所以   

在Matlab中建立M文件如下

%exp2.m

A;b;

[R’，R]=chol（A）;

X=R\（R’\b）

III）QR分解

对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形

式，即：

A=QR

方程   A*X=b   变形成    QRX=b

所以    X=R\（Q\b）

上例中   [Q,R]=qr（A）

X=R\（Q\B）

在Matlab中建立M文件如下

%exp3.m

A;b;

[Q，R]=qr（A）;

X=R\（Q\b）

2．求线性齐次方程组的通解（A*X=0）

在Matlab中，函数null用来求解零空间，即满足A•X=0的解空间，实际上是求出解空

间的一组基（基础解系）。

在Matlab中建立M文件如下

%exp4.m

formatrat     %指定有理式格式输出

A;b=0;

r=rank（A）;

bs=null（A，‘r’）;%一组基含（n-r）个列向量

%k,k,……,k

%X=k*bs（:

1）+k*bs（:

2）+……+k*bs（:

n-r）方程组的通解

pretty（X）      %让通解表达式更加精美

3   求非齐次线性方程组的通解（A*X=b）

非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。

因此，步骤为：

第一步：

判断AX=b是否有解，（利用基本思路的第一条）

若有解则进行第二步

第二步：

求AX=b的一个特解

第三步：

求AX=0的通解

第四步：

AX=b的通解为：

AX=0的通解加上AX=b的一个特解。

在Matlab中建立M文件如下

%exp4.m

clearall

A;b;                   %输入矩阵A,b

[m,n]=size（A）;

R=rank（A）;

B=[Ab];

Rr=rank（B）;

formatrat

ifR==Rr&R==n           %n为未知数的个数，判断是否有唯一解

x=A\b;

elseifR==Rr&R

x=A\b                  %求特解

C=null（A,r）          %求AX=0的基础解系,所得C为n-R列矩阵，这n-R列即为对

%应的基础解系

                       %这种情形方程组通解xx=k（p）*C（:

P）（p=1…n-R）

elseX=Nosolution!

   %判断是否无解

end

第3章线性方程组的迭代解法

3.1实验目的

理解线性方程组计算机解法中的迭代解法的求解过程和特点，学习科学计算的方法和简单的编程技术。

3.2 概念与结论

1.    n阶线性方程组

如果未知量的个数为n,而且关于这些未知量x1,x2,…,xn的幂次都是一次的（线性的）那末,n个方程

a11x1+a12x2+…+a1nxn=b1

┆┆┆

（1）

an1x1+an2x2+…+annxn=bn

构成一个含n个未知量的线性方程组,称为n阶线性方程组。

其中,系数a11,…,a1n,a21,…,a2n,…,an1,…,ann和b1,…,bn都是给定的常数。

方程组

（1）也常用矩阵的形式表示,写为

Ax=b

其中,A是由系数按次序排列构成的一个n阶矩阵,称为方程组的系数矩阵,x和b都是n维向量,b称为方程组的右端向量。

2.n阶线性方程组的解

使方程组

（1）中每一个方程都成立的一组数x1*,x2*,…,xn*称为式

（1）的解,把它记为向量的形式,称为解向量.

3.向量X数的三种常用X数

4．矩阵的四种常用X数

5.谱半径

设nn阶矩阵A的特征值为i（i=1,2,3……n）,则称

（A）=MAX|i|

1in

为矩阵A的谱半径.

矩阵X数与谱半径之间的关系为:

（A）||A||

6.严格（行）对角占优阵A

如果矩阵A=（aij）满足

n

|aii|>|aij|i=1,2,……n,

j=1,ji

则称方阵A是严格（行）对角占优的.

7.收敛定理

对任意初始向量x（0）及任意右端向量g，由迭代x（k+1）=Bx（k）+g产生的迭代向量序列{x（k）}收敛的充要条件是谱半径（B）<1

8.收敛判别条件

判别条件1：

若||B||<1,则迭代x（k+1）=Bx（k）+g对任何初始向量x（0）都收敛.

判别条件2：

如果A为严格对角占优阵，则其Jacobi迭代和Seidel迭代对任何初始向量x（0）都收敛。

判别条件3：

如果A为对称正定阵，则其Seidel迭代对任何初始向量x（0）都收敛。

9.迭代法的误差估计

若||B||<1，则对迭代格式x（k+1）=Bx（k）+g　有　　

3.3程序中Mathematica语句解释

a*matrix数a与矩阵matrix相乘

matrix1+matrix2矩阵matrix1和矩阵matrix2相加（注意矩阵的大小相同）

matrix1.matrix2矩阵matrix1和矩阵matrix2相乘（注意矩阵乘法的规则）

Transpose[matrix]求矩阵matrix转置

Inverse[matrix]求矩阵（方阵）matrix的逆

DiagonalMatrix[list]使用列表list中的元素生成一个对角矩阵.

IdentityMatrix[n]生成n阶单位矩阵

Max[x]求向量x中元素的最大值

3.4方法、程序、实验

解线性方程组的迭代法是将线性方程组Ax=b化为等价线性方程组

x=Bx+f

再由矩阵迭代格式

x（k+1）=Bx（k）+f

构造向量序列{x（k）}来求线性方程组解的。

如果得出的向量序列{x（k）}收敛至某个向量x*，则可得该向量x*就是所求方程组Ax=b的准确解.

线性方程组的迭代法主要有Jocobi迭代法、Seidel迭代法和超松弛（Sor）迭代法。

1.Jocobi迭代法

1）Jocobi迭代法的构造过程

假设aii0，依次在第i个方程解出xi,i=1,2,,n并令

cij=-aij/aii（ij），gi=bi/aii

就得到如下Jocobi迭代格式：

x1（k+1）=c12x2（k）+c13x3（k）++c1nxn（k）+g1

x2（k+1）=c21x1（k）+c23x3（k）++c2nxn（k）+g2

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

xn（k+1）=1x1（k）+2x2（k）++（n-1）xn-1（k）+gn

若令

则有Jocobi迭代的矩阵格式：

x（k+1）=BJx（k）+gJ

BJ称为Jocobi迭代矩阵。

Jocobi迭代可以写成如下紧凑格式：

在给定初始迭代向量x（0）后就可以进行Jocobi迭代求解了。

2）Jacobi迭代算法

1.输入变量个数n、初值向量x（0）、迭代精度eps、系数矩阵A、常数项b和迭代最大次数nmax

2Fori=1,2,…,n

2.1如果|aii|

3.BjE-D-1A

4.gjD-1b

5.Fork=1,2,…,nmax

5.1xBj.x0+gj

5.2如果||x-x0||

6.如果||x-x0||>eps,输出迭代失败，终止。

3）Jacobi迭代法程序

Clear[a,b,x];

nmax=500;

n=Input[“线性方程组阶数n=”];

a=Input["系数矩阵A="]；

b=Input["常数项b="]；

x0=Input["输入迭代初值向量x0"]；

eps1=0.000001;

eps=Input["输入精度控制eps="]；

Do[If[Abs[a[[i,i]]]

If[t1==1,

Print["Jacobi迭代法失效"],

d=DiagonalMatrix[Table[a[[i,i]],{i,1,n}]];

d1=Inverse[d];

bj=IdentityMatrix[n]-d1.a;

gj=d1.b;

Do[x=bj.x0+gj;

err=Max[Abs[x-x0]];

Print["x=",x//N,"i=",i,"err=",err//N];

If[N[err]

{i,1,nmax}];

If[err>=eps,Print["迭代失败"]]

]

说明本程序用于求线性方程组Ax=b的解。

程序执行后，先通过键盘输入线性方程组阶数n、系数矩阵A、常数项b、迭代初值向量x0和输入精度控制eps，程序即可给出每次迭代的次数和对应的迭代向量序列x（k）,其中最后输出的结果即为所求的根。

如果迭代超出500次还没有求出满足精度的根则输出迭代失败提示，如果出现主对角线元素aii=0给出Jacobi迭代法失效提示。

程序中变量说明

x0：

存放初始向量和迭代过程中的向量x（k）

x:

存放迭代过程中的向量x（k+1）

nmax:

存放迭代允许的最大次数

err:

存放误差||x-x0||

t1:

临时变量

注：

迭代最大次数可以修改为其他数字。

4）例题与实验

例1.用Jacobi迭代法解如下线性方程组

5x1+2x2+x3=-12

-x1+4x2+2x3=20

2x1-3x2+10x3=3

要求误差||x（k+1）-x（k）||<10-4，并用取不同初值的方法实验观察迭代收敛的情况。

解：

执行Jacobi迭代法程序后在输入的四个窗口中按提示分别输入:

3、{{5,2,1},{-1,4,2},{2,-3,10}}、{-12,20,3}、{0,0,0}、0.0001

每次输入后用鼠标点击窗口的“OK”按扭，得如下输出结果：

x={-2.4,5.,0.3}i=1err=5.

x={-4.46,4.25,2.28}i=2err=2.06

x={-4.556,2.745,2.467}i=3err=1.505

x={-3.9914,2.6275,2.0347}i=4err=0.5646

x={-3.85794,2.9848,1.88653}i=5err=0.3573

x={-3.97123,3.09225,1.96703}i=6err=0.113286

x={-4.03031,3.02368,2.02192}i=7err=0.0685705

x={-4.01386,2.98146,2.01316}i=8err=0.042216

x={-3.99522,2.98995,1.99721}i=9err=0.0186374

x={-3.99542,3.00259,1.99603}i=10err=0.0126367

x={-4.00024,3.00313,1.99986}i=11err=0.0048186

x={-4.00122,3.00001,2.00099}i=12err=0.00312067

x={-4.0002,2.9992,2.00025}i=13err=0.00102319

x={-3.99973,2.99983,1.9998}i=14err=0.000625697

x={-3.99989,3.00017,1.99989}i=15err=0.00034136

x={-4.00005,3.00008,2.00003}i=16err=0.000155236

x={-4.00004,2.99997,2.00003}i=17err=0.000106099

x={-4.,2.99997,2.}i=18err=0.0000414468

此结果说明迭代18次，求得误差为err=0.0000414468的近似解，最后显示的近似解向量为x={-4.,2.99997,2.}，它表示所求解为

x1=-4，x2=2.99997，x3=2。

本题的准确解为

x1=-4，x2=3，x3=2。

如果将如上输入的初值改为{21,-18,30}，执行Jacobi迭代法程序后得如下输出结果：

x={-1.2,-4.75,-9.3}i=1err=39.3

x={1.36,9.35,-0.885}i=2err=14.1

x={-5.963,5.7825,2.833}i=3err=7.323

x={-5.2796,2.09275,3.22735}i=4err=3.68975

x={-3.88257,2.06642,1.98375}i=5err=1.39703

x={-3.62332,3.03749,1.69644}i=6err=0.97106

x={-3.95428,3.24595,1.93591}i=7err=0.330963

x={-4.08556,3.04347,2.06464}i=8err=0.202475

x={-4.03032,2.94629,2.03015}i=9err=0.0971858

x={-3.98455,2.97734,1.98995}i=10err=0.0457716

x={-3.98893,3.00889,1.99011}i=11err=0.0315451

x={-4.00158,3.00771,2.00045}i=12err=0.0126504

x={-4.00318,2.99938,2.00263}i=13err=0.00833246

x={-4.00028,2.99789,2.00045}i=14err=0.00289753

x={-3.99925,2.99971,1.99942}i=15err=0.0018145

x={-3.99977,3.00048,1.99976}i=16err=0.000770763

x={-4.00014,3.00018,2.0001}i=17err=0.000375926

x={-4.00009,2.99992,2.00008}i=18err=0.000261658

x={-3.99998,2.99994,1.99999}i=19err=0.000107579

x={-3.99997,3.00001,1.99998}i=20err=0.0000714045

从计算结果可以看到虽然所取的初值不同，但迭代总是收敛相同的结果。

考察本题系数矩阵可以发现它是严格对角占优矩阵，由收敛判别法，Jacobi迭代对任意初值都是收敛的，我们通过实际计算验证了这个结果。

例2.通过Jacobi迭代序列观察用Jacobi迭代法解如下线性方程组

x1+2x2+2x3=-12

-x1+4x2+x3=20

2x1-3x2+x3=3

的收敛性。

解：

执行Jacobi迭代法程序后在输入的四个窗口中按提示分别输入:

3、{{1,2,2},{-1,4,1},{2,-3,1}}、{-12,20,3}、{0,0,0}、0.0001

每次输入后用鼠标点击窗口的“OK”按扭，得如下输出结果：

x={-12.,5.,3.}i=1err=12.

x={-28.,1.25,42.}i=2err=39.

x={-98.5,-12.5,62.75}i=3err=70.5

x={-112.5,-35.3125,162.5}i=4err=99.75

…………………………………………….

{246.719,-120.176,696.719}i=8err=763.5

x={-1165.09,-107.5,-850.965}i=9err=1547.68

x={1904.93,-73.5303,2010.67}i=10err=3070.02

x={-3886.28,-21.4355,-4027.45}i=11err=6038.12

x={8085.77,40.2917,7711.26}i=12err=11972.1

x={-15515.1,98.6279,-16047.7}i=13err=23758.9

x={31886.1,138.141,31329.1}i=14err=47401.2

x={-62946.5,144.247,-63354.7}i=15err=94832.5

x={126409.,107.068,126329.}i=16err=189683.

x={-252883.,25.0775,-252494.}i=17err=379292.

x={504925.,-92.4305,505845.}i=18err=758339.

…………………………………….

通过观察迭代过程中的误差是不断变大的特点可以知道本题的Jacobi迭代序列是不收敛的，因此，本题线性方程组不能用Jacobi迭代法求解。

这个实验说明并不是每个线性方程组都能用Jacobi迭代法求解。

2.    Seidel迭代

1）Seidel迭代的构造过程

为了加快收敛速度，同时节省计算机的内存，对Jocobi迭代作如下的改进：

每算出一个分量的近似值，立即用到下一个分量的计算中去，即用迭代格式：

x1（k+1）=c12x2（k）+c13x3（k）++c1nxn（k）+g1

x2（k+1）=c21x1（k+1）+c23x3（k）++c2nxn（k）+g2

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

。

xn（k+1）=1x1（k+1）+2x2（k+1）++（n-1）xn-1（k+1）+gn

如上迭代可以写成如下紧凑格式：

这样所得的迭代法就称为Seidel迭代法。

利用Ax=b及A=L+D+U,其中D为对角矩阵,L,U分别为严格下，上三角矩阵．则有Seidel迭代法的矩阵形式为

x（k+1）=Bsx（k）+gs其中：

Bs=-（L+D）-1U，gs=D-1b

Seidel迭代矩阵为Bs=-（L+D）-1U。

在给定初始迭代向量x（0）后就可以进行Seidel迭代求解了。

Jacobi迭代和Seidel迭代格式可表述为统一形式:

2）Seidel迭代算法

1.输入变量个数n、初值向量x（0）、迭代精度eps、系数矩阵A、常数项b和迭代最大次数nmax

2.Fori=1,2,…,n

2.1如果|aii|

3.Fori=1,2,…,nmax

3.1Fori=1,2,…,n

3.2如果||x（k+1）-x（k）||

4.如果||x-x0||

3）Seidel迭代法程序

Clear[a,b,x];

nmax=500;

n=Input[“线性方程组阶数n=”];

a=Input["系数矩阵A="]；

b=Input["常数项b="]；

x0=Input["输入迭代初值向量x0"]；

eps1=0.000001;

eps=Input["输入精度控制eps="]；

x=x0;

Do[If[Abs[a[[i,i]]]

If[t1==1,

Print["Seidel迭代法失效"],

Do[Do[u1=Sum[a[[i,j]]*x[[j]],{j,1,i-1}];

u2=Sum[a[[i,j]]*x0[[j]],{j,i+1,n}];

x[[i]]=（b[[i]]-u1-u2）/a[[i,i]],

{i,1,n}];

err=Max[Abs[x-x0]];

Print["x=",x//N,"k=",k,"err=",err//N];

If[N[err]

{k,1,50}];

If[err>=eps,Print["迭代失败"]]

]

说明本程序用于求线性方程组Ax=b的解。

程序执行后，先通过键盘输入线性方程组阶数n、系数矩阵A、常数项b、迭代初值向量x0和输入精度控制eps，程序即可给出每次迭代的次数和对应的迭代向量序列x（K）,其中最后输出的结果即为所求的根。

如果迭代超出500次还没有求出满足精度的根则输出迭代失败提示，如果出现主对角线元素aii=0给出Jacobi迭代法失效提示。

程序中变量说明

x0：

存放初始向量和迭