24小时热门版块排行榜    

查看: 2372  |  回复: 50
本帖产生 1 个 数学EPI ,点击这里进行查看
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

xjw0413

铜虫 (初入文坛)

[交流] 【讨论】预调件共轭梯度法(PCG) 已有4人参与

有限元计算经常碰到大型稀疏矩阵,由于此类线性方程组通常条件数是比较大的,方程组的性态不好,所以最好用迭代方法求解,比方说是预调件共轭梯度法,但此方法在选择预调件矩阵时似乎没有一个同一的标准,大多推荐的是采用incomplete LU decomposition做为预调件矩阵。incomplete LU decomposition的计算方法似乎又有很多种。
1. incomplete LU decomposition 的计算时间应该比 LU decomposition要快速的多吧,不然直接用LU decomposition不就解出来了吗,又何必再来PCG迭代呢?
2. 采用PCG方法的前提应该是系数矩阵对称、正定吧,因为其原理是一个相当于势函数的东西取极小值。那对于非正定的系数矩阵能求解吗,我构造了几个非正定的,有的似乎是能够收敛到正确结果的。

希望各位虫用解答和讨论。
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

The life I want, there's no short cut.
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)

★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
lovibond(金币+1): 鼓励交流 2011-06-18 10:34:45
1、还有改良条件数、数值稳定性的考虑
2、系数矩阵非正定的话,迭代过程中应该出现残量上升的情况,能够求解只能算是运气,个人以为实在迭代到负方向之前程序就已经终止。非正定的情况下通常用另一种Krylov子空间迭代法——GMRES求解。
2楼2010-09-06 22:58:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)

★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
lovibond(金币+1): 鼓励耐心交流 2011-06-18 10:35:05
引用回帖:
Originally posted by xjw0413 at 2010-09-07 08:04:33:

“个人以为实在迭代到负方向之前程序就已经终止”,你的意思是说,设定判断标准,如果残量变大就终止程序并报错吗?
我试一试你推荐的那个方法。
谢谢你了!

在特定精度下终止程序实际上就是在一定的子空间内取得了解的近似。我的意思是你得到看似正确的解可能是因为这个子空间还没有包含任何系数矩阵的负曲率方向,这种情况下能得到满足精度要求的一个近似解似乎也很自然。一旦侦测到一个负的方向,迭代序列很可能因此发散。具体的表现是什么,还是需要比较严谨的推导。总之这种情况下用CG得到的结果是没有说服力的,这个不是改变终止准则就能修正的。
4楼2010-09-07 16:06:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)

★ ★
小木虫(金币+0.5):给个红包,谢谢回帖
小雨萌萌(金币+1): 谢谢回帖 2011-06-17 14:34:05
引用回帖:
Originally posted by mastergxm at 2011-06-14 04:37:26:
我有一个别人的程序,里面用到一个scale coefficient,这个系数是在系数矩阵里找到了一个最小值,然后把系数矩阵里的所有数除以这个最小值,然后进行计算,请问这是什么意思呢?盼赐教

我猜测是系数矩阵的元素太小,为了减小数值误差而做的处理。这个应该跟preconditioning无关。
6楼2011-06-15 03:32:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)

★ ★ ★
小木虫(金币+0.5):给个红包,谢谢回帖
nono2009(金币+2): 鼓励应助 2011-07-01 07:10:32
引用回帖:
Originally posted by mastergxm at 2011-06-15 02:57:22:
你说的这种可能我同意,刚在书上看到有这种减少误差的方法。
但是我又仔细看了一下, 不是简单里除以一个最小的系数。下面图1是PCG方法迭代格式,是从一本书上粘过来的。根据我对代码的反推,发现P逆应该是图 ...

通常意义下的PCG,preconditioner只要求对称正定,这里的CC(i)只要都是正数,显然P可以用作preconditioner。但是需要知道的是,对于不同的线性方程组,preconditioner只能依照系数矩阵来做具体的选择。这里为什么这么用,针对这个代码要解决的问题,你可以看看P的应用有没有达到这样的效果:
(1) 使系数矩阵的特征值聚集;
(2) 改善条件数。

因为P是个简单的对角矩阵,除非系数矩阵也近似于相关的对角矩阵,否则,我更偏向于(2)。

另外,你可以让P=I(单位矩阵),这样PCG就变成标准CG了,可以对比一下,这个P有没有改善你的计算结果。
8楼2011-06-18 00:29:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)

★ ★
小木虫(金币+0.5):给个红包,谢谢回帖
nono2009(金币+1, 数学EPI+1): 鼓励应助 2011-07-01 07:11:01
引用回帖:
Originally posted by mastergxm at 2011-06-18 02:58:08:
解答的很专业。我会试一下
我用到的系数矩阵是一个五对角矩阵,用了这个P,我纳闷的是如何求得这个P。书上讲PCG法时说这个P要采用ILU分解方法,最后取LU的乘积,没有更具体的。有没有详细讲的?

