（计算数学专业论文）关于电子结构计算中的一类非线性特征值问题.pdf

中文摘要 摘要 电子结构计算中产生的非线性特征值问题通常采用自洽场( s c f ) 的迭代算 法，并结合电荷密度混合技巧改善收敛行为。本文的第一节先描述非线性特征 值问题的形式，并回顾s c f 算法以及电荷密度混合技巧。第二节求解了一个电 荷密度混合技巧中引出的密度矩阵存在性问题。接下来应用向量值函数微分的 相关理论，第三节提出一种n e w t o n 型的电荷密度混合技巧。第四节把特征值问 题看成是非线性方程组，对其直接进行n e w t o n 迭代。数值结果表明，新方法能 很好地解决该物理模型，并且在某些方面比其他算法有优势。 关键词：s c f ( 自洽场) ；电荷密度；电荷密度混合；n e w t o n 型混合；n e w t o n 迭代 中图法分类号：0 2 4 1 7 一i i 英文摘要 a b s t r a c t s c f ( s e l f c o n s i s t e n tf i e l d li sw i d e l yu s e dt oh a n d l et h en o n l i n e a rd g e n v a l u ep r o b - l e mf r o me l e c t r o ns t r u c t u r ec a l c u l a t i o n at e c h n i q u ec a l l e dc h a r g ed e n s i t ym i x i n gi si n t r o d u c e dt oi m p r o v et h ec o n v e r g e n c ew i t hs c f i nt h i sp a p e r , t h en o n l i n e a re i g e n v a l u e p r o b l e mi ss h o w e da n dt h es c fi t e r a t i o na n dc h a r g ed e n s i t ym i x i n ga r er e v i e w e di n s e c t i o n1 ap r o b l e mo f t h ee x i s t e n c eo f t h ec h a r g ed e n s i t ym a t r i xf r o mc h a r g ed e n s i t y m i ) 【i n gi ss o l v e di ns e c t i o n2 o n t h eb a s i so f t h et h e o r yo fd e r i v a t i v eo f v e e t o r - v a l u e d f u n c t i o n ，i ns e c t i o n3a n e w t o nt y p ec h a r g em i x i n gt e c h n i q u ei si n t r o d u c e d i ns e c t i o n 4 ，t h en o n l i n e a re i g e n v a l u ep r o b l e mi sv i e w e da sn o n l i n e a re q u a t i o n sa n dn e w t o ni t e r - a t i o ni su s e dt os o l v ei t n u m e r i c a lr e s u l t ss h o wt h a t t h en e w a l g o r i t h m sw o r kw e l lo n t h ep h y s i c a lm o d e la n da r eb e r e rt h a no t h e ra l g o r i t h m si ns o m ea s p e c t s k e yw o r d s ：s c f ；c h a r g ed e n s i t y ；c h a r g ed e n s i t ym i x i n g ；n e w t o nm i x i n g ；n e w t o n i t e r a t i o n c h i n e s el i b r a r yc l a s s i f i c a t i o nn u m b e r ：0 2 4 1 7 一i 一 前言 前言 在电子结构计算的k o h n - s h a m 模型 6 ，8 ，1 0 】中，需要在波函数满足正交性 的条件下极小化总能量。使用平面波离散化和l a g r a n g e 乘子法，能量极小化问 题转化成为一个非线性特征值问题( 细节可参考 1 4 】) 。 处理该非线性特征值问题和线性特征值问题有些区别，对于线性特征值 问题的处理方法可参考 2 ，3 】，而对于该非线性特征值问题，目前一般采用 s c f ( s e l f c o n s i s t e n tf i e l d ) 的迭代算法。在s c f 的每一步，先由近似的电荷密度 产生一个矩阵，再求解这个矩阵( 该矩阵称为h a m i l t o n ) 最小的若干个特征值 以及相应的特征向量，用所得的特征对计算下一步电荷密度的近似，这样循环 下去直到收敛。然而，s c f 迭代有时候不能保证收敛。在实际应用中，s c f 每一步的h a m i l t o n 更新之前，先对当前近似的电荷密度做一定的“改善”，通 常是把当前的电荷密度跟前一步或者前几步的电荷密度做某种混合，然后再生 成h a m i l t o n ，求解线性特征值问题。该技巧称为电荷密度混合。实践表明，在 s c f 中进行合适的电荷密度混合，能有效地改善s c f 的收敛性。常用的电荷密 度混合技巧有简单混合 5 】、p u l a y 混合 9 1 、b r o y d e n 混合 4 等。 然而，关于电荷密度混合技巧合理性的相关理论比较缺乏，很多混合技巧 也有其特别的假设条件，使得该技巧的应用有所限制。 本文第一节通过引用已有的结果，介绍非线性特征值问题的形式，并且回 顾了现有的s c f 迭代算法以及电荷密度混合技巧。第二节从电荷密度混合技巧 中引出一个密度矩阵的存在性问题，并通过构造法找出了该问题的个解。第 三节我们通过应用向量值函数微分的相关理论，在一定条件下把该算法看做是 对电荷密度的不动点迭代，提出了一种n e w t o n 型的电荷密度混合技巧。第四 节，我们同样利用向量值函数微分理论，把该非线性特征值问题转化以后看做 是非线性方程组，并对该非线性问题进行n e w t o n 迭代法的探索。数值实验表 明，新方法能够有效地求解一类特殊的物理模型。 一一 论文独趔性声踢 本论文是我个人在导师指导下避符的研究工作及取褥的研究成果论文孛 除了特另g 加以标浪和致谢的地方外不包含其他入或其落机构已经发表或撰鬻 过鹃嵇究残果茭她懑恚对本砑究的癞发秘度骰妁贾献均墨在论文中作了露糖 魏声疆并表示了落意 作者按名：基圣玺日期：兰监j 工 论文使孀授投声讶 举人完全了解簸旦大学有关保留、使用学位论文的觌宠，即：学校有权保 留送交论文的复印糌允许论文被查阕耩诺阕：学校可以公布论文的全部或部 努蠹容。霹浚象蘩彩鼙、臻露或其它复爨手段保存瓷文绦密兹逡文在麓蜜嚣 遵守此规定 谗者签名：基查尘二导师签寝：杰兰! 盈一目期：鲨j 二5 ：! 谗考签名：丛芷望导师簿寝：圆兰! 型一目期：垫! ! 二望：! 第一节问题形式以及s c f 算法 第一节问题形式以及s c f 算法 本节，我们先描述非线性特征值问题的形式，然后回顾s c f 算法以及电荷 密度混合技巧。 1 1 非线性特征值问题 对于在电子结构计算的k o h n s h a m 模型【6 ，8 ，1 0 ，需要在波函数满足j 下交 性的条件下极小化总能量。通过使用平面波离散化和l a g r a n g e 乘子法，把能量 极小化问题转化成为如下的非线性特征值问题( 参看文献【1 4 】) 其中p 表示电荷密度，满足 h ( p ) x = x a k x x = i k ， p d i a g ( x x + ) ，( 1 2 ) 这里d i a g ( ) 的含义是，当z 是一个向量时，d i a g ( x ) 表示以z 元素为对角线的对 角矩阵；当a 为矩阵时，d i a g ( a ) 表示把a 的对角元取出来做成的向量。h ( p ) 为k o h n s h a mh a m i l t o n 算予，nx 阶矩阵x 是离散化的波函数，k 表示电子 个数，n 是空问自由度，札= d i a g a 1 ，k ，厶表示后阶的单位矩阵。 h ( p ) 的具体形式【1 4 茭j 剐2 j 1 三+ + e ，w l w ；+ d i a g ( 跏) + d i a g 砒 ( 1 3 ) 其中l 是离散的l a p l a c e 算予，d 。是对角矩阵，表示离子的局部势能，劬 是已知参数，s 是离散l a p l a c e 算子的逆，。( 伽) 表示交换能的导数。显 然， 工+ d 溉+ l 毗嘶是其线性部分，d i a g ( s p ) + d i a g ( u 。( 力) 是其非线性部 分。由于交换能的导数d i a g ( # x c ( p ) ) 形式较为繁琐，d 加、f 蚴叫以及每一项 前面的系数对理论分析没有大的影响，同时为了叙述方便，我们考虑如下的简 化模型 h ( p ) = l + d i a g ( s p ) ，( 1 4 ) 文中的理论分析以及相应的算法可以自然地推广到( 1 3 ) 式。 第一节问题形式以及s c f 算法 1 2s c f 算法与电荷密度混合技巧 我们用月一来表示矩阵a 的转置，用a 来表示a 的共轭转置，a ( ) 表示迭 代过程中第i 步得到的近似。 目前主要采用s c f 类的迭代方法求解问题( 1 1 ) 。s c f 迭代算法本质上是一 个不动点迭代。给定x 的一个初始近似x ( o ) ，可以得到最初的h a m i l t o n 日( 1 ) = l + d i e 。g ( s p ( x ( o ) ) 然后求解h 0 ) 最小的七个特征值对应的近似特征向量x ( 1 ) ，再形成新的 h a m i l t o n 日( 2 ) ，以此类推。在第i 步迭代检查p 卜1 ) 和p ( x ( ) ) 的差，当它们 之间差足够小的时候算法终止。算法框架如下： 算法1 1 ：s c f 迭代 0 ) 给定，最大迭代步数 。，阀值g ，x ( 0 ) ， 1 ) 初始化i = 1 ，用x ( o ) 计算得到p 嚣以及日( 1 ) ， 2 ) 当i z 。时 3 )求解线性特征值闯题日( ) x = x 舭的最小k 个特征对a 和x ( n ， 4 )计算j ，2 = d i a g ( x ( ) x ( 小) ， 5 )计算彪和p 璺的差。如果差值小于阀值，程序终止。 6 )甜1 ) = 锟， 7 )用4 1 形成日( 件1 ) ， 8 )i = i + 1 ，回到2 ) 。 s c f 迭代是求解非线性特征值问题( 1 1 ) 的基本方法。在实际应用中，s c f 迭 代往往不收敛，或者是收敛很慢。为了改善收敛性，在更新日( 件1 ) 之前把当前 得到的电荷密度锟做一些修正，期望修正为比原来“更好”，作为嚣更 新仃( 件1 ) ，这个技巧称为电荷密度混合。实践表明，s c f 和电荷密度混合的结 合跟单纯的s c f 迭代相比，有着更好的收敛性。 最常见的混合技巧是简单混合( 参考【5 1 ) p l 1 = ( 1 一) 耀+ 口矗也，0 o 1 ( 1 5 ) 应用简单混合到s c f 中得到 算法1 2 ：简单混合s c f 迭代 一2 一 第一节问题形式以及s c f 算法 o ) 给定三，最大迭代步数 。，阀值，x ( 0 ) ，a 1 ) 初始化i = l ，用x ( o ) 计算得到p 0 以及h ( 1 ) ， 2 ) 当 z 。m 。时 3 )求解线性特征值问题爿【) x = x a 的最小南个特征对a 和x ( o ， 4 ) 计算p 沈= d i a g ( x ( ) x ( 日) ， 5 )计算p 和趔的差。如果差值小于阀值，程序终止。 6 ) p 警1 ) = ( 1 一口) 西：+ 口p 曼。， 7 )用p ”形成日( 件1 ) ， 8 )i = i + 1 ，回到2 ) 。 算法1 2 与算法1 1 的区别在第6 步。 除了简单混合，还可以选择其它方法计算虞”。1 9 8 0 年p u l a y 假设r 纠是 关于p 的线性函数，提出了p u l a y 混合【9 】。以拟n e w t o n 方法为基础j o h n s o n 在 1 9 8 8 年提出了b r o y d e n 混合方法 4 l 。另外还有其它的电荷密度混合技巧可参看 文献【1 l ，1 2 】等。 虽然这些混合方法对s c f 迭代有改善作用，但是对于这类方法一直没有合 理性的说明，而且求解的同时也带来了一些难度，例如简单混合罩面q 的选 取、p u l a y 混合的强假设、b r o y d e n 混合要求解一个泛函极小化等。我们尝试解 释电荷密度混合的合理性，并寻求更高效的混合方法。 一3 一 第一二节密度矩阵的存在性问题 第二节密度矩阵的存在性问题 本节我们从电荷密度混合技巧中引出一个密度矩阵的存在性问题，并通过 构造法给出该问题的一个解。 2 1 密度矩阵的存在性问题 记，是礼阶单位矩阵；t r ( a ) 表示矩阵a 的迹，也就是a 所有对角元的 和；忙表示向量x 的p 范数。 在纯粹的s c f 迭代中，每一步得到的电荷密度近似都对应着一个密度矩阵 的近似。在电荷密度混合技巧结合s c f 的算法中，第l 步生成h a m i l t o n 之前， 是把当前得到的电荷密度近似p 勉与前一步或者前几步豹电荷密度作混合，混 合后得到电荷密度近似p 5 ，并保证l l 醛1 f i l = 后。我们关注：对于经过混合 得到的电荷密度近似，是否存在一个密度矩阵的近似戈与之对应? 这就归结为 以下问题： 问题2 1 ：设p = ( p l , p 2 ，肪) t ，0 胁1 ，并且i i p l l l = k ，k 和n 都是自然 数并且后 p 2 ，p 2 十p 4 1 ( 2 3 ) a 。：! ! 卫匕旦丛l 二出 血一p l 砒：! 型坐兰州 p 2 一p 1 - 6 第一二节密度矩阵的存在性问题 胁：丝丝型! ! 型 p 2 一p l 履：尘型! 趔! 型 融一p 1 z = ( 沥，0 ，一俩，同t ，y = ( 0 ，沥，俩，佣t 。 如果( 2 3 ) 式成立，那么 1 ) q 3 + 风= p s ，o r 4 + 屈= p 4 ，0 3 屈= o t 4 尻。 2 ) 0 a 3 ，o r 4 ，风，屈 1 。 3 ) 护可= 0 ，i 2 = j l y 1 2 = 1 。 证明：1 ) 根据d 3 、岛、a 4 、z 4 的定义代入即得。 2 ) 由( 2 3 ) 式以及1 ) 易知0 0 f 3 、风、0 1 4 、屈 4 ，n = 1 ，2 ，竹 ， 砚) 讵是实数，满足佻= k ，其 中2 k n ，0 他 i ( i = 1 ，2 ，n ) ，并且这仃个实数满足假设2 4 。 一7 一 第二苇密度矩阵戆夺在 生勰撩 那么 讹 龌n 可以分成四缰，弛为g 、q 、岛、q ，满足以下条件 1 ) g 3 、q 巾备只有一个数： 2 ) 记s ；为q o = 1 ，2 ，3 ，4 ) 中所商数的和。 0 曼8 1 奄一1 ，0 s 2 ，魏，8 4 惫一i ，8 2 - t - s 4 1 g 4 ) i t , l l ：设y ，满足 叭1 。z ；m 1 = r a i n 一蕃嘲，；e y 一 。 避p 都么n x y 枣至少存秀个数。黪不然，无论n y 审骞一个数或者n y 跫空 集，由y 的宛义以及他 k 一1 。 令 馈= x l ，岛= 瓢) ，岛= ，岛一 姚k y g ， 记蔻& g l ，2 ，3 ，鸯审麝有数瓣移那么s l 蠡一1 ，孬茭| l 藏毒忿+ 善魄 l ，这与y 的定义矛盾。易见0 鬟s 2 ，s 3 ，s 4 1 ，并殿肖 s 1 + 8 4 = z l + 啻k 一1 ， s 2 + s 4 = m i l ， 德， 这样，我们裁得捌了一个努缓g g = l ，2 ，毒4 ) ，涛怒1 ) 2 ) 。 ( 2 ) 如果( 1 ) 不成立，即对任意的9 豫) 培y ，x l + y k 一1 ，特别的有 茹1 + 玑 奄一1 第l + y t - - i 戋n 搿2 轨。 此时，所有数排列为 挈1 伽篓蜘曼犰霉2 z l 令岛= 。2 ，c 4 = u d 。设r 满筵 勋+ y i k 一1 令c 1 = 轧t ，势 u x l ，一 y l ，汝，势 ，郑么岛擎空，蚕粼r 一0 ， e 送雨靛育搿l 十嚣 l ，矛露e 这榉我常j 就存 缸= 1 t 8 1 + s 4 = z l + 驰 k - 1 ，8 2 + 8 4 他 1 t = r i e y 这样我 | 、j 载我搿了满是1 ) 耸筋分缝。 综上，弓f 理得证。一 在引理2 8 中，令奄= 2 ，对予绦定的p = ( p t ，以，p o ) t ，例f l = 2 ，利用引 理2 8 豹避号，我囊把豹珏令分曩分成1 ) 2 ) 豹q g = 1 ，2 ，3 ，萄。 餐焉这个分缀形式我 f 】绘凄矽，弼2 ) 戆一今勰。 命题2 9 ：使用上面的记号，设a = c p 叁i ，岛一 毋) ；n 句- - 2 + 。以下两个问题 的解等价： ( 1 ) 求铴，m ，热，热，覆得 ) 0 嘞，o - 4 ，热，蕊。 i i ) q 3 + 鹰一s 3 ，锄+ 反= 8 4 ，蝴恁= 锄蕊。 i i i ) x t s f = 0 ，1 2 = i l y 1 2 = i ， 这里茹= ( 佤，0 ，一俪，同t ，y 。( 0 ，佤瓶，俑t 。 ( 2 ) 求口3 ，a 4 ，蕊，厦，使得 ) 0 2 以及给定的p = ，p 2 ，肪) t ， i p l l , = 七，利用 引理2 8 的记号，我们把p 的n 个分量分成g ( i = 1 ，2 ，3 ，4 ) ，并且满足引理2 8 的1 ) 2 ) 。 命题2 1 0 ：使用上面的记号，设g = c # 名1 ，岛= 孝) n 诂- 件21 。以下两个问题 的解等价： ( 1 ) 求n 3 ，m ，岛，屈，使得 i ) 0 0 ( 3 ，0 ( 4 ，风，伪。 i i ) + 岛= s 3 ，+ 尻= 8 4 ，a 3 岛= 劬厦。 i i i ) z t 耖= 0 ，i i x l l 2 = 七一1 ，i l y l l 2 = 1 ， 这罩z = ( 何，o ，一雁，同t ，y = ( o ，厄，瓶，倜t 。 ( 2 ) 求0 3 ，a 4 ，风，俄，使得 i ) 0 0 3 ，0 1 4 ，风，角。 一1 0 第_ 二节密度矩阵的存在性问题 i i ) 0 3 + 侥= 8 3 ，c 4 十风= 8 4 ，a 3 风= m 风。 = 1 ， 0 ，0 ，0 ，一 磊，同t ， 一，佰面，瓶，厕t ，训e r 。 j，、n z 证明：由8 1 = c 以及s 2 = 可知，以上两个问题都在求解 = 1 t = ，+ 1 口3 + 尻= 8 3 a 4 + 屈= 8 4 a 3 + i f 4 = k 一1 8 1 尻+ 尻= 1 一s 2 口3 风= a 4 尻 0 2 的情况类似。特别地，当k = 1 时，设x = z ，对 应的特征值是a ；当k = 2 时，设x = ( z ，可) ，对应的特征值分别是a l 、a 2 ，这 里z ，y 是n 维向量a 当k = 1 时，显然 1 8 i 弋跳 2l 蕊丽 当k = 2 时，把x 一0 ，y ) 按列分开考虑，我们有 卜( ( 鬈) t ( 苗) t ) 其中鬈以及篱都是n 阶矩阵。类似地 如、 咿( | | j 其中器和雾是n 阶矩阵。从而 ，= ( v ，) t v 9 = ( ( 筹) t ( 苗) t ) 一1 3 一 如、 噩1 ， 8 0l 第三节n e w t o n 型混合技巧 = ( 鬈) t 骞+ ( 筹) t 历o y 现在我们考虑j 的具体形式。由于，( x ) = d i a g ( x x + ) ，针对j 的第一部分 袅。岛= 1 时 笪：odiajg(：xx*)：2出ag(z)ox2 刁f 一2 砌a g _ ， k = 2 时， ( v ，) t = 2 ( d i a g ( 茹) d i a g ( y ) ) 引理3 1 ：设带参数q 的对称矩阵4 ( 理) 的特征值与相应的特征向量 分别为a ) ，z ) ，如果在n o 的某领域中a ) 、a ) 、z ) 关于a 可 导，( 伽) 、( 咖) 、：f ( a o ) 存在，并且a 重数为1 ，x * x = 1 在咖的该领域 中恒成立。那么在口= n o 处，有 一= 一( a a ，) 4 z 其中“4 ”的定义如下：对于n 阶矩阵b ，如果存在n 阶可逆矩阵p ，使得 尸- 1 b p = ( ：吕) ， 这里c 是可逆矩阵，那么 、f0 剖一尸io：。) 一 该引理的细节以及证明可以参看文献 7 】。 根据特征值单重的假设以及特征值问题的形式，上述引理条件满足。具体 地，当k = 1 时， 筹一( n ) 咿为 其中 z ：掣z ：( s d i a g ( z ) ) t 口口 一1 4 一 第三节n e w t o n 型混台技巧 。= 2 时v 9 ：( 奏) ：( 二 ：二a i i ) 4 h x ) 3 3n e w t o n 型混合算法 综上，我们可以建立n e w t o n 型混合的s c f 算法。 算法3 3 ：n e w t o n 型混合s c f 迭代 给定三，最大迭代步数 。；，阀值s ，x ( o ) ， 初始化i = 1 ，用x ( o ) 计算得到矗0 以及日( 1 ) ， 当i 0 是参数，l 和s 如前定义。我们应用n e w t o n 型混合的s c f 算法到 这类简化模型，并跟无混合和简单混合的s c f 算法作比较。 下列命题在l 和卢的一个比较强的假设下得到特征值单重的条件。 命题3 4 ：假设l 的特征值为九，i = 1 ，2 ，n ，从小到大排列为0 a l a 2 k a k + 1 k 。卢是n 维非负向量，满足0 西1 ，z = 1 ，2 ，n ，并且e t 声= k ，e 是元素全为l 的扎维向量。如果 m i n ( a + l a t ) 0 2 i a ； 从而对任意的指标j o ，1 ，后 ， a i + x 一= 知+ 1 + 岛十1 一一乃。啦n ( a 件l 一九) 一2 n 0 这表明，矩阵l + 3 d i a g ( s p ) 的最小七+ 1 个特征值都是币的并且两两不同。即 该矩阵的最小后个特征根都是单重的。 注35 ：命题3 4 只是从理论上保证，工的最小七个特征值没有重根，并且卢适 当小的时候，h a m i l t o n 矩阵最小尼个特征值是单重的假设满足。由于对于实际 的l ，p 和p ，日( 几乎总没有重特征值，在后面的数值试验中，我们没有限制 口的选取。 3 4 2 试验结果 在所有数值试验中，我们取尼= 1 0 ，初始值x ( o ) 随机生成。三维l a p l a c e 算子的离散网格划分为1 2x1 3 1 4 ，即矩阵l 的阶为 = 2 1 8 4 。在简单混合的 s c f 算法中，我们始终取o l = 0 1 。 定义残差e ( 0 = 0 以。一崩肌在三种算法中，我们采用相同的收敛准则 g ( ) 2 的情况类似。在以下 讨论中，当南= 1 时，设x = z ，对应的特征值是入，即b = a ；当= 2 时，设x = ( 而可) ，这罩z ，y 是佗维向量，z ，y 对应的特征值分别是a l 、a 2 ，即 a d = ( a l ，a 2 ) t 。 当k = 1 时。( 4 1 ) 中的f c “，g c ，从而t c ( “+ 1 ) “，同时 z c 。1 ) “，我们对( 42 ) 进行n e w t o n 迭代，此时的j a c o b i 矩阵j 的形式为 j ：f 磊o f 磊, g f1 瓦两 接下来我们给出，中的各项偏导数。 ( 1 ) 瓦o f c ”， a f o l x o ( d i a g ( s p ) z ) a a z 瓦2 瓦十瓦一一百 ( 4 3 ) 其中 瓦。l ， 坐攀型：d i a g ( s p ) + 2 s 。( x x t ) ， 挲：地 这里“o 表示h a d a m a r d 乘积，j 是扎阶单位阵。 ( 2 ) 丽o f “， ( 3 ) 筹c 1 “， o f引一a z ) 两2 督2 哪 磐一o ( z * x - 1 ) 2 ， 毡o x 一 一2 0 第四节直接n e w t o n 算法的探索 ( 4 ) 筹c ， 育o g ：o a a 一 这样，我们就得到了n e w t o n 迭代的j a c o b i 矩阵( 4 3 ) ，并且j c 扣+ 1 ) 。( 件1 ) ， 此时第i 步的迭代格式为 ，( ( z ( 件1 ) 一z ( ) = 一t ( z ( ) ，“5 ) 当后= 2 时。我们先计算各个偏导数的形式。记v 取表示f 关于x 的 j a z o b i 矩阵。 0 ) v e x ， 仿照前一节的做法，把x = ( z ，y ) 按列分开考虑，利用k = 1 已有的结果，就 v c l x ，x = ( 。三。：”) ， 同理我们可以得到 v c 出a g c s 力x k = ( 出a g 。( s 。p 出芝葛力) + z ( ；) 。( 荔：) ， ， 可c 酬x = ( 之譬) ， 这里u = 1 ，2 ) 是一个对角元全为的扎阶的对角矩阵。 ( 2 ) v f a 。， 一2 l 一 l 可 。 ，i一， 一 i l d 取 v 魄 可 p 第四节直接n e w t o n 算法的探索 ( 4 ) v g d ， v g x = 2 铲 0 1 x n) v g d = 0 2 2 令 扛( v v 败e xv f o v g 。) ， 力 v 败d 。 。 按照上面各个偏导数的形式，可得d c ( 槲2 ) ( 轨十2 ) 。 把z ( 件1 ) 一z ( 旬c ( 叶1 ) 2 的列进行重排，记 + 1 ) - z ：f ，如幻1 6 a lj a 2 其中如，6 y c ”“，以l ，6 a 2 c 。设 妒k 睁) l ， ， 同理，把t ( z ( o ) c ( ”+ 1 ) 。2 的列仿照( 4 8 ) 进行重排，设重排后的结果为 2 ( z ( 0 1 c ( 轨+ 2 ) “， ，由( 4 7 ) 定义，那么此时的n e w t o n 迭代格式为 4 3n e w t o n 迭代法 ) y = 于( z ( ) ) ( 4 9 ) 根据上面的讨论，我们可以给出求解( 4 2 ) 的n e w t o n 迭代法。由于当k 1 的时候，算法的应用涉及一些另外的处理方法，我们在此不做详细讨论，这罩 仅给出当k = 1 时算法的形式。 算法4 1 ：求解非线性方程组( 4 2 ) 的n e w t o n 迭代法( 七= 1 ) 一2 2 第四节直接n e w t o n 算法的探索 o ) 初始化z ( ，a ( ， 1 ) f o r i = 1 ，2 ，一， 2 )计算( z ( ) ) ， 3 )检查收敛性， 4 )计算j a c o b i 矩阵，( ”， 5 )从j ( ( z ( 件1 ) 一z ( ) = 一t ( z ( ) 中解出z ( 件1 1 ，a ( 件， 6 1e n d f o r 对于算法4 1 ，我们有 命题4 2 ：如果日只含线性项，那么算法4 1 就相当于反幂法。并且当日对称正 定时，算法41 收敛。 证明：当南= 1 时，由于日跟x 无关，那么 箬= 只篆= 喵筹= 矿，筹= 。， 此时的j a c o b i 矩阵 j = ( 二 在第i 步，n e w t o n 迭代格式为 了) ( 篡t x 。( 0 ) ( 篇二= 一( 筑并) ， 比较等式两边的第一个元素可得 z ( 件1 ) ：a “十1 ) h 一1 z ( ) 更新z ( 件1 ) 的方向为h x ( 0 ，也就是说，此时的算法4 1 相当于反幂法。 接下来我们说明，当日对称正定的时候，了可逆，即算法4 1 收敛。 设t ，c n “，a c ，使得 j ( ：) = ( 暴- - 。x ) ( ：) = 。 如果口= 0 ，那么根据上式可得h v = 0 ，由h 对称萨定可得t ，= o ；如果 一2 3 第四节直接n e w t o n 算法的探索 a 0 ，根据上式可得口= 0 ，进而有 h 口一。z = 0 = u = a h 一1 z 号z t 甜= a 2 t h 一1 z = 0 由于z 0 ，上式与日是对称正定矛盾。 4 4 数值例子 我们把算法4 1 应用到以下简化物理模型中( k = 1 ) h ( p ) x = 尬 ，z 一1 = 0 p = d i a g ( x x + 、 其中 h ( p ) = l + ( f i a g ( s p ) 工是三维离散l a p l a c e 算子，对两个不同的n 做实验，网格划分为1 0 1 0 1 0 以及1 2 1 3 1 4 ，即礼= 1 0 0 0 以及竹= 2 1 8 4 。在算法4 1 的第5 步，我们用 直接法求解线性方程组。初始向量由随机生成。定义残量e r r ( ) = i i h ( p ( o ) z ( 日一 a ( ) z ( ；) 忆当e r r ( o 1 0 1 0 时算法终止。 对于绝大部分初始值，算法都是收敛的。图4 1 是其中两个例子，由收敛曲 线可见，算法达到了二次收敛。 ( a ) n = 1 0 0 0 ，七= 1 i 璺| 4 ，1 收敛曲线 选十t 步量 ( b ) n = 2 1 8 4 ，k = 1 注4 3 ：虽然对于绝大部分随机生成的初始值，从结果上看该算法都能达到二次 收敛，但是对于包含非线性项的h a m i l t o n ，算法4 1 中j 的有可能在迭代过程 一2 4 第四节直接n e w t o n 算法的探索 中奇异，并且当，奇异的时候，算法未收敛到真实解。在实际应用中，我们可 以通过一定的技术处理使得j 在迭代过程中保持非奇异，目前该算法的细节还 在研究中。 我们可以构造一个初始值，使得j 是奇异的，但是该初始值并不是非线性特征 值问题的解。 例4 4 ：设竹= 2 ，k = 1 ，z = ( x l ，x 2 ) t r 2 ， 三= ( 二1 1 ) 一= ；( ；：) 考虑非线性特征值问题 h ( x ) = l + 2 d i a g ( s p ) ， h ( x ) x 一妇= 0 ， ( 4 1 0 ) ，o = 1 酢，x - a x = 舭萎鬻裂五= 。 j a c o b i 矩阵j 为 ；r 之蓼- 3 + 8 砰 4 x l x 2 3 - - 3 x l 6 十2 x + 4 x l 一3 a + 8 x ；一3 2 2 2 x 2 0 容易得到，特征对( 牙，譬) 使得j 矩阵奇异。但是，( 曩警) 并不是原非线性特征 值问题( 4 1 0 ) 的解，而( 牙，2 ) 才是原问题的特征对。这说明，从理论上看，使得 j 奇异的近似值有可能不是原问题的解。 4 5 小结 本节我们把非线性特征值问题的转化看做是非线性方程组，利用向量值函 数微分的相关理论，给出了n e w t o n 迭代法的j a c o b i 矩阵计算方法以及迭代格 式，并建立了k = 1 情形的n e w t o n 迭代法。通过简单的数值实验，从结果上 看，算法达到了二次收敛。 一2 5 一一 叁耋塞坚 参考文献 【1 1a l a i l d vk w c h ua n d pl a n c a s t e r , d e r i v a t i v e so f e i g e n v a l u e s a n d e i g e n v c o t o r s o f m a t r i x 铀锄o n s 锄j 0 t l m a l 可m n 蚋x 舢n 岫s 矾d a p p l i t t o n 1 4 ：9 0 3 9 2 6 , 1 9 9 3 1 2 1z b a i ，j d e m m l , j d o n g a r r a , a r u h e ，a n dh v a nd e r v b m , e d i t o r s t e m p l a t e s f o r s o l u t i o no f a l g e b r a i c e i g e n v a l u e p r o b l e m s ：a p r a t i c a l g u i d e s i a m , p h i l a 蛐h i a , 2 0 0 0 , 3 1g g o l u b a n d c v a nl o a n m a t r i x c o m p u t a t i o n t h ej o h n s h o p k i n s u n i v e r s i t y p r e s s , b a l t i m o r e , t h i r de d i t i o n , 1 9 9 6 4 1d d j o h n s o n m o d i f i e db r o y d e n sm e t h o df o ra c c e l e r a t i n gc o n v e r g e n c ei ns c l f - c o u s i s m n tc a l c u 1 a t i o n s ，p h y s i c a l r e v i e w b 3 8 ：1 2 0 8 7 1 2 8 1 3 ，1 9 8 8 【5 】g pk e r k e r e f f i c i e n ti t e r a t i o ns c h e m ef o rs e l f - c o u s i s t o mp s e u d o p o t e n t i a lc a l c u l a t i o n s e h y 耐c a l r e v i e w 最2 3 ：3 0 8 2 3 0 8 4 1 9 8 1 【6 】gk r e s s ea n dj f m t t l m n l l e r e f f i c e m yo f 曲i n i f i o t o t a l e n e r g yc a l c u l a t i o n s f o r m e t a l s a n d 蝴n i c o n d u c t o r s u s i n g a p l a n e - w a v e b a s i ss e t c o m p u t a t i o n a l m a t e r i a l s s c i e n c e 。6 ：1 5 - 5 0 1 9 9 6 川c d m e y e ra n dg ws t e w a r t d e r i v a t i v e sa n dp e r t u r b a t i o n so fe i g e n v e c t o r s s i a m j o u r n a lo f n u m e r i c a l a n a l y s i s , 2 5 ( 3 ) ：6 7 9 - 6 9 1 ，1 9 8 8 i s m c p a y n e ，m pt e t e r , d c a l i 越ta a r i a sa n dj d j o a n n o p o u l o s i t e r a t i v em i n i m i z a t i o n t e c h n i q u e sf o ra bi n i t i ot o t a l - e n e r g yc a l c u l a t i o n s ：m o l e c u l a rd y n a m i c sa n dc o n j u g a t eg r a d i e n t s r e v i e w so f m o d e r np h y s i c s 6 4 ( 4 ) ：1 0 4 5 1 0 9 7 ，1