线性方程组的直接解法

利用列主元求解，调用： format long A=[420 210 140 105; 210 140 105 84; 140 105 84 70; 105 84 70 60]; b=[875 ; 539 ; 399; 319;]; [x,U]=GaussColumnMain(A,b) 得到：
则可知 x 的近似解为：
0 2 x1 x 2 x1 6 x 2 2 x3 4 x 3 x 8 x 2 3 1
2 4 5
调用： A=[2 -1 0; 1 6 -2; 4 -3 8; ]; b=[2 ; -4 ; 5;]; [x,U]=GaussOrder(A,b) 程序运行结果如下：
（n） A( 2) ....... A A = A(1)
` 在第 k-1 步结束，矩阵 A( k ) 将被构造出来；它的样子如下所示，其中用线框 起来的第 k 行和仅仅前面有线的第 k 列说明由消元过程产生的结构
(k ) 我们的任务是描述怎样从 A( k ) 得到 A( k 1) ， 为了在第 k 列的主元素 a kk 下面产生 0，
通过解方程
Ae ( 0) r ( 0)
得到。此即为迭代细化。 假如方程已用高斯消去法解出，由于舍入误差，我们不能期望所得的结果是 精确解，我们用 x ( 0) 表示它，然后计算 r ( 0) ， e ( 0 ) ，并用下面三个式子计算 x (1) ：
r ( 0 ) b Ax ( 0) ( 0) (0) Ae r x (1) x ( 0 ) e ( 0 )
运行时间 0.040390 秒
为了区别顺序高斯消元法和列主元高斯消去法，我们修改原方程组的系数为： A=[0 -1 2; 1 6 -2; 4 -3 8; ]; 调用： A=[0 -1 2; 1 6 -2; 4 -3 8; ]; b=[2 ; -4 ; 5;]; [x,U]=GaussOrder(A,b) 得到：
重复过程可以得到更精确地解 x ( 2) , x (3) ...
二、 数值实验
2.1 计算机配置： 内存：2G 操作系统：win 7 Cpu ：intel core i3，2.13GHz Matlab 版本：7.8.0.347 （R2009a ） 2.2 顺序高斯消去法： 我们利用 matlab 编制 M 文件，求解如下方程组：
我们把它下面的行减去第 k 行的倍数。 第 i=1,2， …， k 行不改变。 所以公式是：
( k 1) aij (k ) aij 若i k (k ) (k ) (k ) (k ) aij (aik / a kk )a kj 若i k 1和j k 1 0 若i k 1和j k
可见， 共迭代了 2 次就满足 1e-16 的精度要求，可见迭代细化能保证方程组的解 满足较高的精度要求，且收敛速度也较快。 最终解为：
x1 0.999999999999998 x 1.000000000000005 2 x3 1.000000000000014 x 4 0.999999999999979
x x (0) A1 (b Ax (0) ) x (0) e (0)
其 中 e (0) A1 (b Ax ( 0) ) 称 为 误 差 向 量 . 对 应 于 近 似 解 的 残 差 向 量 是
r (0) b Ax (0) ,它是可以计算的。我们不计算 A( 1) ，但是向量 e (0) A1r (0) 可以
x1 1.000000000000022 x 2 0.999999999999745 x 1.000000000000615 3
x4 0.999999999999600
以此为初始值进行迭代细化，调用： A=[420 210 140 105; 210 140 105 84; 140 105 84 70; 105 84 70 60]; x=[ 1.000000000000022; 0.999999999999745; 1.000000000000615; 0.999999999999600;]; b=[875 ; 539 ; 399; 319;]; [x1,R,E,iter,error]=IterRefine(A,b,x,1e-16,2000) 其中：x1：为最终满足精度要求的解 R：存放每一步迭代结果， E：存放每一步迭代误差 Iter：为共进行的迭代次数 Error：为迭代成功与否的标志，1 标志成功 A，b:为方程组的系数 精度要求为：1e-16 最大迭代次数为 2000 最终得到：
（1-1）
若 det(U)≠0，即 u ii ≠0（i=1,2,…n）,则（1-1）有唯一解，可以从最后一个方 程得到：
xn y n / u n
代入倒数第二个方程可得
xn1 ( y n u n1,n xn ) / u n1,n1
通常设已求得 xn , xn1 ,, xi 1 , 则方程组的第 i 个方程可得：
线性方程组的直接解法初步探讨
魏栋 20103950
一、 线性方程组的直接解法
1.1 顺序 Gauss 消去法 Gauss 消去法是求解线性方程组较为有效的方法。它主要包括两个操作：消 元和迭代。 消元即是将线性方程组转化为与其同解的上三角方程组；回代是指通 过上三角方程组逐个解出方程组的未知数。 1.1.1 消元过程： 为正式地描述高斯算法过程，我们把它理解为一个连续 n-1 步主步产生的如 下矩阵系列：
说明此时顺序高斯消去法不适用了 利用列主元法，调用： tic; A=[0 -1 2; 1 6 -2; 4 -3 8; ]; b=[2 ; -4 ;
5;]; [x,U]=GaussColumnMain(A,b) toc; 得到：
说明此时列主元高斯消去法仍适用 2.4 迭代细化： 我们编制 matlab 程序，选取方程组为： A=[420 210 140 105; 210 140 105 84; 140 105 84 70; 105 84 70 60]; b=[875 ; 539 ; 399; 319;];
三、 附录
3.1 顺序高斯消去法的程序： GaussOrder.m 3.2 列主元高斯消去法的程序： GaussColumnMain.m
3.3 迭代细化的程序： IterRefine.m
(k ) 当消元过程进行到第 k 步时，按顺序高斯消去法，应该选取 a kk 为主元，使
得第 k 个方程不变。但注意到，此时，若该主元为 0 或是其绝对值很小时就会出 现上述问题，所以我们选择第 k 列的元素 ai(kk ) （i=k,k+1,…,n）中绝对值最大的 元素为主元，即令：
(k ) (k ) a rk max aik k i n
xi yi
j i 1
u
n
ij
xj / u ii （i=n-1,n-2,…1）
1.2 选主元： 由于顺序高斯消去法存在一个明显的缺点，当有一步的消元过程中，主元是 0 时，消元过程就无法进行下去了，并且当主元很小时，由于舍入误差的影响， 也可能引起很大的误差，从而使高斯消去法失效。 解决问题问题的方法之一是重新选取主元，可以有行尺度高斯消元法，列尺 度高斯消元法，或是全主元高斯消去法。这里我们以列主元为例进行介绍：
然后，我们取 U= A( n ) 并定义 L 为：
(k ) (k ) aik / a kk 若i k 1 lik 1 若i k 0 若i k 1
这里 A=LU 是矩阵 A 的标准高斯分解，L 为单位下三角阵，而 U 为上三角阵。 注意，当消元过程的任何一个主元素为 0 时将中断。
2 4 5
调用： tic; A=[2 -1 0; 1 6 -2; 4 -3 8; ]; b=[2 ; -4 ; 5;]; [x,U]=GaussColumnMain(A,b) toc; 得到：
可得最终解仍为：
x1 0.6200 x 2 0.7600 x 0.0300 3
如果 a r(kk ) =0，说明矩阵不可逆，即方程组的解不确定；否则，当 r>k 时，则其在
(k ) 增广矩阵中交换第 k 行和第 r 行，使 a r(kk ) 成为主元（此时变成了 a k ，然后再按 k ）
Gauss 消去法进行消元运算，此即为列主元高斯消去法。 1.3 迭.1.2 回代过程： 设原方程为： AX=Y 假设按上述步骤消元后得到 U，此时解方程组即为解：UX=Y 式中： u1n u11 u12 u1,n 1 y1 y u 22 u 2,n 1 u 2n 2 U ，Y y u n 1,n 1 u n 1,n n 1 u nn yn
可得方程组的解为：
x1 0.6200 x 2 0.7600 x 0.0300 3
运行时间为：0.015229 秒
2.3 列主元高斯消去法 我们仍利用 matlab 求解方程组：
0 2 x1 x 2 x1 6 x 2 2 x3 4 x 3 x 8 x 2 3 1