matlab追赶法解三对角方程组_数值计算方法 第三章 线性代数方程组的直接解法（1）...

写在章前：求解线性方程组的数值方法大体上可分为直接法和迭代法两大类.。

直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解法，又称为精确法；迭代法则是采取逐次逼近的方法，亦即从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列，其极限才是方程组的精确解，只经过有限次运算往往得不到精确解。

在这一章，我们将主要介绍解线性方程组的一些基本的直接法。

一、Gauss消去法

1、三角形方程组的解法

①下三角方程组解法

考虑

，这里

已知，

未知，且

非奇异。则我们可由上往下解方程且不断代入下面的方程。

这种方法叫前代法。

for 

运算加减乘除次数为：

②上三角方程组解法

同样的考虑上三角方程组

,我们采用

回带法

for 

对于一般线性方程组

如果我们能够将 A 分解为：

，即一个下三角阵

与一个上三角阵

的乘积，那么原方程组的解便可由下面两步得到：

(1) 用前代法解

得

(2) 用回代法解

得

所以对于求解一般的线性方程组来说，关键是如何将 A 分解为一个下三角阵 L 与一个

上三角阵 U 的乘积，这正是我们本节的目标。

2、高斯变换

欲把一个给定的矩阵 A 分解为一个下三角阵 L 与一个上三角阵 U 的乘积，最自然的做法便是通过一系列的初等变化，逐步将 A 约化为一个上三角阵，而又能保证这些变换的乘积是一个下三角矩阵。

①算法思路：

我们先回顾一下初等矩阵：

图源水印

矩阵左乘表示行变换，右乘表示列变换。

我们这样分步考虑，首先找到一个下三角矩阵能把第一个列向量化到正确的形式。

这很简单。从左乘表示行变换来看，我们横着读矩阵

。从保证变换矩阵为下三角矩阵的角度来看，为了使变换矩阵为下三角矩阵，我们需要使变换矩阵的每行的后半部分的元素为0，也就是说，第

行变换后的线性组合不能包含

行之后的行向量。

从变换结果来看，我们需要变换后除第一个分量外都为0。

综合起来就是，用前

行的线性组合消掉第

行。

很自然也很简单的方法是，只用第一个元素消掉后面的所有元素。（按理说应该有很多种组合方法。但是为了复合矩阵容易求，我们希望变换矩阵非零元素最好集中在一列，这样在求多个下三角矩阵的乘积时能保证不做非零重叠位置的加法。）

用

消掉后面的n-1个元素。这样

的第一行和第一列就化完了。

我们只需看剩下的n-1阶的部分，同上可用

消掉第二列后面的n-2个元素。

以此类推。写出过程中的初等变换矩阵：

我们要求的是

,于是需要把对A的变换乘逆搞到右边去。

初等矩阵的逆和复合十分好求：

得到

分解为：

②算法的实现：

首先考虑存储。由于每次消元都要先计算消元的系数，这时需要给这个数分配内存。考虑到消元法矩阵的下矩阵都将被消为0，无用，我们可以利用这些位置存储系数。

举个栗子：

分解

算法代码：

for 

③适用情况

通常称 Gauss 消去过程中，当前的

为主元素，也就是。显然，当且仅当主元均不为零时，算法才能进行到底。那么自然要问：给定的矩阵 A 满足什么条件，才能保证所有主元素均不为零？这一问题可由下面的定理回答。

定理1：主元素均不为零的充分必要条件是 A 的前 k 阶顺序主子阵都是非奇异的。

证明：k=1时显然成立；

设用归纳法，设k-1时成立，则k阶时，k阶顺序主子式具有这样的形式

所以

即

不为0的充要条件是

不为0；由归纳法得证。

（主要证明思路就是进行的这种初等行变换不改变行列式的值，且是下三角，对角线元素均不为0和行列式不为0是充要条件）

将这个定理与前面的讨论相结合，就得到了下面的矩阵三角分解存在的充分条件.。

定理2：若

的顺序主子阵

均非奇异，则存在唯一的单位下三角阵

和上三角阵

使得

④算法比较：

cramer：

n+1个矩阵的行列式与n次除法。每个行列式有n！个元素，每个元素参与n-1次乘法

Gauss：

次乘除法，

次加减法

3、其他的三角分解

简单介绍：

其实三角分解分为两种。如果存在单位上三角阵L和下三角阵 U使得

，则称之为

的 Doolittle 分解；如果存在上三角阵 L 和单位下三角阵 U 使得 ，则称之为的

的 Crout 分解。以上我们只介绍Doolittle分解。

在大型稀疏方程组中，最常见的是带状方程组，其系数矩阵是带状矩阵，非零元素仅集中在对角线附近的带状区域内。

矩阵A是带状矩阵，上半带宽为3，下半带宽为2，总带宽为6

特别的，当上下带宽相等时我们A称方程组为等带宽方程组。总带宽为3的等带宽方程组我们又叫三对角方程组。

算法优化思想：

1、矩阵的存储

带状矩阵的压缩存储：带方程组的系数矩阵是压缩矩阵，应压缩存储。可把对角线转为行或列存储。

①对角线转为行，列标不变的压缩存储：把红色部分向中间靠，直到对角元对齐

压缩前后矩阵的下标的对应关系：

②对角线转为列，行标不变的压缩存储：把蓝色部分向右压齐

压缩前后矩阵的下标的对应关系：

2、以三对角矩阵为例介绍一次算法具体改进（追赶法）：

宏观理解：如果对带状矩阵进行三角分解，所得两个三角矩阵保持带状结构，我们就可以对算法进行优化。

保带状结构定理 ： 设

为上半带宽为

，下半带宽为

的带状矩阵。其

阶顺序主子式均不为零。则

有唯一的三角分解式

，其中

为下半带宽为

