张瑞
中国科学技术大学数学科学学院
rui@ustc.edu.cn |
- 在实际应用中,经常会遇上大型的稀疏矩阵。
- 直接法会破坏矩阵的稀疏结构。
- 这样,直接法在应用中,就需要的存贮和的计算量。
这对于大型的稀疏矩阵非常不适用。
迭代法
需要构造一个收敛到解的向量数组。
- 首先,将变为等价的的形式。
取,其中可逆。则
- 然后,取得到迭代序列。称为迭代矩阵。
设线性方程组的根为,即
则
这样
可以看到,当时,。
定理 1.
对任意初值序列收敛,当且仅当。
定义 1.
若矩阵满足,则称为收敛矩阵
定理 2.
矩阵为收敛矩阵,当且仅当
推论 1.
对任意初值迭代格式都收敛的充要条件是。
问题. 若,点列有可能收敛吗?
推论 2.
若有,则迭代格式对任意初值都收敛。
与非线性问题的迭代方法相比较,
构造的线性方程组的迭代序列的收敛性,可以与初值的选取无关。
Jacobi 迭代
对于线性方程组
变换为等价形式
产生迭代格式
写成分量形式为
算法
1、输入系数矩阵A和向量b,和误差控制eps
2、x1={0,0,…..,0} , x2={1,1,…..,1} //赋初值
3、while( ||x1-x2||>eps) {
x1=x2; // x1 是第k步,x2是第k+1步
for(i=0;i<n;i++) {
x2[i]=0;
for(j=0;j<i;j++) {
x2[i] += A[i][j]*x1[j]
}
for(j=i+1;j<n;j++) {
x2[i] += A[i][j]*x1[j]
}
x2[i]=-(x2[i]-b[i])/A[i][i]
}
}
4、输出解x2
- 可以看到,每一步迭代是的乘除,的加减,共次的运算。
- 在每一步迭代,并不会改变中元素的值
- 若矩阵中有大量0元(稀疏矩阵),这些0元并不会随着迭代的计算而改变,
而且,0元并不需要参与运算。
- 可以看到,矩阵中的所有非0元参与1次乘除和1次加减,就是一步迭代的计算量。
- 因而,对于稀疏矩阵,一步迭代的计算量与非0元的个数一致。
将格式写成矩阵形式
迭代矩阵可以写成
若记
则Jacoib迭代的迭代矩阵可以写为
Jacobi迭代可以写成
- 迭代收敛的充要条件是。这个条件对于对一个任意形态的阶矩阵使用并不方便。
- 可以给出一些方便使用的充分条件
定理 3.
矩阵满足如下条件之一,则Jacobi迭代收敛。
- 是严格行或列对角占优阵
- 满足
证明. 迭代矩阵为
利用迭代矩阵的,可以证明2。
用迭代矩阵的,可以证明严格行对角占优时,Jacoib迭代收敛。
当严格列对角占优时,
由
即证
注.
若有,满足,则有
即与有相同的特征值
注.
为严格行对角占优阵,则的Jacobi迭代收敛。
由
可证明Jacobi迭代收敛。
例 1. 计算如下矩阵的Jacobi迭代的迭代矩阵,并判定对任意初值的收敛性
解. 矩阵的Jacobi 迭代矩阵为
特征多项式,特征值
矩阵的Jacibo迭代矩阵
特征多项式为,可得
Gauss-Seidel迭代
在Jacibo迭代中,使用新计算的分量来代替前一步的分量,得到格式
两边取极限后,可以看到仍然满足原方程
写成分量形式为
- 可以看到,计算量与Jacobi迭代是一致的。同样对于稀疏矩阵,一步迭代的计算量与矩阵非0元个数是一致的。
算法
1、输入系数矩阵A和向量b,和误差控制eps
2、x1={0,...,0}, x2={1,1,…..,1} //赋初值
3、while( ||x1-x2||>eps) {
x1=x2; //x1 是第k步,x2是第k+1步
for(i=0;i<n;i++) {
t=0.0
for(j=0;j<i;j++) {
t += A[i][j]*x2[j]
}
for(j=i+1;j<n;j++) {
t += A[i][j]*x2[j]
}
x2[i]=-(t-b[i])/A[i][i]
}
}
4、输出解x2
记
则Gauss-Seidel迭代的可以写成
进一步有
Gauss-Seidel迭代的迭代矩阵为
Gauss-Seidel迭代收敛的充要条件是。
这个条件在使用中很不方便。
定理 4.
若满足如下条件之一,则Gauss-Seidel迭代收敛:
- 为严格行或列对角占优阵
- 为对称正定阵
证明. 若为行或列对角占优阵,则它的Gauss-Seidel迭代的迭代矩阵的特征多项式为
可以看到,当时,仍然是一个对角占优阵。
这样,就有。
也就是说,若,则有。
即
例 2. 给出如下矩阵的Jacobi和Gauss-Seidel迭代的迭代矩阵
解. Jacobi迭代矩阵
特征多项式,特征值是
,。
得到谱半径为
Gauss-Seidel迭代矩阵
谱半径为
可以看到Gauss-Seidel迭代的谱半径比Jacobi迭代的谱半径小,因此收敛更快。
Gauss-Seidel迭代的迭代矩阵,可以直接由公式来计算。
或者用如下方法:
Gauss-Seidel迭代的分量形式为
将上式中的表达式代入其它的式子,
将上式中, 的表达式代入其它的式子,
最终,可以得到
取初值,计算,得到
Jacobi:
Total Steps: 25
Result: [1.0000000e+00 1.1628402e-09 1.0000000e+00]
迭代矩阵:
[[0. 0. 0.75 ]
[0.33333333 0. 0. ]
[0. 0.33333333 0. ]]
特征值: [-0.21839512+0.37827144j -0.21839512-0.37827144j 0.43679023+0.j ]
谱半径: 0.43679023236814957
---------------------------------------
Gauss-Siedel:
Total Steps: 9
Result: [ 1.00000000e+00 -9.69033742e-11 1.00000000e+00]
迭代矩阵:
[[ 0. 0. -0.75 ]
[ 0. 0. 0.25 ]
[ 0. 0. -0.08333333]]
特征值: [ 0. 0. -0.08333333]
例 3. 计算如下矩阵的Gauss-Seidel迭代的迭代矩阵,并判定对任意初值的收敛性
解. 矩阵的 Gauss-Seidel 迭代矩阵为
特征多项式,特征值
矩阵的Gauss-Seidel迭代矩阵
特征多项式为,特征值
松弛迭代
记
为迭代格式的修正量,将迭代格式看成是前一步值加上一个修正量的过程。
在修正量上乘上一个控制因子,可以更好的控制收敛的速度。
这样,就可以得到新的迭代格式,这个格式更灵活
其中的是某个格式得到的。
取为Jacobi迭代格式,则有
称为Damped Jacobi Method。
取为Gauss-Seidel迭代格式,则有
称为松弛迭代 (Successive Over-Relaxation method)。
松弛迭代算法
1、输入系数矩阵A,向量b,松弛因子omega,和误差控制eps
2、x1={0,...,0}, x2={1,1,…..,1} //赋初值
3、while( ||x1-x2||>eps) {
x1=x2; //x1 是第k步,x2是第k+1步
for(i=0;i<n;i++) {
t=0.0
for(j=0;j<i;j++) {
t += A[i][j]*x2[j]
}
for(j=i+1;j<n;j++) {
t += A[i][j]*x2[j]
}
t=-(t-b[i])/A[i][i]
x2[i]=(1-omega)*x2[i]+omega*t
}
}
4、输出解x2
定理 5.
若SOR方法收敛,则
证明. 当SOR方法收敛,则,所以有
另外,
即有
定理 6.
若是对称正定阵,则SOR方法当时收敛。
证明. 设是的任一特征值,是对应的特征向量,则
两边对求内积
可得
由是正定的,则,即是正定的,且。
记,则有
则
当时,
即。因此SOR方法收敛。
注.
类似,可以证明,若对称正定,Jacobi迭代收敛 正定
程序作业
编写Gauss-Seidel迭代与松弛迭代的程序,并且计算
给出迭代步数及最佳松弛因子。
- 取初值向量为,误差阈值
- 打印出根
- 给出Gauss-Seidel迭代需要的迭代步数
- 对,给出SOR迭代需要的迭代步数
输出:
根为:
0.03
0.10
...
Gauss-Seidel迭代的迭代步数: 100
SOR的迭代步数:
0.02, 300
0.04, 320
...
最佳松弛因子是: xxx
例 4. 对线性方程组
- 写出Jacobi迭代和Gauss-Seidel迭代的分量形式
- 判断(1)中两种方法的收敛性,如果均收敛,说明哪一种收敛更快
- 取初值,用Gauss-Seidel迭代法计算2步
本节读完
例 5. 对线性方程组
- 写出Jacobi迭代和Gauss-Seidel迭代的分量形式
- 判断(1)中两种方法的收敛性,如果均收敛,说明哪一种收敛更快
- 取初值,用Gauss-Seidel迭代法计算2步