共轭梯度法的预处理
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
x0 0, r0 b Ax0 , z0 M 1r0 , p0 z0 ,
对 k 0,1, 直到收敛完成。
k (rk , rk ) / (p k , Ap k ), xk 1 xk k pk , rk 1 rk k Apk ,
zk 1 M 1rk 1 ,
下三角矩阵，则预处理矩阵C为：
1 C ( D CL ) D 1 / 2 (2 )
其中 (0 1) 是参数，而且有 1 A CC T ( D CL ) D 1 ( D CL )T (2 )
二、算法
（1）共轭梯度法算法 令 x0 0, r0 b Ax0 , p0 r0 对 k 0,1, 2,直到收敛完成
ຫໍສະໝຸດ Baidu
C=chuli1(A); A1=C*A*C; b1=C*b;
ep=10^(-2); [x,n,tol]=cg1(A1,b1,ep,C);
function C=chuli1(A) D=diag(diag(A)); C=D^(-1/2);
C=chuli(A,0.009); A1=inv(C)*A*inv(C'); b1=inv(C)*b;
SSOR分裂
x= 7.85971307544587 0.422926408295008 -0.0735922390240462 -0.540643016894626 0.0106261628540363 m= 4 tol = 6.567520686375984e-05
方法
迭代 次数 5
xk [7.85971307544587,0.422926408294999 -0.0735922390236594,-0.540643016894632 0.0106261628540375] [7.859713075445865;0.422926408295008; -0.073592239024046; -0.540643016894626;0.010626162854036] [7.859713075445860;0.422926408295008; -0.073592239024046; -0.540643016894626;0.010626162854036]
方法 共轭梯度法 迭代 次数 5 xk 误差
[1;0.200000000000000;0.030000000000001;0. 3.8959661649 003999999999562;5.000000000000054e-04] 39117e-06 [1;0.200000000000000;0.030000000000000;0. 2.7755575615 004000000000000;5.000000000000000e-04] 62891e-17
~
~ ~
~
一、知识回顾
（3）预处理共轭梯度方法 由于预处理矩阵的构造方式很多，所以就形成了多种预处理 方法，总的来说，大致可归为如下几种途径： 1、取预处理矩阵（也称预优矩阵）为A的一个小带宽部分 （如三对角或对角线部分）； 2、通过A的分裂，尤其是线性稳定迭代法中的矩阵 A的分裂 构造预处理矩阵（如矩阵分裂法）； 3、通过A的各种近似分解得到预处理矩阵（如不完全分解）； 4、通过矩阵A的多项式构造预处理矩阵； 5、子结构、区域分裂、EBE预处理途径。
x= 7.85971307544587 0.422926408294999 -0.0735922390236594 -0.540643016894632 0.0106261628540375 n= 5 tol = 6.745229690695623e-07
（2）预处理共轭梯度法代码
以对角预优矩阵为例 clc;clear A=[0.2 0.1 1 1 0;0.1 4 -1 1 -1;1 -1 60 0 -2;1 1 0 8 4;0 -1 -2 4 700]; b=[1 2 3 4 5]‘; %SSOR此处程序改为
对角预优矩阵
1
SSOR分裂
1
[1;0.200000000000000;0.030000000000000;0. 004000000000000;5.000000000000001e-04]
1.2688408343 39045e-17
五、参考文献
【1】蔡大用，白峰杉.现代科学计算.北京：科学出版社，2000 【2】Richard L. Burden,J. Douglas Faires 著.数值分析.冯烟利，朱海燕译.北 京：高等教育出版社，2005 【3】曾喆昭，黄创霞，周富照.数值计算方法与应用.北京：科学出版社， 2013
误差 6.7452296906 95623e-07
共轭梯度法
对角预优矩阵
4
6.5675206863 75984e-05
SSOR分裂
4
4.7317536575 62764e-04
此矩阵并没有显示出预处理的优势，所以我们又取矩阵 A=diag([1 10 100 1000 10000])，输出结果如下
k (rk , rk ) / (p k , Ap k ), xk 1 xk k pk , rk 1 rk k Apk ,
(rk , rk ) k 1 , (rk 1 , rk 1 ) pk 1 rk k 1 pk ,
结束。
二、算法
（2）预处理共轭梯度算法 令
function C=chuli(A,w) D=diag(diag(A));L=tril(A,-1); C=(D-w*L)*D^(-1/2);
function [x,n,tol]=cg1(A,b,ep,C) r=b;p=r;n=0;x=zeros(length(b),1); while n<10 alpha=(r'*r)/(p'*A*p);r1=r;
Span ( p1 , p2 ,, pk )
当 k n 时， xn 一定是方程 (1)的精确解，这是共轭梯度法的出 发点。
一、知识回顾
2、预处理共轭梯度法 （1）预处理共轭梯度法的引进 共轭梯度法理论上是一种直接方法，即按照这个算法得 到的 xn 应当是方程 Ax b 的精确解。但由于其数值的稳定 性不好，实际上由于舍入误差的影响 xn 不可能满足原方程。 如果把它看成一种迭代法，当A条件数较大时，它收敛的很 慢。因此引进了改进共轭梯度法的手段——预处理共轭梯度 法。
共轭梯度法的预处理
设 A 为 n 阶对称正定矩阵, 考虑线性方程组
Ax b
其中x 是未知向量,b是右端已知向量。
(1)
一、知识回顾
1、共轭梯度法 如果我们对于固定的k ，在 R n 空间中取 k 个线性无关的向 量
p1 , p2 ,, pk
使 xk 满足
(x k ) min (x)
k 1 (z k 1, z k 1 ) / (rk , z k ),
pk 1 zk 1 k 1 pk
结束。
三、代码
（1）共轭梯度法代码 • function [x,n,tol]=cg(A,b,ep) • r=b;p=r;n=0;x=zeros(length(b),1); • while n<10 • alpha=(r'*r)/(p'*A*p);r1=r; • x=x+alpha*p; • tol=norm(r); • if norm(r)<ep • break; • end • r=r-alpha*A*p; • beta=(r'*r)/(r1'*r1); • p=r+beta*p; • n=n+1; • end
x=x+alpha*C*p;
tol=norm(r); if norm(r)<ep break; end r=r-alpha*A*p; beta=(r'*r)/(r1'*r1); p=r+beta*p; n=n+1; end
%SSOR此处程序改为
x=x+alpha*inv(C')*p
对角预优矩阵
x= 7.85971307544586 0.422926408295008 -0.0735922390240463 -0.540643016894626 0.0106261628540363 m= 4 tol= 4.731753657562764e-04
一、知识回顾
对角预优矩阵法 如果对称正定线性方程组( 1) 的系数矩阵 A 的对角元素相差 较大, 则可取预优矩阵C-1为
C 1 D1 / 2 , D diag (a11, a22 ,, ann ), A (aij )nn
这是最常用的预处理方阵
SSOR分裂
T 将A分裂为 A D CL CL ，其中D diag (a11, a 22,ann) , C L 为严格
一、知识回顾
（2）预处理共轭梯度法原理 将Ax=b变形为C-1Ax=C-1b。 ~ 令 A C 1 AC T , ~ x CT x, b C 1b 则
Ax b A x b
要求 cond2 ( A) cond2 ( A) 1 cond ( A ) K ( A) 2 当 A 是对称正定矩阵时， n 其中 1 是A的最大特征值，n 是A的最小特征值，K(A)是A的条件数
对 k 0,1, 直到收敛完成。
k (rk , rk ) / (p k , Ap k ), xk 1 xk k pk , rk 1 rk k Apk ,
zk 1 M 1rk 1 ,
下三角矩阵，则预处理矩阵C为：
1 C ( D CL ) D 1 / 2 (2 )
其中 (0 1) 是参数，而且有 1 A CC T ( D CL ) D 1 ( D CL )T (2 )
二、算法
（1）共轭梯度法算法 令 x0 0, r0 b Ax0 , p0 r0 对 k 0,1, 2,直到收敛完成
ຫໍສະໝຸດ Baidu
C=chuli1(A); A1=C*A*C; b1=C*b;
ep=10^(-2); [x,n,tol]=cg1(A1,b1,ep,C);
function C=chuli1(A) D=diag(diag(A)); C=D^(-1/2);
C=chuli(A,0.009); A1=inv(C)*A*inv(C'); b1=inv(C)*b;
SSOR分裂
x= 7.85971307544587 0.422926408295008 -0.0735922390240462 -0.540643016894626 0.0106261628540363 m= 4 tol = 6.567520686375984e-05
方法
迭代 次数 5
xk [7.85971307544587,0.422926408294999 -0.0735922390236594,-0.540643016894632 0.0106261628540375] [7.859713075445865;0.422926408295008; -0.073592239024046; -0.540643016894626;0.010626162854036] [7.859713075445860;0.422926408295008; -0.073592239024046; -0.540643016894626;0.010626162854036]
方法 共轭梯度法 迭代 次数 5 xk 误差
[1;0.200000000000000;0.030000000000001;0. 3.8959661649 003999999999562;5.000000000000054e-04] 39117e-06 [1;0.200000000000000;0.030000000000000;0. 2.7755575615 004000000000000;5.000000000000000e-04] 62891e-17
~
~ ~
~
一、知识回顾
（3）预处理共轭梯度方法 由于预处理矩阵的构造方式很多，所以就形成了多种预处理 方法，总的来说，大致可归为如下几种途径： 1、取预处理矩阵（也称预优矩阵）为A的一个小带宽部分 （如三对角或对角线部分）； 2、通过A的分裂，尤其是线性稳定迭代法中的矩阵 A的分裂 构造预处理矩阵（如矩阵分裂法）； 3、通过A的各种近似分解得到预处理矩阵（如不完全分解）； 4、通过矩阵A的多项式构造预处理矩阵； 5、子结构、区域分裂、EBE预处理途径。
x= 7.85971307544587 0.422926408294999 -0.0735922390236594 -0.540643016894632 0.0106261628540375 n= 5 tol = 6.745229690695623e-07
（2）预处理共轭梯度法代码
以对角预优矩阵为例 clc;clear A=[0.2 0.1 1 1 0;0.1 4 -1 1 -1;1 -1 60 0 -2;1 1 0 8 4;0 -1 -2 4 700]; b=[1 2 3 4 5]‘; %SSOR此处程序改为
对角预优矩阵
1
SSOR分裂
1
[1;0.200000000000000;0.030000000000000;0. 004000000000000;5.000000000000001e-04]
1.2688408343 39045e-17
五、参考文献
【1】蔡大用，白峰杉.现代科学计算.北京：科学出版社，2000 【2】Richard L. Burden,J. Douglas Faires 著.数值分析.冯烟利，朱海燕译.北 京：高等教育出版社，2005 【3】曾喆昭，黄创霞，周富照.数值计算方法与应用.北京：科学出版社， 2013
误差 6.7452296906 95623e-07
共轭梯度法
对角预优矩阵
4
6.5675206863 75984e-05
SSOR分裂
4
4.7317536575 62764e-04
此矩阵并没有显示出预处理的优势，所以我们又取矩阵 A=diag([1 10 100 1000 10000])，输出结果如下
k (rk , rk ) / (p k , Ap k ), xk 1 xk k pk , rk 1 rk k Apk ,
(rk , rk ) k 1 , (rk 1 , rk 1 ) pk 1 rk k 1 pk ,
结束。
二、算法
（2）预处理共轭梯度算法 令
function C=chuli(A,w) D=diag(diag(A));L=tril(A,-1); C=(D-w*L)*D^(-1/2);
function [x,n,tol]=cg1(A,b,ep,C) r=b;p=r;n=0;x=zeros(length(b),1); while n<10 alpha=(r'*r)/(p'*A*p);r1=r;
Span ( p1 , p2 ,, pk )
当 k n 时， xn 一定是方程 (1)的精确解，这是共轭梯度法的出 发点。
一、知识回顾
2、预处理共轭梯度法 （1）预处理共轭梯度法的引进 共轭梯度法理论上是一种直接方法，即按照这个算法得 到的 xn 应当是方程 Ax b 的精确解。但由于其数值的稳定 性不好，实际上由于舍入误差的影响 xn 不可能满足原方程。 如果把它看成一种迭代法，当A条件数较大时，它收敛的很 慢。因此引进了改进共轭梯度法的手段——预处理共轭梯度 法。
共轭梯度法的预处理
设 A 为 n 阶对称正定矩阵, 考虑线性方程组
Ax b
其中x 是未知向量,b是右端已知向量。
(1)
一、知识回顾
1、共轭梯度法 如果我们对于固定的k ，在 R n 空间中取 k 个线性无关的向 量
p1 , p2 ,, pk
使 xk 满足
(x k ) min (x)
k 1 (z k 1, z k 1 ) / (rk , z k ),
pk 1 zk 1 k 1 pk
结束。
三、代码
（1）共轭梯度法代码 • function [x,n,tol]=cg(A,b,ep) • r=b;p=r;n=0;x=zeros(length(b),1); • while n<10 • alpha=(r'*r)/(p'*A*p);r1=r; • x=x+alpha*p; • tol=norm(r); • if norm(r)<ep • break; • end • r=r-alpha*A*p; • beta=(r'*r)/(r1'*r1); • p=r+beta*p; • n=n+1; • end
x=x+alpha*C*p;
tol=norm(r); if norm(r)<ep break; end r=r-alpha*A*p; beta=(r'*r)/(r1'*r1); p=r+beta*p; n=n+1; end
%SSOR此处程序改为
x=x+alpha*inv(C')*p
对角预优矩阵
x= 7.85971307544586 0.422926408295008 -0.0735922390240463 -0.540643016894626 0.0106261628540363 m= 4 tol= 4.731753657562764e-04
一、知识回顾
对角预优矩阵法 如果对称正定线性方程组( 1) 的系数矩阵 A 的对角元素相差 较大, 则可取预优矩阵C-1为
C 1 D1 / 2 , D diag (a11, a22 ,, ann ), A (aij )nn
这是最常用的预处理方阵
SSOR分裂
T 将A分裂为 A D CL CL ，其中D diag (a11, a 22,ann) , C L 为严格
一、知识回顾
（2）预处理共轭梯度法原理 将Ax=b变形为C-1Ax=C-1b。 ~ 令 A C 1 AC T , ~ x CT x, b C 1b 则
Ax b A x b
要求 cond2 ( A) cond2 ( A) 1 cond ( A ) K ( A) 2 当 A 是对称正定矩阵时， n 其中 1 是A的最大特征值，n 是A的最小特征值，K(A)是A的条件数