的单位下三角阵，

为上半带宽为

的上三角阵。

根据此定理，三角分解带宽以外的元素都为0，不必再参与计算。

在数值计算中，如三次样条插值或用差分方法解常微分方程边值问题，常常会遇到求解三对角方程组。我们改进算法得到追赶法。

追赶法的基本思想与Gauss消去法及三角分解法相同,只是由于系数中出现了大量的零，计算中可将它们撇开,从而使得计算公式简化，也大大地减少计算量。

数学思维：一图读懂三对角分解(Doolittle)：

由于A只需消去ai变为上三角,所以只需用对角元消掉对角元下面的一个元素，L与U都保持带状结构

当

时解得

原方程组转化为：

(emmm前为追后为赶，so )

算法实现：存储时位置互不冲突，在A原位置存储即可。

代码就是还没写！在网上先荡了一个hhh

%追赶法解三对角线性方程组

二、选主元的三角分解

大家知道，对于方程组来说，只要非奇异，方程组就存在唯一的解.。然而，非奇异并不能保证其顺序主子阵均非奇异.。因此,非奇异并不能保证 Gauss 消去过程能够进行到底。这样，我们的问题自然便是，怎样修改算法才能使其适应于非奇异矩阵呢？此外，在算法中计算时，位于分母上的主元虽不为零但很小时，由一章绪论我们知道，较小数作为分母，舍入误差会增大。如何解决这些问题呢？

解决方法是：如果出现小的主元，我们就选择一个合适的元素，交换 A 的行和列，将此元素换到主元位置上。选主元的方法有：列主元消去法，全主元消去法等。

①列主元消去法：

（

，

,

）

宏观理解：即首先

进行行交换，把所在列绝对值最大的元素换到主元位置上，相应的

应进行行交换。这时的矩阵满足主元素的要求，再进行

分解。整个过程如下：

• 选主元：在将选主元的一列中，找出绝对值最大的元素，设这个元素为主元
• 交换行；
• LU分解解方程；

数学思维：整个过程中，A的变化为

其中，L为单位上三角矩阵，P为行置换矩阵，U是上三角矩阵

令

,则

此为矩阵A的列三角分解。可用归纳法证明L是一个单位下三角矩阵，是Doolittle分解。

算法实现：

attention：记录行交换的一种方法——用向量

记录，第一行与第三行交换，则记为

,第二行不换，则记为

;也可以直接在原矩阵上做修改。

function

②全主元消去法：

宏观理解：即首先

进行行和列交换，每次选取余下方阵中绝对值最大的元素。相应的

应进行行交换。这时的矩阵满足主元素的要求，再进行

分解。整个过程如下：

• 选主元：在未选主元的方阵中，找出绝对值最大的元素，设这个元素为主元
• 交换行与列；
• LU分解解方程；

数学思维：整个过程中，A的变化为

其中，L为单位上三角矩阵，P为行置换矩阵，U是上三角矩阵

(因为行列置换只在k行k列之后进行，所以这里列置换矩阵Q只是把单位矩阵k行和j(j>k)行交换了一下，是单位下三角矩阵。

令

,则

。令

，则

此为矩阵A的全三角分解。

定理：若

为方阵，则存在n阶排列矩阵P和Q，单位下三角阵L和上三角阵U，则存在唯一的单位下三角阵

和上三角阵

使得

，而且L的所有元素的绝对值不超过1，U的非零对角元的个数等于矩阵A的秩。

算法实现：

function

三、对称正定矩阵Cholesky分解

平方根法与改进的平方根法不仅计算量仅是Gauss消去法的一半，其数值稳定性良好，是求解中小型稠密对称正定线性方程组的好办法。

①一般的平方根法

这种算法的优势在于，当需要针对同一个A解相对不同b的不同x的时候，新增的部分都可以以

复杂度完成。

cholesky分解定理：A为是对称正定矩阵，存在非奇异下三角矩阵L，使

，其中

的对角线元素均大于0。

称为矩阵A的cholesky分解或平方根分解，矩阵L称为A的cholesky因子或平方根因子。

计算步骤：解方程整个过程如下

• 求A的cholesky分解
  ，方程组变为
• 计算下三角形方程组
• 求解上三角形方程组

分解的数学推导：

把

写成矩阵形式，比较两边对应元素

这个过程其实也说明了Cholesky分解的唯一性。

②改进的平方根法

平方根法用到开方运算，为避免开方，我们可以进行如下的分解：

，其中L是单位下三角矩阵，D是对角元素均为正数的对角矩阵。这一分解称作

分解，是cholesky分解的变形

计算步骤：

• 求出
  分解，方程组变为
• 求解下三角方程组
• 求解对角方程组
• 求解上三角方程组

分解的数学推导：

在介绍LDL分解之前首先介绍LDU分解：

在LU的基础上， 如果我们再进一步，引入对角矩阵D， 那么LU分解就变成了LDU分解。

这里仔细想想对角矩阵的作用，正着想，单位矩阵首先左乘了U，相当于对U的行做了倍数不同的乘法。所以逆过来想，用对角矩阵提出行的特定倍数（以对角元为单位提出倍数），可以把U变成单位上三角矩阵。

图源水印

既然有了LDU分解， 那么结合LL分解，LDL分解就比较好理解了。其实LDL分解也能看成LL分解：

数学推导：

（观察公式发现有重复的计算：令

，如果不做代换则乘除运算量与原来差不多，因为虽然避免了开方运算，但是计算每个元时多了相乘的因子。

算法代码：

attention：在编程实现算法时，应考虑节省存储量，可将L的严格下三角元素存储在A的对角线位置上。

代码再补俺先赶进度掌握一下整体结构再抠细节（主要是也没搜到代码）
嗷嗷嗷老师已经咧到第8章了