解线性方程组的基本思想

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

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.mA;b;[L，U]=lu (A);X=U\(L\b)II）Cholesky分解若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即：其中R为上三角阵。

方程A*X=b 变成所以在Matlab中建立M文件如下% exp2.mA;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.mA;b;[Q，R]=qr(A);X=R\(Q\b)2．求线性齐次方程组的通解(A*X=0)在Matlab中，函数null用来求解零空间，即满足A•X=0的解空间，实际上是求出解空间的一组基（基础解系）。

在Matlab中建立M文件如下% exp4.mformat rat %指定有理式格式输出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.mclear allA;b; %输入矩阵A,b[m,n]=size(A);R=rank(A);B=[A b];Rr=rank(B);format ratif R==Rr&R==n % n为未知数的个数，判断是否有唯一解x=A\b;elseif R==Rr&R<n %判断是否有无穷解x=A\b %求特解C=null(A, r ) %求AX=0的基础解系,所得C为n-R列矩阵，这n-R列即为对%应的基础解系% 这种情形方程组通解xx=k(p)*C(:,P)(p=1…n-R)else X= No solution! % 判断是否无解end第3章线性方程组的迭代解法3.1实验目的理解线性方程组计算机解法中的迭代解法的求解过程和特点，学习科学计算的方法和简单的编程技术。

3.2概念与结论1.n阶线性方程组如果未知量的个数为n ,而且关于这些未知量x1,x2, …,x n的幂次都是一次的(线性的)那末, n 个方程a11x1+a12x2+ …+a1n x n=b1┆┆┆(1)a n1x1+a n2x2+ …+a nn x n=b n构成一个含n个未知量的线性方程组,称为n阶线性方程组。

其中,系数a11,…,a1n,a21, …,a2n, …,a n1, …,a nn和b1, …,b n都是给定的常数。

方程组(1)也常用矩阵的形式表示,写为Ax=b其中,A是由系数按次序排列构成的一个n阶矩阵,称为方程组的系数矩阵,x和b都是n维向量,b称为方程组的右端向量。

2. n阶线性方程组的解使方程组(1)中每一个方程都成立的一组数x1*,x2*, …,x n*称为式(1)的解,把它记为向量的形式,称为解向量.3. 向量范数的三种常用范数4．矩阵的四种常用范数5.谱半径 设 n ⨯n 阶矩阵A 的特征值为λ i (i=1,2,3……n),则称ρ (A) = MAX | λi |1≤i ≤ n为矩阵A 的谱半径.矩阵范数与谱半径之间的关系为: ρ (A) ≤ ||A||6.严格(行)对角占优阵A如果 矩阵 A=(a ij )满足n⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛====∑∑=≤≤∞=n n k k k nk n k kx x x x x x x x x x 21122111m ax .211222211121112111111maxmax最大特征值是行范数列范数A A A a a a a a a a a a A a A a A a A T nn n n n n n k kj n j f n k jk n j n k kj n j λλ=⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛====∑∑∑∑===≤≤∞=≤≤|a ii | > ∑ |a ij | i=1,2,……n,j=1,j ≠i则称方阵A 是严格(行)对角占优的.7.收敛定理对任意初始向量x (0)及任意右端向量 g ，由迭代x (k+1) =B x (k) +g 产生的迭代向量序列{x (k)}收敛的充要条件是谱半径ρ（B ）<18.收敛判别条件判别条件1：若||B||<1, 则迭代x (k+1) =B x (k) +g 对任何初始向量x (0)都收敛.判别条件2：如果A 为严格对角占优阵，则其 Jacobi 迭代和Seidel 迭代对任何初始向量x (0)都收敛。

判别条件3：如果A 为对称正定阵，则其 Seidel 迭代对任何初始向量x (0)都收敛。

9.迭代法的误差估计若||B||<1，则对迭代格式 x (k+1) =B x (k) +g 有3.3 程序中Mathematica 语句解释a*matrix 数a 与矩阵matrix 相乘matrix1+matrix2 矩阵matrix1和矩阵matrix2相加（注意矩阵的大小相同）matrix1.matrix2 矩阵matrix1和矩阵matrix2相乘(注意矩阵乘法的规则)Transpose[matrix] 求矩阵matrix 转置Inverse[matrix] 求矩阵(方阵) matrix 的逆)()1(*)()1()(*)(1.21.1o k k k k k x x BB x x x x B B x x --≤---≤--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 迭代法的构造过程假设a ii ≠0，依次在第i 个方程解出x i , i=1,2,⋯,n 并令c ij = -a ij /a ii (i ≠j) ， g i = b i /a ii就得到如下Jocobi 迭代格式：x 1(k+1)= c 12x 2(k)+c 13x 3(k)+⋅⋅⋅⋅ +c 1n x n (k)+g 1x 2(k+1)=c 21x 1(k) +c 23x 3(k)+⋅⋅⋅⋅ +c 2n x n (k)+g 2。

x n (k+1)=c n1x 1(k) +c n2x 2(k)+⋅⋅⋅⋅ +c n(n-1)x n-1(k) + g n若令⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=n n n n n n J g g g g x x x x c c c c c c B 212121221112000则有Jocobi 迭代的矩阵格式：x (k+1) = B J x (k) +g JB J 称为Jocobi 迭代矩阵。

Jocobi 迭代可以写成如下紧凑格式：在给定初始迭代向量x (0)后就可以进行Jocobi 迭代求解了。

2) Jacobi 迭代算法1.输入变量个数n 、初值向量x（0）、迭代精度eps 、系数矩阵A 、常数项b 和迭代最大次数nmax 2 For i=1,2,…,n2.1 如果|a ii |<eps1,则输出“迭代失败”提示并终止3. Bj ⇐ E-D -1A4. g j ⇐ D -1b5.For k=1,2,…,nmax5.1 x ⇐Bj.x0+ g j5.2 如果||x-x0||<eps ,输出解向量x ，终止；否则x（0） ⇐ x6. 如果||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.;eps=Input["输入精度控制eps="]； ni ni j j ii i k j ii ij k i a b x a a x ,,2,11)()1( =≠=+∑--=Do[If[Abs[a[[i,i]]]<eps1,t1=1;Return[],t1=0],{i,1,n}];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]<eps,Break[],x0=x],{i,1,nmax}];If[err>=eps,Print["迭代失败"]]]说明本程序用于求线性方程组Ax=b的解。