（计算数学专业论文）二次参数方程组线性化处理方法.pdf

2 0 0 6 年上海大学硕士学位论文 摘要 本文主要讨论二次参数方程组( a 2 a + a b 十e ) z ( ) = ，的数值求解方法，二 次参数方程组在许多实际应用领域中经常出现，如求解p d e 问题，控制论，结 构力学，q c d i 司题等等。因此如何建立二次参数方程组的有效数值方法是有 重要意义的。特别是当二次参数方程组的系数矩阵阶数很大，且参数的个数 很多时( 一般可达几百个) ，这种重要性就更突出了。 本论文比较全面地论述了二次参数方程组的有关问题。我们介绍了目前 计算二次参数方程组的三种主要方法，线性化方法，级数展开方法和投影方 法。而我们的工作主要是讨论和分析各种不同线性化模式对计算所可能造成 的影响，并结合数值实验加以说明。由于线性化方法最后归结到位移方程组 的计算，因此本论文对位移方程组的数值方法也给予了一章来介绍。全文可 以分成三部分，首先介绍了二次矩阵的谱分析和求解二次参数方程的三种方 法，第二部分着重介绍了位移方程的构造方法，最后详细介绍了二次参数方 程组的不同线性化模式处理，并给出数值实验的结果。 关键词：二次参数方程组，线性化方法，k r y l o v 子空间方法，位移方程组 2 0 0 6 年上海大学硕士学位论文 i i a b s t r a c t i nt h i st h e s i s ，w ed i s c u s st h en u m e r i c a lm e t h o d sf o rs o l v i n gq u a d r a t i cp a r a m e t e r i z e d e q u a t i o n s ( 2 a + b + c ) n ) = f q u a d r a t i cp a r a m e t e r i z e de q u a t i o n sa r i s ei nav a r i e t yo f p r a c t i c a la p p l i c a t i o n s ，f o ri n s t a z l c e 】t h en u m e r i c a lm e t h o d sf o rs o l v i n gp a r t i a ld i f f e r e n t i a l e q u a t i o n s ，c o n t r o lt h e o r y , s t r u c t u r a ld y n a m i c s ，a n dq c dp r o b l e m s s oh o wt o d e r i v ea e f f e c t i v en u m e r i c a lm e t h o di si m p o r t a n t s p e d a l l y ，w h e nt h ec o e f f i c i e n tm a t , r i c ea x el a r g e ， a n dt h ep a r a m e t e ri sh u g e ( af e wh u n d r e d s ) ，t h ee f f e c t i v en u m e r i c a lm e t h o di sn e c e s s a r y i nt h i st h e s i s ，w er o u n d l yd i s c u s sq u a d r a t i cp a r a m e t e r i z e de q u a t i o n s w ei n t r o d u c e t h r e ep r i m a r ym e t h o d s ：l i n e a r i z a t i o n ，s e r i e se x p a n d i n ga n dp r o j e c t i o n w ec o n s i d e rt u r n i n g q p et os h i n e dl i n e a rs y s t e m s ，t os o l v eq u a d r a t i cp a r a m e t e r i z e de q u a n o n s i nt h i sp a - p e r ，w em a i n l yu s el i n e a r i z a t i o nt os o l v eq u a d r a t i cp a r m n e t e r i z e de q u a t i o n s d i f f e r e n tl i n - e a x i s a t i o nc a ni n d u c ed i f f e r e n tn u m e r i c a ls o l u t i o n s p e c i a l l y ，w ea n a l y z et h er e l a t i o no fs n l u - t i o na n dp a r a m e t e ra n dc o e f f i c i e n tm a t r i c e s n u m e r i c a le x p e r i m e n t sr e p o r tt h ed i f f e r e n te l - f e e t i v e n e s so fd i f f e r e n tl i n e a x i z a t i o n t h i st h e s i se a x lb ed i v i d e di n t ot h r e ep a r t sf i r s t l y ，w e i n t r o d u c et h es p e c t r a lt h e o r ya t ) o u tq u a d r a t i cm a t r i xa n dt i n e ek i n d so fm e t h o d su s e dt o s o l v eq u a d r a t i cp a r a m e t e r i z e de q u a t i o n s e c o n d l yw ei n t r o d u c ek r y l o vs u b s p a c em e t h o d s t os l o v es h i f t e dl i n e a rs y s t e m s f i n a l l yw ei n t r o d u c et h o s em o d e st ol i n e a rq u a d r a t i cp a - r a m e t e r i z e de q u a t i o n s n u m e r i c a le x p e r i m e n t ss h o wa l lr e s u l t s k e yw o r d s ：q u a d r a t i cp a r a m e t e r i z e de q u a t i o n ，l i n e a r i z a t i o nm e t h o d ) k r y l o vs u b s p a c e ，s h i f t e dl i n e a rs y s t e m s 原创性声明 本人声明：所呈交的论文是本人在导师指导下进行的研究工作。 除了文中特别加以标注和致谢的地方外，论文中不包含其他人已发表 或撰写过的研究成果。参与同一工作的其他同志对本研究所做的任何 贡献均己在论文中作了明确的说明并表示了谢意。 签名：丝蕴日期：2 竺5 ：! ：! 竺 本论文使用授权说明 本人完全了解上海大学有关保留、使用学位论文的规定，即：学 校有权保留论文及送交论文复印件，允许论文被查阅和借阅；学校可 以公布论文的全部或部分内容。 ( 保密的论文在解密后应遵守此规定) 签名：盔盘叁导师签名：i 煎楚童日期：型：丝： 2 0 0 6 年上海大学硕士学位论文 第一章引言 在大量的科学与工程计算中，如计算流体力学、统计学、结构工程、量子物 理、化学工程、经济模型、航空航天工业、水利发电、气象预报、集成电路模拟、 信号处理和控制、网络排队的马尔克夫链模拟等众多领域中，见 1 ，9 ，2 7 ，2 8 ，3 6 1 ， 经常要遇到二次参数方程组问题 ( ；a + a j b + g 扣= z = z ( b ) ， j = l ，s ( 1 1 ) 的数值求解，其中a ，b 和c 为n n 实( 或复) 矩阵，复参数a j 可在一个比较大的范 围内取值，且这些参数可达几百个，对每一个参数值b 都对应一个解z ( ) 。我 们希望能有效地把所有的解同时计算得到。 上述问题的数值求解的理论研究、算法研制和软件开发是当今计算数学 和科学与工程计算研究领域中的热点之一，国际上的研究工作比较活跃，对 于国内这一部分的工作还不多。我们知道，在实际问题中，矩阵的阶数往往可 以达到几千阶、几万阶、甚至几百万阶，其计算量十分可观。因此，提供一种 有效的方法将十分有意义。 处理二次参数方程组( 1 1 ) ，目前主要有三类方法，足线性化处理方法【9 j ； 二是幂级数( 或有理数) 展开方法【25 ；三是直接投影方法【1 】。 所谓线性化处理方法，即把( 1 1 ) 转成关于参数为线性的方程组 ( k + a j m ) z = d ，j = l ，s ( 1 2 ) 其中，例如 k = ( 三：) ，m = a 一) ，z = ( ：) ，a = ( ：) 这是- - 2 n 阶的方程组，规模比原问题扩大了一倍，且矩阵k 有可能丧失矩 阵a ，b ，c 的结构。然而方程组( 12 ) 的好处在于其可以转换成下列位移方程 2 0 0 6 年上海大学硕士学位论文 2 缮 ( t + j ，) 2 = d ， j = i ，s ( 1 , 3 ) 其中，t m k m 一1 ，享一m z 。这里，我们假漫了矩路m 非奇，即a 和c 非奇f 若a 或e 非 鸯，劐褥褥裂骞效楚瑷，凳【3 甓，褒本论文孛，我霸蕊霰设a ，g 蘩奇) 。我翻知逶， 求解位移方程组有赣有效的方法，舆型的是k r y l o v 予空间方法【1 2 ，1 3 ，1 7 ，2 0 ，2 2 ，3 2 】， 以此可以弥补方程扩大一倍的缺陷；而且，目前有关不精确k r y l o v 子空间方法 豹探索雯是进一步捷化了解位移方程维斡造毽，这样也使褥二次方程维懿线 桎纯魑蠼方法更其旁可行性和有效性。 幂级数展开方滋是通过展开z ( a ) = 曼就n a o ) 4 和右端，( ) = 登， ( a a o p ( 其中抽是菜一参考值) ，进丽魄较( 舻 + a 8 + g ) 墨韪( a 一 。) ：萎五d a 。关 予 幂次颚的系数，求积以获得。妁。就方法静撬点是可戳齄疆，依赖于 酌右 端，( ) 时的情形，其不足足，( a ) 与解z ( a ) 必须在a o 处可以展开成级数，即对解和 右端的光滑性要求很赢，且对远离a o 的参数a ，级数可能丧失了收敛性。因此， 在这些 瓣遥必绥选箨获瓣 。再送行震开。 而赢接投影方法照通过将原求的大规模矩阵问题约化为中小规模的矩阵 问题，辩通过现有的方法求解中小型问题而得到原始问题解的一种方法。首 先要寻找一个低维的予空闻，秀溶爨二次参数阉怒直接投影到这一予空闻，这 样就被约化为低阶豹二次参数阔愆。求解低酚静二次参数闽邀瑟容易很多，但 是问题愿如何在花费不太大的情况下寻找这样一个予空问成为该方法关键所 在。目翦，鸯接投影的k r y l o v 予空问法大多采用生成一个二阶k r y l o v 子空阈孵 方式泉约化豫掏蘧，其髂静过程参冤第二章第二节中直接投影方法的介缨。 对于求解二次参数方程组，以上的三种方法各有它们的优缺点。为了最 终能有效而又比较精确的得到二次方程组的解，就要对具体问题具体分析，选 舞菜一秘方法，或胃淤适当缝合冬秘方法，我爱嫒霞捷兹求黪方式。这是一个 很值得研究的问题。 我们的工作主要是讨论线性化处理方法。在把问题( 1 1 1 转化成线性参数 润蘧( 1 ，2 ) ，进两转化戏位移方程缓( 1 3 ) 数遭程中，褒操持系数落等绥意义下( 游 见第二嚣第一节) ，魉阵k 和m 有各种取法，如此形成了不同鼬线性亿模式。对 2 0 0 6 年上海大学硕士学位论文 3 于不同形式的模式，在求解位移方程组的有效性上会产生不同的结果。我们 试图从矩阵的结构性态上，如对称性，正定性，稀疏性以及参数的选取范围上 给予各模式的一定分析，并结合数值计算结果来给予说明。 我们知道二次参数方程的方法研究还远没有完善。虽然目前已有一些数 值方法可以进行计算，但其理论分析的结果还相当少，特别是二次矩阵 弛+ a 。b + c 的谱分析还没有好的结果。这是比较复杂的问题，涉及到二次特征值 问题。为此，我们在第二章的第一节中给予这方面的一些简单介绍。 我们知道，二次参数方程的线性化处理最后归结到位移方程组的计算，而 位移方程组本身在许多领域中也经常出现。为此，本文对位移方程组的数值 方法也给予了比较全面的论述。我们的论述是按右端项的不同类型进行的。 l 右端项为定常的，即右端不依赖于参数 ，由于k r y l o v 子空问的不变性，采 用k r y l o v 子空间方法求解相当有效，这些方法中包括g m r e s 方法f 1 3 ，3 2 】，f o m 方 法【3 2 ，3 4 】，q m r 方法【3 2 j 以及b i c g s t a b 方法【1 2 】等等。 2 右端项b = 6 ( ，) ，即依赖于参数 ，一般情况下，可采用种子投影方法【1 7 l 。 或在右端有结构6 ( ) = f ( ) 时，其中f 为与a 无关的矩阵，可以采用块k r y l o v 予 空问法1 2 0 ，2 2 1 。 3 右端项b = b ( a ，孔) ，i j ，即右端项不仅依赖于a 还与方程组前面的解 有关，显然只能逐个方程进行求解，此时如何进行预条件处理则显得有意义 了l ，2 z l 。 全文的结构如下：第一章是引言，比较全面的论述了本论文所要讨论的 问题以及有关方法和我们的工作；第二章是二次参数方程组的综述，主要介 绍二次矩阵的谱理论和数值方法；第三章则是位移方程组数值方法的介绍；第 四章是线性化处理模式的比较，也是我们的主要工作。 2 0 0 6 年上海大学硕士学位论文 4 第二章二次参数方程组综述 本章专门针对二次参数方程组进行了概述。首先介绍了二次矩阵的谱理 论知识，见1 3 8 1 。在第二节中介绍了解二次参数方程组目前较主要的四种算法， 包括直接法，线性化方法，幂级数展开法以及直接投影法。 s 2 1 谱理论 谱性质考虑这样的二次矩阵多项式 q ( a ) = a 2 a 十 b + d( 2 1 1 ) 其e e a ，b 和c 为n ”实( 或复) 矩阵。换句话说，q ( a ) 为关于参数a 的二次多项式。 般称这样的矩阵为 矩阵。我们将q ( a ) 的谱记做a ( q ) ， a ( q ) = a c ：d e t q ( a ) = o ， 即是q ( a ) 特征值的集合。 当0 ( a ) 的行列式对于所有的a 并不完全为零时，我们称矩阵0 ( a ) 为正规矩 阵；否则，称为非正规矩阵。所以没有特别说明下面所讨沦的矩阵q ( a ) 都为j 下 规矩阵。 由矩阵的特征多项式( d e t q ( a ) = d e t ( a ) a 2 ”+ 低阶项) 可以看出，当矩阵a 非 奇时，q ( a ) 为正规阵且有2 n 个有限特征值；当矩阵a 奇异时，d e t q ( a ) 的次数r 2 n 且q ( a ) 有r 个有限特征值，我们可为其再添加2 n r 个无限特征值。这里无限 特征值可看作与反多项式a 2 q ( a 一1 ) = f c + a b + a 的零特征值相对应。而一个 正规阵q ( ) 可能有两个不同的特征值对应同一个特征向量。为了解释这一现 象，我们可以例举这样的a 矩阵 q ( a ) = a + l6 a 2 6 a0 2 a6 a 2 7 + 10 00a 2 + i 2 0 0 6 年上海大学硕士学位论文 5 即可等价地看作 a = 060 060 001 可见，该矩阵为正规矩阵。因为 b = l一60 2 70 o00 g = ， d e t q ( a ) = 一6 a 5 + i i a 4 1 2 a 3 + 1 2 a 2 6 a + 1 0 且该矩阵有6 组特征对( 儿，z b ) ，k = 1 ：6 ，如下表所示 可以看出，有五个特征值是有限的( 它们即为q d ) 行列式的根) ，还有一个 特征值为无限的。而且，当a 1 a 2 时，却有z 1 = x 2 。这个例子阐释了这样一 个事实：如果一个正规矩阵q ( ) 有2 n 个不同的特征值，则该矩阵拥有n 个线性 独立的特征向量。另外，当q ( a ) 的系数矩阵为实阵时，0 ( a ) 的谱则关于复平 而的实轴对称，因而其特征值要么为实数要么就以共轭复数对出现。从以上 的例子及描述中，我们可以看出二次矩阵的谱分析是个很复杂的问题，远不 像一次矩阵的谱分析来得简单。这就注定了解二次参数方程也充满了复杂性。 规范形如果有两个行列式为非零常数的a 矩阵e ( ) 和f ( a ) ，使得成立 p 似) = e ( a ) 0 ( a ) f ( ) ， 其中p ( ) 与q ( ) 的维数相同，则将两矩阵视为等价。p ( a ) 行列式等于零时，q ( a ) 的 行列式也为零。且若e ( a ) 和f ( a ) 都不依赖a ，则称p ( ) 与q ( a ) 为严格等价。 s m i t h 定理【1 5 】给出了矩阵q ( ) 的一个规范等价矩阵r ( ) ，满足 q ( a ) = e ( a ) r ( a ) f ( a ) 2 0 0 6 年上海大学硕士学位论文6 其中r ( a ) 为对角矩阵，即r ( a ) = d t 叼( 8 1 ( ) ，e 2 ( a ) ，e n ( a ) ) ，且e ( ) 都为首多项 式，使得岛( a ) 可整除8 州( ) 。该对角矩阵r ( a ) 就被称作q ( a ) 的s m i 恤形或规范形， 且为唯一的。当然e ( a ) 和f ( ) 可取不同值。 形式。我们需要的就是要找一个与q ( ) = a 2 a + a b + c 等价的线性 矩阵k a m 。 q ? 三 = e c 柚c k a m 妒n ， 那么这个2 n 2 n 的线性矩阵k 一 m 就为q ( a ) 的线性化形式，其中e ( ) 和f ( ) 都是 行列式为非零常数f j | c j 2 n 2 n 的肚e 阵。显然，印( a ) 和一a m 的特征值是一致的。 而线性化形式不是唯一的，重要的是要找出好的形式能尽可能反映出q ( a ) 的 对称性和别的特殊的结构属性等等。 一般在实际运用中最常被用到的是以下两种线性化模式： m 二小； 协， 皿 蚓一- b 斗 僻， 加：【。r。j 怛1 4 其中d 可以是任意非奇的”n 矩阵。可以看出在这里只要选取 e c - ，= ( 口：? d 。1 ： ，f c a ，= 二； ， ( a 2 a + a b 十c ) x = ， ( 2 1 5 ) 这样就可以将陔式另记为a t t + b u + c x = ，。也就得到了一个线性方程组 二- - 训b 。u 卜 ：川= 2 0 0 6 年上海大学硕士学位论文7 与l l 中d 取为i 相对应。当然，也可以通过将( 2 1 5 ) 另记为a a “+ a b x + c x = ，得 到l 2 这种线性化模式。 至于l 1 与l 2 两种模式的选择，主要看矩阵a 和c 的非奇性。总的来说，d 一 般取为单位矩阵，或是单位矩阵的倍数矩阵，例如| | 刮，或0 e i i ，。 22 数值方法 下面考虑这样的二次参数方程组 ( a 2 a + a b + c ) x = f ，。= z ( a ) ( 2 2 1 ) 其中a ，b 和c 为n n 实( 或复) 矩阵，a = ，j = 1 ，2 s 。右端项是否依赖a 在下面 的各节中会分别说明。 2 2 1 直接法 在许多实际问题中，我们都会遇到形如( 2 2 1 ) 的二次参数方程。当( 22 1 ) 的 系数矩阵不太大时，我们就可以采用直接法进行计算。 首先线性化该二次参数方程组为线性方程组k a m 形式，再对k 和m 进行 一般s c h u r 分解，得到 w 7 k z = t w r m z = s 其中w 和z 皆为酉阵，t 和s 为上三角矩阵。这样，如果k a m 是以( 2 1 4 ) 的形式 得到的，且a 不是特征值时， m ) ：邛啦【t 一凋一。t 。 m ) ：q _ l ( w ( u 【jj 因此，只要计算t s c h u r 分解，我们就可以只花费o ( ( 2 n ) 2 ) 的浮点运算计算出q 一1 ( a ) ，( ) 因为t a s 是2 n 维的三角矩阵。 而当系数矩阵较大时，转化为位移方程组后，直接求解方法就显得力不 从心了，这时就可以运用目前所知的求解大型位移方程组的方法进行处理， 其中最常用的就是k r y l o v 子空间法，详见222 节。当然，除了利用线性化来求 2 0 0 6 年上海大学硕士学位论文8 解( 22 1 ) 外，还有一些别的思路也可以解决这种问题。例如当( 2 2 1 ) 的右端项依 赖于参数a 时，可以选择对z ( a ) 做关于某个值a o ( 在参数附近) 的幂级数展开的方 式来求解，见2 2 3 节；或是由j 2 4 节中将介绍的直接投影法来求解。 2 2 2 线性化处理方法介绍 我们知道解位移方程有着很有效的方法，所以就自然的考虑到将二次参 数方程换化为位移方程来求解的方法。例如，我们首先将( 2 2 1 ) 等价写为下面 的参数方程 忙小( ：圳= ( 0 ( 2 2 2 ) 这里，我们考虑右端项为定常项，即不依赖a 的情形。上式简记为( k + a m ) z = d 。 再将其转化为位移方程形式，如 ( t + a ，) 拿= d( 2 2 3 ) 其中，t = k m ，= m z 。当然还有各种不同的模式可将( 2 2 1 ) 线性化为位移 方程，具体方法参见本文第四章。 我们知道，解( 2 2 3 ) 有许多已知的方法，常见的是k r y l o v 子空间方法。这里， 我们主要运用的是a r n o l d i 过程和l a n c z o s 过程形成k r y l o v 子空问：利用a r n o l d i 过 程解决t 非对称时的情形；l a n c z o s 过程则可用于解系数矩阵t 为对称时的情形。 在用k r y l o v 子空间迭代法求解时，只要取各方程的初始解都为零，即葡( ，) = o 则 每个方程的初始残向量就相同r 0 ( a ，) = d 一( t + a i ) 0 ( a j ) = d 。从而由每个方程 产生的k r y l o v 子空间也就完全相同墨。( t 十h i ，t o ( ) ) = k m ( t ，d ) 。这样，我们求 解这些位移方程组时，k r y l o v 子空问只要由一个方程计算一次即可，大大减少 了计算量，提高了解方程组的效率。我们看到系数矩阵t 是由原问题( 2 2 1 ) 中 矩阵a ，b ，c 构成的，它的规模比原问题的规模扩大了一倍，而且也可能丧失 了a ，b ，c 原有的较简洁的结构。但是，我们还是可以利用k ，m 的特殊构造 用有效的k r y i o v 子空间方法弥补这些不足f 具体的矩阵构造将在第四章中介绍) 。 而且，正如引言中介绍的，在下面3 2 节中将介绍的不精确k r y l o v 子空间法进一 步丰富了解位移方程组的方法，提高了解题效率，也加大了线性化处埋方法 2 0 0 6 年上海大学硕士学位论文 9 的可行性。 2 2 3 幂级数展开方法介绍 我们可以将( 2 2 1 ) 记为 a ( a ) x ( a ) = f ( a )( 2 2 4 ) 考虑右端项依赖a 时的情形。在a 附近取参考值知满足| f a ( a o ) | | 0 ，分别对。( a ) 和，( a ) 进 行幂级数展开，得到： z ( ) = 瓢一知) ； ( 2 2 5 ) ，( a ) = ( a 一如) ( 2 2 6 ) 并将a ( a ) 写作如下形式： a ( a ) = ( 一抽) i a i ， 其中，a o = 硝a + a o b + c ，a l = b + 2 a o a ，a 3 = a 。则( 2 2 4 ) 就可写作如下形式 萎2 ( a a 。) a i z ( ) = ，( a ) ( 。z r ) 而a ( x ) i a 0 可被记做a ( x ) = a o + 5 a ( a ) 。因为a o = a ( a o ) ，且i l a ( a o ) | | 0 ，所以存 在筇1 。这样我们可以在( 2 24 ) 两边同d j ：乘t a 0 1 ，由于a ( a ) = a o + 5 a ( a ) 故可得 到 ( j + a 0 1 6 a ( a ) ) z ( a ) = a 0 1 ，( a ) 记a ( a ) = ，+ a 0 1 5 a ( a ) ，可见，只要 口陋i 1 5 a ( a ) 】 i( 2 2 8 ) 成立( 其e e o ( a ) 表示a 的谱残量，即定) t o ( a ) = m a 。c a 小 t 表示a 的第i 个特征值) ， 则有 五一1 ( a ) = e ( 一1 ) 2 【4 i l d 4 ( a ) 1 t ( 22 9 ) i = 0 所以，由一( a ) i i al i ( i i | i 表示任意矩阵范数) 可知，当 | | 筇1 5 a ( a ) i | 1 ( 2 2 1 0 ) 2 0 0 6 年上海大学硕士学位论文 1 0 等式( 2 2 9 ) 成立。以上不等式可用于估计。( a ) 的幂级数展开式残量的一个上界。 假定( 22l o ) 的条件是满足的，我们就可以利用下面的式子计算z ( a ) 的幂级数展 开式各项的系数。 五( a ) ( z o + x l ( 一a o ) + ) = a 0 1 ( f o + ( a a o ) + ) ( 2 2 1 1 ) 对( 2 2 1 1 ) 中的每一项进行比较，便可得到如的推导公式： 2 孔= 筇1 ，l 一a f t l a j x i _ j ，i = 1 2 ( 2 2 1 2 ) j = l 其中，q f 当i j o 时值为0 。这样就得到了解z 的幂级数展开式的幂次项系 数，解也就近似求得了。注意到，在上面的表达式中，只需要计算一个矩阵求 逆，即求筇1 。所以我们可以通过对a 。进行l u 分解来求得a i l ，减少计算的运 算量。而且l u 分解也有助于估计i i a 0 1 f | 以便估计出a o 的条件数。当然正如引言 所说，该方法可以较方便的处理右端项依赖a 的情形，但不足是，( a ) 与解。( a ) 必 须在a o 处可以展开成级数，尤其当参数a 数量较多时，参考值a o 的选取就更加 复杂。 2 2 4 直接投影方法介绍 直接投影方法的关键思想就是将原来的大规模矩阵问题直接约化为中小 规模的矩阵问题，再通过现有的解中小规模问题的方法求出问题的解。这里我 们介绍的主要是k r y l o v 子空间的直接投影法，首先寻找个合适的低维k r y l o v 子 空俐。 我们同样针对方程( 2 2 1 ) ，取，：，( a ) ：量 ”，假定解z ( a ) 在a ：o 处可以 展开为 z ( a ) = n 。 t = 0 则通过比较方程组( 22 1 ) 的两边可以得到以下关于r 的关系： 32 = r 2一 r i a 一 汁 川 r r b 口 氟卜卜 一 一 一 e e g = i | = 珊 q n 2 0 0 6 年上海大学硕士学位论文 1 1 如果右端项只包含一项，也就是说i 1 时 = 0 ，则r t 这组向量就被称为二 阶k r y l o v 向量组，它们张成的空间就是一个二阶k r y l o v 子空间，记为风= s p a n r o ， r 1 ，h 一， ，见【1 】。若s p a n q 。) = 以，也就是说q 。为空问的正交基矩阵，则 可以通过q 。将原方程约化为一个新的低维的二次参数方程 o o ( a 2 a m + a b m + c _ m ) z d ) = 6 ( ) = b i a ( 2 2 1 3 ) i = o 其e 0 a 。= q 。t a q 。，巩产q t b q 。，g k = q 磊g q 。，且缸= q 磊，l 。直接投影方法的 步骤如下所示： 1 将原二次参数问题( 2 1 ) 用以上方法约化为新的低维的二次参数方程( 2 2 1 3 ) 。 2 用已知的求解中小规模二次参数方程的方法解出方程( 2 2 1 3 ) 的解。 3 由z = q 。得到原方程( 2 2 1 ) 的近似解。 这里给出了一个定理( m o m e n t m a t c h i n g 定理) ，见【1 】，说明了近似解与真解 有同样的1 1 2 个幂次项系数。 m o m e n t m a t c h i n g 定理：设m ；是近似解的第i 项系数，也就是说( a ) = 嚣m i 小。 则对所有i = 0 ，l ，m 一1 都有m 。= n 。 可以看到直接投影方法的优点是不用扩大方程组的规模，而是直接投影 到一个低维的k r y l o v 子空间。这样可能就比线性化方法存解题过程中需要汁算 的矩阵乘向量的计算量要小，但是该方法在产生这个二阶低维的k r y l o v 子空间 时的花费又比线性化方法来的麻烦。可以说，每种方法都有它们的长处与不 足。所以，z h a o j u n b a i l l y a n g f e n gs u 1 l 也给m 了直接投影方法的一个变体，结合 了线性化的方法，其思想如下。 考虑二次参数方程组 ( a 2 a + a b + g ) $ ( a ) = f 0 + a ，1 ( 2 2 1 4 ) 假定一个新的变量z ( a ) 满足 z ( ) = c 一1 ，0 + a 4 a ) ， 这里c 1 ，0 为z ( a ) 的首顶，将( 2 21 5 ) ) 代入( 2 2 1 4 ) 得到 ( 2 2 1 5 ) a a z ( ) + a b z ( a ) + c z ( a ) = n b c f o ( 22 1 6 ) 2 0 0 6 年上海大学硕士学位论文 1 2 ，一a ( 一：，a 一：。b ) ( ：之) = ( c - 1 。 c 一- l b f o e 一。，。，) c z 卫- ， 可见，所得的方程组右端被转化为不依赖 的常右端，见【1 】。 一般，假设一个二次参数方程组可以被线性化如下 卅憾料一圳) 一。 其中蜘= c 一，0 也就是x 幂级数展开式的首项。且这里z ( a ) 的阶数不需要与。( a ) 相 同o ，、 如果( ：甚；) 是线性方程组c 。z 1 s ，的解，那么z c a ，就一定是二次参数方程 组( 2 2 1 4 ) 的解。这也就表明，f 叫柚1 的第i 项系数对应于的z ( a ) 顶部就等于z ( a ) 的 名( ) 第i 项系数n 。这样就可以通过求解线性方程组( 22 1 8 ) 的a r n o l d i 过程而得到二 2 0 0 6 年上海大学硕士学位论文 1 3 第三章位移方程组的数值方法 在上一章中我们介绍了有关二次矩阵的谱理论及求解二次参数方程组的 数值方法，而我们的工作主要围绕线性化方法来展开的。用线性化方法求解 二次参数方程组最终归结为位移方程组的求解，所以，这一章我们专门针对 位移方程组的数值方法进行论述。 3 1 三种形式的右端项 在本节中，我们将位移方程组依据其不同的右端项形式分成三种类型，分 别给予数值方法的介绍：第一小节介绍了定常右端情形；第二小节是关于右 端项只依赖于参数的情形；第三小节介绍了右端项既依赖于参数又与前面的 解有关的情形；第四小节则介绍了两种预条件处理方法。 3 1 1 定常右端问题 在控制论，q c d i h 题( q u a n t u mc h r o m o d y n a m i c s ) 以及结构动力学中，都会遇 到形如以下的位移方程组 + a ，j ) z ( ) = b ，j = 1 ，5 ( 3 1 1 ) 由于需要求解上百个位移方程，如何建立有效的数值方法是一个感兴趣的问 题。我们知道，k r y l o v 子空间 ( a ，u ) = s p a n v ，a v ，a m - 1 v ) 对于位移具有不变性，即( a ，”) = ( a + a ，”) ，这样在求解位移方程组的过 程中，我们只需要构造次k r y l o v 子空问就可以同时解出所有的位移方程，大 大减少了计算量，因此使用k r y l o v 子空问方法求解( 3 ，1 1 ) 是自然的选择。众所周 知，k r y l o v 子空间方法有长递推和短递捧两种，前者是基于l a n c z o s 过程产生的 2 0 0 6 年上海大学硕士学位论文 1 4 算法，如c g 方法，q m r 方法以及m 一对称l a n c z o s 方法等等，后者是基于a r n o l d i 过 程产生的算法，如g m r e s 方法，f o m 方法等等。 具体随来，利用a r n o l d i 过程( 设a 非对称) 可生成( a ，6 ) 的一组正交基= 【v l ， 。】，且成立- a r n o l d i 分解a v r m = v r m + l 瓦。，由此成立 ( a + a j z ) v 。= + l - y m ( ) ， 其中7 k ( b ) = h m + a j _ 仇，一i m = k ，o 】7 。现在置。( a ，b ) e e 寻找所有方程的近似解 ( 没初值为0 ) ：z 。( b ) = 蜘( ) ，使其残量 r m ( ) = b 一( a + ，) 如。( a ，) = b 一( a + j ) 蜘( ) = b 一+ 1 日m ( ) ( ) = + 1 ( 卢e 1 一再( ) ( 如) ) 与( a ，6 ) 正交：r m ( b ) 上坼。( a ，6 ) ，可得f o m 解： 9 。( ) = ( 如) 。( p e l ) ， 其中，h m ( a j ) 是 h 。( ) 去掉最后一行所得矩阵；或使得 圳z 2 m i 1 1 , o ”瓦( ) 训2 ， 可得g m t l e s 解： ( a ，) 2 ”g 矧r a i t n 。一日t n ( a j ) y l l 2 _ 我们把求解方程组( 3 1 1 ) 酗j f o m g m r e s 方法归纳成下列算法【3 2 ： 算法3 1 1 求解位移方程组 的f o m g m r e s 算法 取初值z o ( ) = o ，j = 1 ，s ，令卢= 恻1 2 ，”1 = 6 伊。 1 对l = 1 ，m 1 ) 计算 = a v f 2 ) “= ( 2 i ，w ) ，w = 一h i l v i ，i = 1 ，f 3 ) h + l ，f = 1 1 w l l 2 4 ) + l = 1 f 2 0 0 6 年上海大学硕士学位论文 1 5 2 对j = 1 ，8 ( i ) f o m ：y m ( ) = ( ) 一1 肛l ( 1 i ) g m r e s ：可玑( b ) = a r g r a “i n m1 | 卢8 1 一日讯( b ) 训2 z 。( b ) = y m ( b ) 我们知道l a n c z o s 过程是一个三项递推过程，比起前述的a r n o l d i 过程，其计 算量大大的减少。然而前提是a 对称。在有些问题中，a 可能没有对称性，但具 有某种意义下的对称，常见的是一种称之为m 一对称的对称。 设m 是对称阵，而a 满足 a t m = m a ，( 3 1 2 ) 称a 是m 一对称矩阵。在这种对称下也可以用l a n c z o s 方法求解方程( 3 1 1 ) 。相应 的m 一对称的l a z l c z o s 算法如下 2 7 1 。 算法3 1 2 求解位移方程组的m 一对称l a n c z o s 方法 取初值z o ( ) = o ，j = l ，s ，令卢= i i b f l ，v l = b z ( 其中f = b t m b ) 。 1 对l = l ，2 l 1 ) 计算蚴= a v f 2 ) 计算”；_ w l 一觑一1 牝l 3 ) 计算啦= ，m ”； 4 ) 计算w ：7 = ”；一a j q 5 ) 计算岛= i i ”? f f m 6 ) v l + l = w ；e l 2 ，g j = 1 一，s 1 ) 形成( b ) = 碥+ ，其中= t r i d i a g ( 卢t ，啦，崩) 2 ) 计算( b ) = t m ( ) 1p e l 3 ) 计算。( ) = ( ) ，其中= 口l 一u 。1 注意n l a n c z o s 在计算近似解。时，按递推修正形式，计算量可减少，但对 于位移方程组来说，则要增加许多内存来存放解向量和辅助向量，故在本文 2 0 0 6 年上海大学硕士学位论文 1 6 中采用了上述形式。 长递推方法在计算量和存储量上会受到限制，为此一般采用重启动形式， 见 1 3 ，3 0 ，a 4 。注意到残量的表达式可写成 r 。( a j ) = 0 ( 卢e 1 一日矗( 1 j ) 。( ) ) 一 。+ l h 。+ l ，。( 1 j ) e 磊v 。( ) 对于f o m 方法，( a j ) 满足= h m ( a j ) y m ( a j ) = f l e l ，故残量 r m ( a j ) = 一v m + l h 。+ l ，。( ) e 三蜘( ) ，( 3 13 ) 即对任意，( ) 都与”。+ ，保持同向，与无关。如此，在重新启动时，可 用_ 0 = + 1 形成新的k r 鲥。”子空间。( 4 + a j i ，_ 0 ) ，其仍具有关于位移如的不变 性。 对于g m r e s 方法，残量( b ) 没有了( 31 3 ) 的形式，而是与有关。如此， 用- o = ( b ) 形成新的k r g z o ”子空间时，将丧失关于位移的不变性，为此，要 修正z 。( ) = ( ) 的限制条件：f | ( a j ) f | 2 = m i n 1 5 e i 一7 k ( ) 。怯而要求 残量r 。( b ) 与r 。强行保持线性相关性， ( ) = 风( b ) ，( 3 1 4 ) 其中r t n 是算法3 11 中g m r e s 解z 。( o ) 对应的残量，与b 无关，如此可以得到( b ) 满 足的关系 9 】： c 酬，训( 黜) 一， 其中 。见下述算法31 3 ，廓则是初始残量r o ( b ) = 卢0 r 0 的相关系数。这样，在重 新启动时，可用_ o = 形成新的k r y l o v 子空i u j k , 。( a + 如j ，- 0 ) ，仍具有关于位 移a t 的不变性。 算法3 1 3 求解位移方程组的重启f o m g m r e s 方法 取初值。o ( b ) = 0 ，j = 1 ，m 固定值m ，令凤= 1 ，卢= i l b l h ，v l = 6 卢。 1a r n o l d i 过程f 见算法3 11 ) 2 ( i ) f o m ：y m ( ) = 且。( ) 。1 肚i ，j = 1 ，s ； ( 1 i ) g m r e s ：y m = a r g m 叫i n 。一日m ( o ) y l l 2 2 0 0 6 年上海大学硕士学位论文 1 7 计算：。= p b l 耳。( o ) 。，r 。= 1 + 1 w 。 扯h 一瑚方砌矾胁，( 裂) 嗍e - 3 对j = l ，s ， z m ( ) = x o ( a j ) + 9 m ( ) 4 , 老：i l r 。( ) 1 1 2 s ，删除收敛的方程，直到所有方程收敛； 否则令。o ( a j ) = z m ( ) ， f o m ：令 l = t k + l ； g m r e s ：令v l = 7 r a l l r 。忆岛= 口。( 如) 。 3 1 2 右端依赖于参数的问题 很多情况下，我们会遇到