注意一点,通常选择preconditioner的时候一个常规的做法就是选择一个系数矩阵的近似,这样的话,条件预优后得到的矩阵就会接近于单位矩阵。

假定使用LU分解,然后用L和U的乘积作为P,那么实际上就是就是用系数矩阵自身作为P,亦即条件预优后的矩阵严格地等于单位矩阵。这是最理想的情况,但是,事实上这么做没有任何意义,因为已经有LU分解,那么可以直接用这个分解求解原方程组,只需要简单的向前向后代入就足够,无需使用CG这样的迭代方法。

这里使用ILU(不完全的LU分解)是因为LU分解代价太大,尤其是对于大规模线性方程组而言。ILU分解得到的L和U作乘积之后可以充当系数矩阵的近似,按照之前所说的,也就是P的一种理想选择。

你可以看看ILU分解的具体做法,Yousef Saad的Iterative methods for sparse linear systems一书里面有专门的章节介绍,有空可以看看。这是本很好的书,在作者的个人主页可以直接下载(http://www-users.cs.umn.edu/~saad/PS/all_pdf.zip)。
10楼2011-06-18 15:21:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)

★ ★ ★
小木虫(金币+0.5):给个红包,谢谢回帖
小雨萌萌(金币+2): 谢谢回帖 2011-07-20 10:25:18
引用回帖:
Originally posted by mastergxm at 2011-06-27 08:44:50:
有另外一个问题,我碰到一个4X4矩阵,是一个准三对角矩阵,即第一行都有值,其它行与一个三对角矩阵都一样。我参考的文献中说可以采用
tridiagonal equation solver  and the Sherman-Morrison formula
方法 ...

把准三对角矩阵分裂为A+u*v^T的形式,A为三对角矩阵,u和v分别是一维向量,然后直接用Sherman-Morrison formula即可。公式见http://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula

话说回来,4*4的矩阵有必要这么大动干戈么?

[ Last edited by saladin983 on 2011-6-27 at 20:10 ]
13楼2011-06-28 02:09:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)


小木虫(金币+0.5):给个红包,谢谢回帖
u 和v是显然的,你可以令u=(1, 0, 0, 0)^T,让v=(0, 0, C_1, D_1)^T。
15楼2011-06-28 15:35:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)


小木虫(金币+0.5):给个红包,谢谢回帖
引用回帖:
Originally posted by mastergxm at 2011-06-28 02:41:43:
4X4矩阵只是我举的一例子,在计算中,这是一个可调参数,也可能是20X20,我只是举了一个例子。
我看过书上介绍上面周期三对角矩阵时,计算U,V的公式。至于我说的这个第一行都有值的这种准三对角矩阵,用Sh ...

你所遇到的准三角矩阵都可以这样处理。
16楼2011-06-28 15:50:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)

★ ★ ★ ★
小木虫(金币+0.5):给个红包,谢谢回帖
nono2009(金币+3): 鼓励应助 2011-07-01 07:11:49
引用回帖:
Originally posted by mastergxm at 2011-06-28 10:46:51:
经过你的提示,我也开了眼,也试着对下图1中的两种准三对角矩阵,计算了一下u,v,但是对于图2中的矩阵也没有了办法,不过似乎图2中的这种矩阵也就不是稀疏矩阵了,用不着 Sherman-Morrison formula了吧?或者选 ...

这是常规技巧,如果你知道拟牛顿法的秩1矫正,应该就能看出来了。

最后这型矩阵似乎在哪里见过,记不得了。对于4*4的矩阵讨论稀疏性没什么意义。就这个矩阵而言,我想应该可以用两次Sherman-Morrison formula解决吧。把矩阵写成A+u_1*v_1^T+u_2*v_2^T的形式,然后逐次使用公式。
19楼2011-06-29 00:47:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

saladin983

铁杆木虫 (正式写手)


小木虫(金币+0.5):给个红包,谢谢回帖
引用回帖:
Originally posted by mastergxm at 2011-07-15 05:35:52:
有另外一个问题,想请教一下,就是下面这个推导过程,这是别人帮我推导的。
主要是涉及变上、下限积分函数的求偏导,他说这个公式是用链式法则得到,可是我大学学链式法则是关于复合函数求偏导,没有学过这个 ...

这个形式上看起来比较复杂,我倒是没有接触过。不过考虑到u也是关于x的函数,这样的结果也很自然。可否这样考虑:
在积分上下限之间取两个常数,这样可以将积分分为三段处理,不知是否能有帮助。或者把整个积分看成一个函数J(u,zeta,h),然后按照链式法则求偏导,似乎是能得到你要的结果。
23楼2011-07-19 05:42:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 xjw0413 的主题更新
普通表情 高级回复 (可上传附件)
信息提示
请填处理意见