（计算数学专业论文）矩阵拟谱计算的若干研究.pdf

摘要 矩阵a 的c - 拟谱是指在2 范数度量下，所有扰动量小于e 的复矩阵特征值组 成的复平面上的一个子集。矩阵拟谱能够给出在扰动影响下谱的一个全局概观， 当矩阵a 具有高度的非正规性时，那么它的一个很小的扰动能够导致矩阵谱的 一个很大的扰动。相对矩阵谱，拟谱能够给出更多的信息，但是计算拟谱是困 难而且高代价的。本文对近来拟谱的研究发展作了一个全而的描述，进而在拟 谱的计算上，既给出了采用格子s v d 方法计算矩阵拟谱的最佳矩形界定区域算 法，也给出了可控精度的自适应的拟谱边界曲线跟踪算法，这两个算法都是新 的方法。本文首次在曲线跟踪方法中借用三次h e r m i t e 插值的方法，从而细化了 拟谱区域的边界曲线。文章给出了充足的数值实验，结合图片，图表给出了详 尽的说明。数值实验表明，文中的的算法及策略是非常有效的。 关键词：拟谱，曲线跟踪，非正规阵，s v d ，矩阵扰动，h e r m i t e 插值 a b s t r a c t t h ee - p s e u d o s p e c t r u mo fa m a t r i xai st h es u b s e to ft h ec o m p l e xp l a n ec o n - s i s t i n go fa l le i g e n v a l u e so fc o m p l e xm a t r i c e sw i t h i na d i s t a n c eeo fa ，m e a s u r e d b yt h eo p e r a t o r2 - n o r m s ot h ee p s e u d o s p e c t r u mg i v e sag l o b a lp e r s p e c t i v eo n t h ee f f e c t so fp e r t u r b a t i o n s w h e na m a t r i xh a sah i g hn o n n o r m a l i t y ，t h e na s m a l l p e r t u r b a t i o no na c a nr e s u l ti na l a r g ep e r t u r b a t i o no ft h es p e c t r u m c o m p a r i n g w i t hs p e c t r u m p s e u d o s p e c t r u mc a ng i v em o r ei n f o r m a t i o nt h a ns p e c t r u md o e s , b u ti ti sd i m c u l ta n de x p e n s i v et oc o m p u t e i nt h i sp a p e r ，w eg i v eac o m p r e h e n s i v ed e s c r i p t i o no nt h el a t e s tr e s e a r c hw o r ko fp s e u d o s p e c t r u m n e x tw e f o c u so n t h ec o m p u t a t i o no fp s e u d o s p e c t r u m ，a n dg i v eap e r f e c tr e c t a n g u l a ra r e aw h i c h c o n t a i nt h ec o r r e s p o n d i n gp s e u d o s p e c t r u m ，a l s op r e s e n ta na d a p t i v ep a t hf o l l o w - i n gm e t h o dw h i c hc a nc o n t r o lt h ec o m p u t a t i o np r e c i s i o n t h et w oa l g o r i t h ma r e a ui l e wm e t h o ds of a ra n dw ef i r s tb r i n gf o r w a r dt h eh e r m i t ei n t e r p o l a t i o ni n t h ep a t hf o l l o w i n gm e t h o dw h e nc o m p u t i n gp e s u d o s p e c t r u m 、v eg i v ea b u n d a n t e x a m p l e sw i t hp i c t u r ea n dt a b l e sa n de x p l a n a t i o no nt h e m e x p e r i m e n t ss h o w t h a tt h em e t h o d sa n da l g o r i t h m sw eg i v ea c h i e v ee m c i e n c y k e y w o r d s ：p s e u d o s p e c t r u m ，p a t h f o l l o w i n gm e t h o d ，h i g hn o n n o r m a l i t y , s v d ， m a t r i xp e r t u r b a t i o n ，h e m i t i a ni n t e r p o l a t i o n 浙江大学研究生学位论文独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发 表或撰写过的研究成果，也不包含为获得逝姿盘堂或其他教育机构的学位或 证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文 中作了明确的说明并表示谢意。 学位论文作者签名：力多芽子 签字日期： 矽谚年乡月石e l 学位论文版权使用授权书 本学位论文作者完全了解逝鎏盘堂有权保留并向国家有关部门或机 构送交本论文的复印件和磁盘，允许论文被查阅和借阅。本人授权逝鎏盘堂 可以将学位论文的全部或部分内容编入有关数据库进行检索和传播，可以采用影 印、缩印或扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者签名： 加帮平 导师签名： 签字日期：动。子年月占日 签字日期： n，f 致谢 转眼已经到了写致谢的时候了。欣喜之余，却是百般滋味在心头。就像一 个长途跋涉的旅人，一路翘首，到了终点，再蓦然回首，却是感慨良多。在玉泉 生活的两年，仿佛白驹过隙，来不及回味，也来不及反省，转眼就要离开了。但 即便是这短短两年的研究生生活，也将会是一个永远融化不掉的情节。文章千 古事，得失寸心知。薄薄的文稿，凝聚了成功的喜悦，失败的茫然，还有太多人 的鼓励和期望，变得格外厚重。 一颗感恩的心，在此感谢所有给我批评给我帮助爱我关心我的人。首先衷 心感谢我的导师张振跃教授。感谢您一直以来家人般的关心和呵护，感谢您的 宽容和厚爱，也感谢您您每次讨论班孜孜不倦的教导和在百忙之中对论文的修 改润色。您渊博的学识、严谨求实的治学态度、踏实做事诚实做事的风格让我 受益匪浅。师恩凝重，无以回报，唯有极尽所能以求不负您的期望。感谢浙江大 学图形图像研究所给我们提供了一个良好的研究环境，为我们的学习提供了便 利；感谢数学系的领导和老师在学习生活上给予的关心。感谢实验室的师兄师 姐。尤其是裘渔洋、卢俊峰、赵凌潇、何睿、曹沛林、徐丽娟等师兄师姐的指导 和关心；感谢蔡小吴和张永亮一直以来的交流，特别是这最后在一起的学习时 光；感谢师弟潘建民、赵科科的帮助和建议；感谢你们，和你们朝夕相处的两年 让我的学习研究室生活变得格外丰富。感谢我亲爱的朋友们，特别是我的室友， 和你们度过了快乐的时光。千言万语，尽在无言。感谢我的家人。特别要感谢我 的父母，感谢你t r - - - 十多年来的含辛茹苦和无尽的支持，你们的爱支撑着我一 路走过来，也将激励着我一路走下去。 我来时，桂花正香；我走时，樱花已落。在西子湖畔最美的季节里写下这一 段告别的文字来纪念这里曾经的人、曾经的事，然后开始我人生新的征程。再 次谢谢你们! 加帮平于浙大玉泉欧阳纯美楼 第一章绪论 1 1 背景知识 线性模型的传统分析都是基于特征值的，并且对于数学，自然科学，工程学 中的交叉问题，这种分析也是非常成功的，尤其对于自适应矩阵，谱分析的效果 更是无与伦比，因为它们的特征向量能够构成一组正交基。特征值技术也成功 应用于声学，建筑学分析，量子力学，流体力学，数值分析上。 当碰到一个矩阵或算子的特征向量不能构成一组正交基时，计算时必须小 心谨慎，这一点已得到广泛的认同。我们称这类矩阵或算子具有非正规性，此 时这些矩阵或算子特征值的重要性就不再那么明显，它们具有行为( b e h a v i o r ) k 的多样性| 9 1 。 给定一个矩阵a ，众所周知，有 l i ml i a n | | 1 n = p ( a ) r i 其中p ( a ) 表示矩阵4 的谱半径。特别的，当p ( a ) o ke 0 e ( 2 6 ) 第二章拟谱的研究理论与计算方法7 对于图( 1 1 ) 中的矩阵g ，这个估计是非常精确的：当e = 1 0 - 8 时，有p c ( g ) = 1 0 3 2 1 ，由不等式( 2 6 ) 得到的一个下界接近于3 2 1 0 6 。因此，对于非正规矩阵， 拟谱能够告诉我们矩阵或算子的瞬态性质，而谱不能。关于拟谱的其它行为及 更多的讨论可以参考2 6 1 。 正如文献【2 】中所说，只有去计算，去研究拟谱，你才能真正理解拟谱，理解 拟谱的行为。而且不难发现，对拟谱的研究大都是建立在其计算基础上的，而 且拟谱起源于数值计算，首先我们来看看拟谱计算的发展。 2 2 拟谱计算的发展 从数值计算的角度看，只有拟谱的定义( 2 4 ) 便于计算，拟谱计算离不开计 算最小奇异值。最初计算拟谱的方法为格子s v d 方法【2 j 2 ：在复平面上的选定 区域划分网格，计算各个网格点z 对应的矩阵a z ，的最小奇异值，然后根 据口。i n ( z i a ) 是否大于e 来近似地确定拟谱边界，并用m a t l a b 的c o n t o u r 函数 命令绘山。事实证明，到目前为止，格子s v d 方法也是最基本，最简单，最稳定 的方法。然而，格子s v d 方法的一个明显的缺点就是，对于大量的网格点需要重 复s v d 分解计算对应矩阵的最小奇异值，计算代价大。减少拟谱的计算量，有 如下三种途径： 1 减少矩阵阶数。 2 提高矩阵最小奇异值的计算效率。 3 界定尽可能小的拟谱区域，从而减少网格点数。 文献 3 】从减小矩阵阶数的角度出发，用a r n o l d i 迭代产生的小规模的 h e s s e n b e r g 矩阵代替原矩阵来估计拟谱。理论上看，a 特征值中并非每个对确 定拟谱都有影响力，影响最大的应该是一些外围的特征值。这就是a r n o l d i 迭代基本可行的原因，大家知道a r n o l d i 迭代更容易逼近矩阵两端的特征值。不 过这种估计是粗糙的。 文献 4 】指出，当a 具有h e s s e n b e r g 形式或是三角形式时，与稠密矩阵相 比，计算o r m i n ( z i a ) 的速度能够提高十倍。因此，在计算网格点处的最小奇异 值之前，首先把a 酉相似于这种形式，比如说采用s c h u r 分解，对于大量的网 格点，计算速度提高明显。 第二章拟谱的研究理论与计算方法8 在拟谱范围的多种估计方法中，数值域( f o v ) 方法是目前最可靠的估计 拟谱区域的方法。然而数值实验表明，就是这样估计出来的区域也是过大的。 上一章中已提到拟谱就是复平面上的些闭区域的集合，如果能够得出拟谱 的边界曲线或是连续的估计拟谱边界上点，就省去不必要的网格计算。对此，文 献【3 0 】给出了一种曲线跟踪的方法，仅仅画出拟谱的边界曲线，其中边界曲线定 义如下： a a 。( a ) 一 z ：o m i n ( z i a ) = e ，z c )( 2 7 ) 缺陷跟踪算法的计算量与跟踪步长密切联系，文献1 0 1 根据上一步的步长和当 前预估点的误差来估计下一个预估点的步长，即得到了变步长的曲线跟踪的方 法。而文献 1 1 】中采用了格子的边界跟踪法，即在拟谱边界附近采用格子方法， 不过在格子移动上的判断过程比较繁琐。 文献 2 】从算子或矩阵的谱，特征向量出发引出了拟谱。拟谱起源于数值计 算。文中对拟谱的发展及在拟谱方面作出的工作给出了详细的概述，是一本初 始了解学习拟谱的很好的一篇读物。在拟谱计算方面，文中还提到了区域排除 法：即当计算得到一个网格点处对应矩阵的最小奇异值时，通过比较此奇异值 与的大小，我们可以判断此网格点周围一个局部区域实在c - 拟谱内，还是在e 一拟 谱区域外。后来这方面的文章有【2 9 1 。接着文章考虑把原矩阵投影到低维空间， 用所得的低维矩阵的拟谱估计原矩阵的拟谱。最显著的为a n o r d i 方法，当然这 在近似估计局部区域的拟谱上还是可行的。如果预先把原矩阵三角化，实验表 明，在计算最小奇异值时可以显著的提高计算效率。计算最大奇异值已有很成 熟的方法，对于最小奇异值的计算，自然想到应用计算逆矩阵的最大奇异值方 法，即逆迭代法。文献【2 】对这些方法进行了总结，接着讨论了并行算法的拟谱 计算。当矩阵阶数上千时，k r y l o v 子空间迭代不适为一个很好的方法。然而数值 实验表明，实际计算并不理想。文中提到在迭代中应该充分利用邻近点的信息， 比较好的一个方法就是连续性方法【4 】。 前面提到，计算矩阵拟谱要计算大量的最小奇异值，特别是对于大规模矩 阵，代价是巨大的。事实上，如果采用s v d 分解，不同的网格点处计算最小奇 异值是独立的，当划分区域后，自然可以采用并行计算。迭代法在计算邻近点 的最小奇异值收敛很快，如果已知左上端的网格点处的最小奇异组，则可以同 时由两个c p u 分别计算下边和右边网格点处对应矩阵的最小奇异组，每计算依 次，都有两个分叉点。只要合理分配，就可以很好的完成计算【5 】f 6 】【7 j 。 第二章拟谱的研究理论与计算方法 9 2 3拟谱基本理论及应用 正是在拟谱计算的过程中，人们发现了拟谱的一些特殊性，它的确能够很 好的解释一些仅仅依靠谱而不易解释的现象，进而促使了拟谱在各个方向上的 发展。如流体力学，泛函分析，控制理论，微分，积分方程，量子力学，马尔可夫 链，数值线性代数等等【1 4 】 由前面我们已知，谱是拟谱的特例。文献f 1 5 1 从熟识的特征值定理出发，阐 述了相应的拟谱定理，每一个特征值定理都是相应的拟谱定理中当e = 0 时的 特例。通过这1 6 个推广的拟谱定理，我们更加深入的了解到拟谱的重要性，即拟 谱具有比谱更多的信息。 对于非正规阵，文献f 1 6 1 较深入的探讨了拟谱的增长，分析了拟谱计算中拟 谱部分坍塌的临界情况以及坍塌( 两个或多个拟谱部分从分离到连接合并的情 形) 前拟谱的增长，其中也分析了拟谱计算中的一些有多重特征值问题。 1 7 】一 文继续讨论了用矩阵拟谱来分析幂增长中的问题，特别的，分析了两矩阵有相 同的拟谱时，由已知的结论1 2 i a i i b i i 2 ，进而探讨| i a n i i l i b 竹i i ，佗2 的 变化情况。 在拟谱边界曲线追踪的计算中，我们知道，当曲线越光滑，尤其是拟谱区域 为凸区域时，计算效果就会好些。 8 】告诉我们矩阵小拟谱的增长是有l i p s c i t z 性 质的，而且每个拟谱部分都是凸区域，并给出了小拟谱增长中的一些性质。事 实上，在拟谱的增长中，当两个或更多个拟谱部分结合坍塌成一个部分时，才出 现凹域。我个人认为，所有这些都是建立在实验结果的观察上，并在猜想下去 验证这些结论，这也就是说拟谱的计算要先于拟谱的理论。拟谱的引入的正是 由于其在非正规矩阵或算子的应用，对于线性算子，微分算子的拟谱研究比较 多，可参考 1 2 1 1 3 。 拟谱的应用进一步带动了拟谱的发展与研究，特别是在数值稳定性和控制 领域方面。考虑多项式矩阵：p ( 入) = a m a m + a m 一1 入m q + + a 1 a + a o ，通 过拟谱我们能够看到当其系数矩阵扰动时多项式矩阵特征值的敏感性。f 1 8 1 定 性的分析了多项式拟谱的边界及其联通部分，并给出了基于预报校正技术的连 续性计算拟谱边界曲线的方法。1 9 1 进一步探讨了多项实矩阵在扰动下特征值 的稳定性，并结合控制系统中的例子给与了验证。f 2 0 1 中的数值例子告诉我们， 当矩阵a 极度不正规时，c h e b y s h e v 多项式矩阵p ( a ) 的双纽线接近于矩阵a 的拟 谱。拟谱在控制理论中的稳定性分析中有很多的应用，考虑最简单的控制系 第二章拟谱的研究理论与计算方法1 0 统，士= a x ，其中z 为状态向量。理论上，当且仅当矩阵a 的所有特征值有非正 实部时，系统才稳定。而往往不同的初始状态或是给一个扰动时，即使矩阵a 的 所有特征值实部都小于零，也会发生不稳定状况。这时候，对于极不正规矩阵， 矩阵拟谱是一个非常有用的工具，这方面的文章有2 1 1 f 2 2 1 。另外，在工程和科 学的计算中，处理微分一代数方程系统的初值问题的一个行之有效的方法就是波 形松弛算法。也称动力学迭代法，而波形松弛算子通常是高度非正规的，同样， 应用谱概念来研究算子和迭代法的特性时就会遇到困难，拟谱在解决这些问题 时也有很好的应用【2 3 】f 2 4 j 【2 5 l 。 文献 9 】是最近的一本关于拟谱方面比较全面系统的书籍，对拟谱的历史发 展，拟谱的行为，拟谱的计算与应用及前人的一些工作给出了总结。文中作者 例举了大量的例子说明了拟谱比谱的重要性。这些在前面的介绍中都或多或少 的提到过。在文中，作者也提到，由拟谱定义2 2 ，很多人都误认为拟谱的研究就 是特征值的扰动理论。而作者一再强调拟谱最大的应用就是非正规矩阵或算子 的行为。另外，类似于广义特征值和拟谱的定义2 3 ，非方形矩阵的拟谱f 3 3 1 也纳 入到人们的研究范围。 2 4最小奇异值的计算 正如上述所说，拟谱的计算离不开最小奇异值的计算，方法有很多种， 如s v d 分解，转化为逆矩阵应用a r n o l d i 方法，迭代求解。文献 4 】特别提到应用 连续性的方法计算最小奇异值，即利用矩阵的最小奇异组的连续性，有如下算 法： 算法2 4 1 计算小规模矩阵b = z i a 的最小奇异值僧次计算时，钆是b 的左 单位奇异向量的一个预先给定的近似向量) ，作矩阵b b h 的逆迭代 1 估计矩阵j e 7 的右奇异向量：v b _ 1 让，通过求解b v = u 得到。 2 0 1 = i l l v l l ，v = v l l v l i 3 估计矩阵b 的左奇异向量：u b 日钉，通过求解b h u = u 得到。 彳0 2 = 1 1 1 让1 1 ，u = u l l u l l 5 当满足i 盯l 一0 2 l s i g c r 2 条件n 寸( s i g 蔓s 设定精度夕，令最小奇异值盯= 1 l i b _ 1 u l l 并退出计算，否则返回依次执行上述步骤。 第二章拟谱的研究理论与计算方法 当矩阵规模较大时，逆迭代法计算效率就会降低。假如我们计算得到乏，一 a 的最小奇异值和相应的奇异向量面，现在要计算矩阵b = z i a 的最小奇异 值，这里z 接近于乏。则可以考虑对矩阵b b n 应用初始向量为云的k 步的l a n c z o s 算 法f 1 1 。因为硅要近于乏，则k 维的k r y l o v 子空间在b 的最小左奇异向量上就会非常 富裕。用矩阵语言，我们有 v h b b h v ：t 这里y 为仇忌的酉阵( 列为l a n c z o s 向量) ，t 为k k 的三对角的实矩阵( 忌 0 和y r ，则矩阵( a 一( z + i y ) i ) 有一个奇异值e ，当且仅 当硅( ) 非空，并且z 五( 秒) 。 证明矩阵a 一0 + i y ) i 有一个奇异值e 当且仅当a 。( z ，y ) 奇异。令 j = 埘 由于 j a 。扛，可，= a - ( x + 匆) j a 何一- 。z e i 一诃，j 】= e c i 可，一z j ， 因而等价地，b , ( i y ) 一x i 奇异。故c 。( i y ) 有实特征值z 。 _ 进一步，我们可以证明如下定理： 定理3 3 给定 0 和y r ，如果咒( 可) 非空，则存在z ，使得o m i n ( a 一( z + i y ) i ) = e ，并且z 五( 可) 。 证明任取z o 置( 可) ，由定理3 2 可得，矩阵a 一( z o + i y ) i 有一个奇异值e 。 如果口m i n ( ( z o + i y ) i a ) = e ，结论自然成立。如果n ( ( z o4 - i y ) i a ) m a x ( 1 l a l t ，i i c l l ) 时，( z ) e ，因而必有z 使得 ( z ) = e 。再次应用定理3 2 可知：z 五( ) 。 类似地，给定z r ，按照定理3 2 中的方法分离y ，我们有下述类似的定 理。 定理3 4 给定e 0 和z r ，则矩阵a 一( z + i y ) i 有一个奇异值e ，当且仅 当k ( z ) 非空，并且匆k ( z ) 。 定理3 5 给定e 0 和z r ，如果k ( z ) 非空，则存在可，使得仃m i 。( a 一( z + i y ) i ) = e ，并且i y k ( z ) 。 3 3最小覆盖矩形的算法 我们将用两分法计算该最小矩形区域q + = 【a ，b + 】【c + ，矿】。其基本想法 是：首先确定人f ( a ) 的一个矩形估计q o = 【a o ，b o 】【c o ，d o 】，不妨就假设这个矩 第三章拟谱范围的矩形界定1 5 形预估计即为q f o y ，然后单边地逐次校正矩形的四个端点值。下面以逐次校正 a o 为例，说明其具体过程。 令q o = m i n r e ( a ) ia a ( a ) ) ，显然：a o a 口o ，并且展( q o ) 有纯虚 特征值。如果展( 知) 有纯虚特征值，结合定理3 3 得a 4 = a o 。不妨设反( q o ) 没有 有纯虚特征值，令z = a o + ，n o 。如果鼠( z ) 有纯虚特征值，则令a 1 = x ，o 1 = q o ， 否则a 1 = a o ，q 1 = z 。一般地我们有f a 七，a 老】，使得鼠( n 七) 没有但鼠( q 七) 有纯虚 特征值考虑b ( 堡学) 的纯虚特征值。并令 陬扎扣 搿 如果最( 篮笋) 没有纯虚特征值 如果b e ( 监笋) 有纯虚特征值 根据有限覆盖定理， a + = l i m a k 类似地，可以证明b = l i mb k = l i m 像，其中阮= m a x r e ( a ) l 入a ( a ) ) c o ，晶= m a x i m ( , k ) l 入入( a ) ) d o ， 酬= ：辫2 】 如果g ( 学) 没有实特征值 如果g ( 学) 有实特征值 以及 g k + l , d k + l ，= 藩苗1 翥蓁薹 攀凄喜冀基 3 4 算法的收敛性与稳定性分析 由于算法采用的是二分法的思想，显然算法是收敛的。故只需要证明 最佳覆盖矩形算法的稳定性，不妨只证明左边界a 七是收敛的。回到算法本 身，当确定了上下界( n 七，口d 时，当反( 盐产) 有纯虚特征值时，则由定理3 2 可 第三章拟谱范围的矩形界定1 6 得，8 + 【a k ，盟笋】。而当尻( 啦笋) 没有纯虚特征值时，假设o + a k ，啦2 j 1 ， 则说明在直线z = 竖笋左侧有一个拟谱部分并且完全在此直线左边。由 于业笋 e 。于是，由矽 的连续性，有： 妒( 徇) = l i m 矽( ) e 而由锄s 知，矽( 劲) e ，因此矽( 徇) = e 。这表明矽在s 为常数，即矽( z ) = e ，z s 。这显然不成立。因此假设不成立，即有s 的内点，不妨仍记为徇，使 得妒( 约) = m i n z s 砂( z ) 。 现证明矽( 徇) = 0 。记z o i a 的相应于最小奇异值矽( 徇) 的左右单位奇异 向量分别为u m i n 和i n 若矽( 徇) 0 ，则由( z o i a ) i n = 矽( 细) u m i n 知， 妒( 劲) 一1 = 秽黑n ( z o z a ) 一1 u m i n 考虑解析函数 g ( z ) = h i n ( z z a ) u m i n 对于z s ，我们有 1 9 ( z ) i l i ( z z a ) 一1 | 1 2 = 矽( z ) 一1 矽( 知) 一1 = l g ( z o ) 1 这表明钿是解析函数g 在有界闭域s 上的一个内部最大值点，这与解析函数 的性质矛盾，因此矽( 徇) = 0 。 由于q o 为特征值的最小实部，故由上述定理得，不可能有一个拟谱部分完 全在其左边，从而不可能有一个拟谱部分完全在其左边直线z = q 七左边。 第四章自适应步长的拟谱边界曲线跟踪 4 1曲线跟踪方法 格子s v d 方法不但计算了拟谱部分外围的点，而月还计算了拟谱区域内 部的所有点，这样所需要的计算量就会非常大，因为我们需要的仅仅是拟谱的 边界，矩阵a 的e 拟谱边界o a 。( a ) 定义如( 2 7 ) 所示。我们定义函数： 砂( z ) = o m i n ( z i a ) ，名c 由此计算拟谱边界即是在复平面上寻找满足妒( z ) = e 的曲线。不难看出函数 妒( z ) 是一个实值非负函数，其零点就是a 的特征值。如果将z 表示为z = x + i y ， 则我们可以使用记号矽( z ) = 矽( z - t - i y ) = 矽( z ，可) ，关于函数矽( z ，y ) 的光滑性有 如下的定理【3 0 】。 定理4 1 令( z + i y ) c a ( a ) ，那么矽在( z ，y ) 的邻域是实解析的，如果 矽( z ，y ) 是( z + i y ) i a 的单奇异值，那么关于矽的梯度，有： v 砂( z ，y ) = ( r e ( v h m i n u 。i n ) ，j l ( q 蒜n u m i n ) ) ( 4 1 ) 这里v m i n ，u m i n 是( x + i y ) i a 的分别对应于矽( z ，y ) = i n 的左右单位奇异向 量。 我们简记曲线跟踪算法为p f 算法，描述如下： 算法4 1 1 p f 算法基本思想( p r e d i c t o r - c o r r e c t o rp a t hf o l l o w i n g ) 优计算得到o a 。( a ) 上的一个初始点名1 ，循环执行如下步骤j ，反 j 决定一个方向p k c ，且物0 = 1 ，以及一个步长靠，计算得到一个方 向m 上的预报值磊= z k 一1 + t k p k 。 2 利用预测值磊，确定一个矫正方向c k ，l l l i = 1 和一个矫正步长以来计算 矫正值z k = 磊+ 以。 第四章 自适应步长的拟谱边界曲线跟踪 1 8 图4 1 ：p f 算法示意图 p f 算法示意图如4 1 所示。预报方向的一个合理的取法是边界曲线在z 七一1 处 的切向，预报步长的选取可以随意，但为了不使预报点偏离拟谱边界曲线 太远，应取适当小，但是过小的7 忌将会增加绘制边界曲线的计算量。可以 考虑p 惫的垂直方向作为作为矫正方向，此时可以通过一维搜索矽( 氛+ o c k ) 来 确定p + ，使得矽( 缸+ o * c 七) e 。由于要求矽一e 的零点，因此也可以考虑通 过n e w t o n 方法【3 0 】来逐步得n z 惫。在预报点氛处，计算移的梯度向量作为矫正方 向，作一步n e w t o n 迭代可以得到新的预报点氛。 珏磊卅k c k 缸一茬 ( 4 2 ) 其中( u m f f m v r a i n ) 为矽在磊处的最小奇异三元组。可以重复迭代得到理 想的矫正点魂。由于我们采用逆迭代法完全可以方便得到最小奇异组，从而 上述n e w t o n 迭代非常简单，而且预测方向可选为梯度的垂直方向。而且如采 用上述一维搜索矫正点，不断增加计算量，而且繁琐，所以选择n e w t o n 法比 较自然。以下讨论中取c 詹，p k 分别为相应点的梯度方向和预测方向，数值上完 全可以用一个复数代替一个2 维方向向量，取p 七= i c k ( i = j ) 。对于步长的 第四章 自适应步长的拟谱边界曲线跟踪 1 9 选取，文献 3 0 】没有仔细讨论，数值中取t k = m i n o 1 ，吾d i s t ( z k 一1 ，人( a ) ) ) ，其 中d i s t ( z k 一】，人( a ) ) ) 表示矫正点讯一1 到矩阵a 的谱集合的距离。 在数值例子中，不难看出当拟谱边界曲线形状变化越快，按照上面方法选 择步长，则矫正点的相对误差越大，即使可以通过不断的n e w t o n 迭代得到高 精度的矫正点，但是往往相邻的两个矫正点相距较大，即连线几乎不能反映这 区间真实的拟谱边界。因此一个自然的想法就是：在曲线变化平缓的时候，预测 步长应该适当拉大，反之，步长可适当减小。文献1 1 0 1 就探测过步长自动选择的 算法，其算法中前后步长有依赖性。事实上，当计算精度确定后，步长的选择应 只依赖曲线本身，而这个本身就反映在下文中所谈到的曲线形状的变化率上。 4 2 自适应的拟谱边界曲线跟踪 由上面的推测，步长的选择应该完全依赖曲线的形状变化，而拟谱边界曲 线形状变化的快慢，在数值上可以用函数矽( z ，可) 的梯度变化率来刻画，因此我 们试图从这方面出发，希望能自适应的选择步长。 令函数h o - ) = g ( z k 一1 + t p k ) 一e ，由泰勒展式，有 忽( 丁) = ( o ) + 尼7 ( o ) 7 + 妻厂”( o ) 丁2 + o ( 17 1 3 ) i 扫- t p k 为9 ( z ) 在z k 一1 点的切线方向，所以危( o ) = 0 。当r = 亿时，舍去后面的3 阶 无穷小，有 9 ( 稚) 一e g ( z k 一1 ) 一e + 去九”( o ) 如果要求矫正点的相对精度为e 1 ，即丛掣 e 1 ，于是可知当满足下面等式的死 三i ( o ) i = ( e 2 - - e 1 ) e ( 4 3 ) 就是一个满足精度要求的最大步长。 下面考虑如何估计九”( 0 ) ，我们首先取一步长氕，得到一个临时预测点 掣d = z k 一1 + 氕m ，则 7 ( 氕) = r e ( p k 也) ，其中d k = u m h 伽u m i n ，u m i n ，v m i 竹分别 第四章自适应步长的拟谱边界曲线跟踪 为霹d j a 最小奇异值对应的左右奇异向量。注意到( o ) = 0 ，因而可以用如 下差商来来估计( o ) 。 ”( 0 ) 圭 将其代入( 4 3 ) 式，可得： 九( 氕一1 ) 一九( o )r e ( p k d k ) 氕一1 0氕一1 t k = 2( 4 4 ) 上述步长估计式中仅需计算一次最小奇异组，代价很小。当得到预估点之 后，只需要通过n e w t o n 迭代矫正预估点使之达到精度要求。算法捕述如下： 算法4 2 1 自适应的拟谱边界曲线跟踪算法 优计算得到o a 。( a ) 上的一个初始点z l ，对于k = 2 ，3 ，循环执行步 骤j ，2 。 j 试探一个预报值2 配= 讯一1 + 彳 k p k ，计算2 k d d 处的奇异三元组，按照等 式似纠重新选择步长，计算新的预测点氛。 2 由公式“砂汁算矫正点，循环e 叫亡伽迭代直至达到精度要求，七= k + 1 尽管( 4 4 ) 式中预估涉及到前一步步长死一1 ，但这并不说明t k 依赖于一1 。 实际上，仅一1 的不同将导致r e ( p 南d k ) 的改变，而商素不会有较大的改变， 这是因为它是h ”( 0 ) 的一个很好的估计，而h ”( o ) 仅与边界曲线的形态有关。这一 点与3 0 1 有很大的不同，所以说步长的选择是自适应的。后面的实验结果表明， 上述算法是可行的，步长的选择是合理的，为了画出更精细的拟谱边界曲线，邻 近矫正点间的越近越好，正如前面所说，步长过短计算复杂度剧增。如下我们 试图用插值的办法来寻找曲线s 来代替端点为z 七，z 后+ 1 的线段。 4 3三次h e r m i t e 插值 一般的，在光滑性方面，样条插值是一个不错的选择，但是样条插值中，需 要解很大的方程组。事实上，在上面的算法中，我们能够知道在拟谱边界曲线 在每个矫正点处的切向量，而这正满足分段三次h e r m i t e 插值的要求。 第四章 自适应步长的拟谱边界曲线跟踪 2 1 设三次h e r m i t e 插值曲线为s ( z ) ，则s ( z ) 在每个子段k ，z 洲】都是三次式， 要求 s ( x i ) = 纨，s ( 婉) = 反， i = 0 ，1 ，礼 则s 3 ( z ) 有表达式： s ( x ) = 妒叭t x - - x i ，、玑+ 妒“t x - - x ij x 玑+ 1 + 见讥( x - 心x i 、y , + 吃矽1x - 佻x ) y ，, + 1 ，h i = x i + l x ix i z x i + i ， ( 4 5 ) 其中： ， l 妒o ( z ) = ( z 一1 ) 2 ( 2 z + 1 ) ，妒1 ( z ) = x 2 ( 一2 x + 3 ) i ( z ) = z ( z 一1 ) 2 ， 妒1x ) = x 2 ( z 一1 ) 事实上，每一个拟谱部分的边界都是封闭曲线。数值上，我们是根据切方 向计算节点处的斜率的，当已知的节点处的斜率比较大时，计算上述仇偏差较 大，而且影响插值点的计算。为了避免较大的偏差，在这区段上，我们把曲线看 成z 是可的函数，这样，z = 雪( y ) 有形式： ( ) 讪( 警) 咖( 警) 坼，地咖( 警) 机妒，( 警) ， ，h i = 阢+ 1 一犰y i y 玑+ 1 ， ( 4 6 ) 通过式4 5 ，4 6 我们能够很好的解决数值上潜在的过大偏差。正如三 次h e r m i t e 插值局部保形一样，在数值实验一章，我们会对插值效果进行较详细 的分析。 第五章数值实验 5 1最小矩形区域界定 这里我们选取的是k a h a n 矩阵，其形式如下： k n = 1 一cc 一c c ss c - - s cs c s 2 - - 8 2 c- - 8 2 c ；j 8 n 一2 - - 8 n 一2 c 8 n 1 其中，8 n 一1 = o 1 ，8 2 + c 2 = 1 。我们考虑n = 3 2 ，当e = o 1 时，蚝2 拟谱的最佳 矩形包含区域如图5 1 中小矩形所示，作为对照，图中外围的矩形是根据传统 的f o v 方法界定的。 012 图5 1 ：矩阵k a h a n ( 3 2 ) 的拟谱的矩形区域界定比较图 3 2 1 o 叫 t 喝 第五章数值实验 由图5 1 可以看出，界定出的矩形区域是最小的包含拟谱的矩形区域，而且 在面积上，我们的方法大约是f o v 方法界定的区域的1 5 ，从而对于格子s v d 方 法来说，计算量可以降低4 倍。其实，也并不是所有的矩阵都有这么好的效果， 譬如说对于下面将要提到的矩阵g 6 4 ，改进的矩形面积仅为f o v 方法界定的矩 形面积的1 1 倍。从另一方面看，计算的过程也就是给出拟谱的最大最小横纵坐 标，这在拟谱的分析上也是重要的。另一方面，我们曾想到用矩形束来包含拟 谱区域，用这些矩形的相交区域的边界来估计拟谱区域，然而这样的方法在效 率并不可行。 5 2自适应的曲线跟踪实验 由拟谱的性质知，当矩阵a 为实矩阵时，a 的拟谱为关于实轴对称，为了表 示一般性，我们还是采用绘制整个拟谱曲线的方法。与f 1 0 1 中不同的是，我们没 有用奇异值分解的方法计算最小奇异组，而采用更快的逆迭代法来计算。为了 统一，我们取e 2 = 0 3 ，而且在计算时，每计算得到下一步迭代步长h 七时，我们 令h k m i n ( h k ，丢d i s t ( z 七一1 ，a ( a ) ) ) ，从下面的数值实验可以看出，这样选取还是 有必要的。 这里我们选取上面提到的k a h a n 矩阵和6 4 阶的g r c a r 矩阵g 6 4 ，其形式 如下： g 6 4 = 111 111 一l 1 1 1 1。 1 1 1 11 11 在下面的数值例子中，计算最小奇异值组采用了逆迭代法【4 l ，图( 5 2 ) 画出了矩 阵风4 和g 6 4 的e = 1 0 ，1 0 ，1 0 ，1 0 5 各四条拟谱边界曲线。 第五章数值实验2 4 k a h a n ( 6 4 )