一种有效的广义特征值分析方法

一种有效的广义特征值分析方法 Ξ 一种有效的广义特征值分析方法 刘寒冰 龚国庆 刘建设 ( )吉林大学南岭校区交通学院 ,长春 ,130025 摘 要 提出了一种适合于自适应有限元分析中求解广义特征值问题的多重网格方法. 这 种方法充分利用了初始网格下的结果 ,通过插值或最小二乘拟合技术来得到网格变化后的新的 近似特征向量 ,然后由多重网格迭代过程实现对结构广义特征值问题的求解 . 在多重网格迭代 的光滑步中 ,选择了收敛梯度法以提高其收敛率 ;在粗网格校正步中 ,则导出了一种近似求解特 征向量误差的方程 . 这种方法将网格离散过程和数值求解过程很好地相结合 ,建立了一个网格 细分后广义特征值问题的快速重分析方法 ,与传统有限元方法相比较 ,具有计算简便 、计算量少 等特点 ,可以作为结构动力问题自适应有限元分析的一种十分有效的工具 . 0 引言 关键词 多重网格迭代 ,广义特征值问题 ,自适应有限元 用自适应有限元方法分析结构动力问题时 ,有限元网格经常会发生变化 ,若用传统的有 () 限元方法 如子空间迭代法去分析自适应过程则存在一定的缺陷 ,因为传统的有限元方法 () 包含的两个过程 即离散过程和数值求解过程是相互独立的 ,两者无相互作用 ,这就造成了 () () 很大的浪费 : 1没有充分利用初始网格下的结果 ; 2当网格逐步加密时 ,传统有限元方法 的计算量会很快增大 ,当网格很多 ,自由度数很大时 ,计算量会更大 ,但收敛效果不一定好 , 这种缺陷对于 动 力 问 题 会 更 明 显. 尽 管 目 前 对 于 特 征 值 问 题 发 展 了 一 些 有 效 的 算 法 , 如 1 ,2 3 Lanczos 方法、Arnoldi 方法等 ,以上的缺陷仍然存在 . 多重网格方法则弥补了传统有限元方法之不足 ,它将离散过程和数值求解过程相结合 , 利用初始网格下的结果 ,通过多重网格迭代 ,可以达到简化计算过程 ,节省计算量的目的. 目 前 ,多重网格技术已经发展到了可以应用于实际工程问题 ,商用软件也已出现 ,线性和非线 性多重网格求解器建立了高效的求解偏微分方程的方法 . 用多重网格方法求解特征值问题 也引起了广泛的关注 . 下面回顾一下多重网格技术在特征值问题方面的研究状况 . 传统的多重网格方法求解特征值问题是将此问题转换成非线性问题 ,然后利用非线性 4 () 多重网格求解器 ,如“全近似方法”full approximation scheme进行求解. 这种求解器应用于 非 5 ,6 线性问题通常是有效的,但对于线性特征值问题 ,这种技术未必总是有效的 ,因为它没 有充分利用特征值问题的一些特性 ,而且 ,非线性多重网格特征求解器通常要对特征簇进行 特殊的处理 ,这又进一步使程序复杂化了. 此外 ,此技术的另一个缺陷就是对粗网格解的近 似程度要求较高. Ξ () 国家自然科学基金 19872029资助. 2002203205 收到第 1 稿 ,2002205217 收到修改稿. ( ) 第二种比较流行的方法是采用一个外部特征值求解器 ,如 Rayleigh 迭代 RQI,它要求 7 解一个移位的刚度矩阵构成的线性方程组,而多重网格迭代则作为一个内部迭代求解器 , 从而构成一个内2外迭代过程 . 使用移位的刚度矩阵的目的主要是使外部迭代过程具有快的 8 收敛性 ,但要注意的是 ,由多重网格迭代求解的线性方程组有时几乎是奇异的. 与此类似 () 的技术如经典的逆迭代 通常是逆子空间迭代作为外部求解器 ,而多重网格迭代仍作为内 9 部迭代求解器. 这种逆迭代技术的收敛性要比上面提到的移位的迭代法的收敛性要慢 ,但 不会遇到求解近乎奇异的线性方程组的麻烦 . 10 ,12 第三种目前有些学者研究得较多的一类方法叫做“多重网格预条件梯度法”.“预条件是一种改进矩阵条件数的技”术. 这类方法一个最明显的特征是 ,它不包含内2外迭代两 () 个过程 ,多重网格迭代仅作为一个“预条件器”preconditioner,它既不是一个完整的线性求 解器 ,也不是一个特征求解器 . 这类方法中较有前景的一种方法叫做“局部优化预条件收敛 13 ,14 15 梯度法.” 与之类似的方法还有 ,预条件最速下降法 、预条件逆迭代法、预条件子空间迭 代法,等等. 本文将提出一种基于插值或最小二乘拟合技术的多重网格方法 ,并通过收敛梯度法来 改善多重网格迭代的收敛性 ,目的在于建立一个网格细分后的快速重分析方法 ,即从初始解 得到重新划分网格后的结构特征值的简便算法 . 通过数值算例与传统的广义特征值求解方 () 法 如子空间迭代法所得到的结果进行了比较 ,说明了该方法的有效性和实用性 . 1 理论推导 在进行有限元动力分析时 ,首先通过有限元方法建立结构粗网格下运动微分方程 ()[ M ]{ u? } + [ C ]{ u } + [ K]{ u} = { F} 1 ( ) s , t 这里 , [ M ]是系统质量矩阵 , [ C ]是粘性阻尼矩阵 , [ K]是刚度矩阵 ,{ F}是外力矢量 ,{ u} ( s , t) 是节点位移矢量 ,{ u }和{ u? }分别 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示节点速度和加速度. 对于无阻尼的固有振动 ,特征值问题 ( λ) () 2 [ K] - [ M]{ u} = {0}l - 1 l - 1 l - 1 l - 1 λ式中 , [ K]和 [ M]分别表示初始网格下的刚度阵和质量阵 ,和 [ u] 表示相应的特 l - 1 l - 1 l - 1 l - 1 征值和特征向量. 有许多成熟的方法可以求出初始网格下的特征频率和特征向量. 1. 1 作为延拓的线性插值和最小二乘拟合 在求得初始网格下的特征向量{ u} 以后 ,通过线性插值 l - 1 x y x 可以得到细网格下的特征向量{ u} . 如图 1 描述了最常用的初 l ( ) ( ΩΩ网格和细网格 ,它们分别是网步为 h, h和 h, l - 1 l l - 1 l - 1 l y x x y y ) h= 2 h, h= 2 h. 为表达方便 ,以下假 h的矩形网格 ,且 l - 1 l l - 1 l l x y x y ΩΩ设 h= h= 2 h , h= h= h ,即和是正方形网格 ,粗 l - 1 l l - 1 l - 1 l l () 网格的步长是细网格的两倍. 对粗网格进行细化 ,设点 a 0 ,0, ΩΩ图 1 初网格 和细网格 l - 1 l() (() ) ΩΩΩΩb 0 ,2 h, c 2 h ,0, d 2 h ,2 h??,这里和分别 表l - 1 l l - 1 l 示步长为 2 h , h 的正方形初网格和细网格 . 对于已知的{ u} ,定义 l - 1 第 4 期 刘寒冰等 : 一种有效的广义特征值分析方法 ?421 ? (){ u} = [ P ]{ u} 3 l l - 1 式中 { u} = { uu u} uuuuu l a1u d b 2345c { u} = { u u} uu l - 1 ad bc T Π2 1Π2 1Π4 11 0 0 0 0 0 Π2 10 1Π2 1 0 1Π4 0 0 0 P ] = [ 1Π4 0 0 0 10 Π2 1 1Π2 0 1Π2 0 1Π4 1Π2 0 0 0 0 1 矩阵 [ P ]称为粗网格向量的延拓 . 除了通过线性插值可以由粗网格特征向量{ u} 得到细网格特征向量{ u} 外 ,也可以 l - 1 l 通过最小二乘拟合技术来得到细网格的特征向量. ( ) ( ) 设已知矩形区域内 N ×M 个粗网格点 x, y, n = 1 ,2 , ?, N ; m = 1 ,2 , ?, M 上的函 n m l - 1 数值 u ,则由最小二乘拟合多项式 ( ) x , y p q i - 1 j - 1 ( ) = ax y4 u ( )ij x , y?? i = 1 j = 1 ( ) 可以求得矩形区域内任意一点处的函数值 u . 式中 ai = 1 ,2 , ?, p ; j = 1 ,2 , ?, q为待 ( ) x , yij 定系数 ,它可以根据最小二乘原理来得到 . p , q 分别表示拟合多项式中 x , y 的最高次数加 1 ,即 x , y 的最高次数分别为 p - 1 和 q - 1. 由插值或拟合得到的细网格下的特征向量{ u} ,再根据瑞利商 l T { u} [ K]{ u} l l l ()λ5 A= l T { u} [ M]{ u} l l l λA可得细网格下的近似特征值 . l 1. 2 光滑过程 对于细网格下的广义特征值问题 ()6 λ[ K]{ u} = [ M]{ u} l l l l l () () λA若以 4式和 5式得到的近似特征向量 和特征值{ u}作为初始的迭代特征对 ,经过迭代 l l () λA过程 ,可以得到细网格下的特征方程 6的解{ u} . 事实上 ,由于 和{ u} 是细网格下较好 l l l 的近似解 ,因此 ,在实际计算中经过少数几次迭代就可以了 . 对于迭代过程 ,既可以采用传统 的 Gauss2Seidel 迭代或雅可比迭代 ,也可以采用其它迭代方法 ,这里采用的是收敛梯度法 . 由 于 Gauss2Seidel 迭代或雅可比迭代的收敛性主要由迭代矩阵的谱半径来 决定 郑伟家庭教育讲座全集个人独资股东决定成立安全领导小组关于成立临时党支部关于注销分公司决定 ,有时即便 [ K] l 是对称正定的 , Gauss2Seidel 迭代或雅可比迭代也不总是收敛的 ,而收敛梯度法是求解稀疏 线性方程组 ,特别是对称正定方程组的最有效的方法之一 ,而且 ,其矩阵与向量相乘的次数 也是很少的 ,所以 ,选用了此方法 . () 若令方程 6的右端项为 λ(){ b} = A[ M]{ u} 7 l l l () 则求解方程 6的收敛梯度法的迭代过程为 算法 1 . 1 CGM (λ) 输入 : m 个起始迭代特征对 A, u, j = 1 ,2 , ?, m , j 表示特征对的阶次 . j j λA计算 : b=[ M]{ u} , j = 1 ,2 , ?, m j l l j j () ()00 () 1令 : d = r= b jj j () 2迭代 :对 i = 0 ,1 , ?,直到收敛 T ( ) ( )ii ( )i r r () α3= Tj ( )( )i i d [ K] d l j j( )( )( ) ( )i + 1 i ii () α= u 4u + d jj jj)( )( )( ( )i i i i + 1 () α5r = r- [ K] d l jj j j T ( ) ( )i + 1i + 1 ( ) r r i + 1 j j () β6=T j ( ) ( )ii rr j j( ) ( ) ( ) ( )i + 1i + 1i + 1i () β7d = r+d j j j j ()8 结束 () λ输出 :新的特征向量{ u} ,再根据瑞利商 5式得到新的特征值 , j = 1 ,2 , ?, m j j 此迭代过程可以简写为 ( )( )i i ( )i ()8 { u } =C{ u } (λ) A jj j ( )i 上述迭代中 ,有矩阵和向量相乘的只有第 3 步和第 5 步的[ K] d ,所以计算量很少. l j 1. 3 粗网格校正过程 对于粗网格校正 ,需要计算残量 ()9 λ( ) { d} = [ K] - [ M]{ u} l l l l l 通过对细网格的限制 ()10 { d} = [ r ]{ d} l - 1 l 可以得到粗网格下的残量{ d} . 限制最简单的取法是平凡单射 r ( ) ()ΩΩrl - d1 inj11 = d, Π x ?< )( inj l l - 1 xl - 1( )lx 尽管这种限制容易作出 ,但它还是有一些缺点 ,为了使限制取得更安全些可参照文献 [ 16 ] . 为了得到粗网格的特征向量误差 ,设 3 3 Δ Δλλλ()= A- 12 {u} = { u} - { u } ,l - 1 l - 1 l - 1l - 1 l - 1 l - 1 3 3 λλ,{ u }分别表示初始网格下精确的特征值和特征向量 ;A ,{ u 式中 }分别表示初始 l - 1 l - 1 l - 1 l - 1 ΔλΔ 网格下近似的特征值和特征向量 ;,{u}分别表示特征值和特征向量的误差 . l - 1 l - 1 则残量 λ( ) { d} = [ K] - A[ M]{ u} = l - 1 l - 1 l - 1 l - 1 l - 1 3 3 (λΔλ) Δ ( ) () K] - [ + [ M]{ u } + {u}=l - 1 l - 1l - 1 l - 1 l - 1 l - 1 3 3 3 ( λΔλ) [ K] - [ M]{ u } - [ M]{ u } +l - 1 l - 1 l - 1 l - 1 l - 1 l - 1 l - 1 3 ( λΔΔλΔ) () 13a [ K] - [ M]{u} - [ M]{u}l - 1 l - 1 l - 1 l - 1 l - 1 l - 1 l - 1 3 3 ( λ) 略去高阶小量 ,并由 [ K- [ M]{ u } = {0} ,可推得 l - 1 l - 1 l - 1 l - 1 3 3 3 λ) ( λ) Δ (λ ()A { d} ? - - [ M]{ u } + [ K] - [ M]{u}13b l - 1l - 1 l - 1 l - 1 l - 1 l - 1 l - 1 l - 1 l - 1 3 3 λΔλ令 =A 和由{ u } = [ r ]{ u } ,可得粗网格下特征向量误差{u } l - 1 l l - 1 l l - 1 第 4 期 刘寒冰等 : 一种有效的广义特征值分析方法 ?423 ? ( λΔ (λλ) ) ()[ K] - A[ M]{u} = A- A[ M] [ r ]{ u} + { d} 14 l - 1 l l - 1 l - 1 l - 1l l - 1 l l - 1 () 14式可简写为 Δ (){u} =15 Q{ d} l - 1 l - 1 l - 1 1. 4 误差的延拓 () 采用和 3式相似的形式 ,通过对粗网格插值或拟合处理 ()16 Δ Δ {u} = [ P ]{u} l l - 1 Δ 可以得到细网格下的特征向量误差{u} . l 1. 5 精确校正( 细网格校正) 3 精确校正为 ()Δ 17 { u } = { u } - {u } l l l 3 T 3 { u }[ K]{ u } l l l 33 λλ最后由瑞利商 =得到细网格下的特征值 . l l 3 T 3 { u } [ M ]{ u } l l l λ( 以上介绍的过程是从 l - 1 层到 l 层的二重网格迭代 ,由于 l 层的特征方程 [ K- l l ( λ) ) [ M]{ u} = {0}与 l - 1 层的特征方程 [ K- [ M]{ u} = {0}具有相同的形式 , l l l - 1 l - 1 l - 1 l - 1 因此可以递推地定义算法 . 在 0 层的方程必须是精确求解 ,但 0 层对应于最粗的网格 ,因此 方程组的个数最少 ,计算量小 . 这样所得的 l 层迭代称为多重网格迭代 . 结构固有特征问题的多重网格方法为 : ΩΩΩ设整个自适应过程中产生了一系列相互嵌套的网格 ,, ?,,相应的总刚系列为 1 2 l K, K, ?, K,总质量系列为 M, M, ?, M,若要求解第 l 层特征方程 ,即 1 2 l 1 2 l ( λ) ()Ω 在上 18 [ K] - [ M]{ u} = {0} ,l l l l l 这时 ,多重网格迭代算法可描述如下 : ()算法 1 . 2 结构广义特征值问题的多重网格方法 MGGEP 0 0 λ, u输入 :初始网格下的特征对 , j = 1 ,2 , ?, m . m 表示所要求的特征对的阶次 j j ()1 迭代 :对 i = 1 ,2 , ?, l ( ) ( ) ( )iT ii { u } [ K ] { u } ( )( )( )i i - 1 i j j ()2 λ{ u } = [ P ]{ u } ; A= )( )( ( ) i i iT j j l { u } [ M ]{ u } j j ( ) ( ) ( )iT ii { u } [ K ] { u } ( )( )( )i i i j j ( )i ()3 λ} ; = { u } = C(λ { u ( ) ( )( ) A) iT i i j jlj { u } [ M ]{ u } j j ( )( )( )( )( )i i i i i ) ( ()λ4 d } = [ K - [ M ]{ u } { j j j )( ( )i i - 1 () 5{ d } = r{ d } inj j j ( ))( i - 1 i ()6 Δ{u } = Q{ u } l - 1 j j ( )( )i i - 1 ΔΔ()7 {u } = [ P ]{u } j j ( ) ( )3 iT ( )3 i i { u } [ K ] { u } )( )( ( )( )i i 3 i 3 i j j ()8 λΔ = { u } = { u } - {u } ; ( ) ( )3 iT ( )3 i jj j j i { u } [ M ]{ u } j j ()9 结束 3 ( i) 3 ( i) (λ) 输出 : l 层细网格下的多重网格特征解 ,{ u }, j = 1 ,2 , ?, m j j 说明 :上述迭代过程中 ,第 3 步和第 6 步均采用了收敛梯度法提高其收敛率 ; 第 2 步和 第 7 步中的延拓 [ P ] 具有相同的形式 ,它既可通过线性插值得到 ,也可由最小二乘拟合得 到 . 由于多重网格迭代的初始特征对是细网格下较好的近似值 ,因此 ,在用收敛梯度法进行 迭代求解线性方程组时 ,只需很少的迭代次数就能得到很好的近似解 ,因此 ,计算量不大 . 2 数值算例 下面以一四边简支的矩形薄板为例 ,说明多重网格方法与传统有限元方法相比较所得 11 2 μ 到的效果 . 其长度 l = 5. 0m ,宽度 b = 2. 5m ,厚度 t = 0. 25m , E = 2. 067 04 ×10NΠm,= 0. 3 , 3 3 γω ω = 7. 833 44 ×10kgΠm. 此 结 构 前 四 阶 固 有 圆 频 率 的 精 确 解 分 别 为 = 767. 108 ; = 1 2 ωω1 2271373 ;= 1 994. 481 ;= 2608. 168 ;单位均是 radΠs. 3 4 ( 表 1 是对图 2 的网格加密后 ,生成图 3 的网格时 ,用多重网格方法求得的细网格下 图 () ω) ω3的固有圆频率 和用子空间迭代法直接求得的细网格下 图 3的固有圆频率 之间的 m s 比较. 表 1 两种不同方法求得的前三阶固有圆频率的比较 1 2 3 i 运行时间Πs ω 736. 899 1 120. 546 1 811. 033 0 ω 757. 241 1 190. 451 1 917. 619 3. 08 s ω 757. 248 1 191. 341 1 967. 255 2. 20 m ω(( ) ) : 粗网格 4 ×2下的固有频率 radΠs,如图 2 . 0 ω(( ) ) : 用子空间迭代法直接求得的细网格 8 ×4下的固有频率 radΠs,如图 3 . s ω(( ) ) : 用多重网格方法求得的细网格 8 ×4下的固有频率 radΠs,如图 3. m ()()图 2 粗网格 4 ×2 图 3 细网格 8 ×4 ( ( ) ) 由表 1 可以看到 :用多重网格方法 ,利用粗网格 图 2下的解所得到的细网格 图 3的 ω() ω解 与用传统有限元方法 子空间迭代法直接得到的细网格下的解 相比较 ,结果十分 m s 接近 ,说明两种不同方法的收敛性是相近的 ,但在计算时间上 ,多重网格方法约少 0. 9 s. (() ) 表 2 是细网格 8 ×4下 如图 3用两种不同方法得到的前三阶固有频率下特征向量的 z 比较. 表中{ u}表示第一阶固有频率下 ,各节点在横向即 z 方向上的幅值的含义类推 . 为了 1 z z 节省篇幅 ,{ u} 、{ u}特征向量中的角度分量均略去 ,只给出横向的幅值 ,同时 ,由于各边界 2 3 节点处的幅值均为零 ,因此 ,在表 2 中也不列出 . 第 4 期 刘寒冰等 : 一种有效的广义特征值分析方法 ?425 ? 表 2 前三阶固有频率下特征向量的比较 子空间迭代法 多重网格方法 节点号 z z zz z z { u}{ u}{ u} { u} { u} { u}1 2 3 1 2 3 7 - 0. 048 7 0. 076 2 - 0. 082 1 - 0. 027 1 0. 070 3 - 0. 052 6 8 - 0. 068 8 0. 107 8 - 0. 116 1 - 0. 038 3 0. 099 4 - 0. 074 4 9 - 0. 048 7 0. 076 2 - 0. 082 1 - 0. 027 0 0. 070 3 - 0. 052 6 12 - 0. 089 8 0. 107 8 - 0. 062 8 - 0. 050 1 0. 099 0 - 0. 038 2 13 - 0. 127 1 0. 152 5 - 0. 088 8 - 0. 070 8 0. 014 0 - 0. 054 0 14 - 0. 089 8 0. 107 8 - 0. 062 8 - 0. 050 1 0. 099 0 - 0. 038 2 17 - 0. 117 5 0. 076 2 0. 034 0 - 0. 065 4 0. 069 7 0. 027 5 18 - 0. 166 1 0. 107 8 0. 048 1 - 0. 092 5 0. 098 6 0. 038 9 19 - 0. 117 5 0. 076 2 0. 034 0 - 0. 065 4 0. 069 7 0. 027 5 22 - 0. 127 1 0. 000 0 0. 088 8 - 0. 070 8 - . 000 0 0. 065 6 23 - 0. 179 8 0. 000 0 0. 125 6 - 0. 100 1 0. 000 0 0. 092 7 24 - 0. 127 1 0. 000 0 0. 088 8 - 0. 070 8 0. 000 0 0. 065 6 27 - 0. 117 5 - 0. 076 2 0. 034 0 - 0. 065 4 - 0. 069 7 0. 030 6 28 - 0. 166 1 - 0. 107 8 0. 048 1 - 0. 092 5 - 0. 098 6 0. 043 3 29 - 0. 117 4 - 0. 076 2 0. 034 0 - 0. 065 4 - 0. 069 7 0. 030 6 32 - 0. 089 8 - 0. 107 8 - 0. 062 8 - 0. 050 1 - 0. 099 0 - 0. 033 6 33 - 0. 127 1 - 0. 152 5 - 0. 088 8 - 0. 070 8 - 0. 014 0 - 0. 047 5 34 - 0. 089 8 - 0. 107 8 - 0. 062 8 - 0. 050 1 - 0. 099 0 - 0. 033 6 37 - 0. 048 7 - 0. 076 2 - 0. 082 1 - 0. 027 1 - 0. 070 3 - 0. 049 1 38 - 0. 068 8 - 0. 107 8 - 0. 116 1 - 0. 038 3 - 0. 099 4 - 0. 069 5 39 - 0. 048 7 - 0. 076 2 - 0. 082 1 - 0. 027 1 - 0. 070 3 - 0. 049 1 (() ) () 表 3 是图 3 的细网格 8 ×4进一步加密后 ,形成更细的网格 16 ×8如图 4,用两种不 (ω) 同的方法计算得到的前四阶固有圆频率的比较 . 表 3 中的 是以网格 4 ×2作为粗网格用 0 (ω() ) 多重网格方法求得的细网格 8 ×4下的固有圆频率 ,如图 3 ;是以细网格 8 ×4作为粗网 格m () ω用多重网格方法求得的更细网格 16 ×8下的固有圆频率 ,而 是用子空间迭代法直接 s () 求得的更细网格 16 ×8下的固有圆频率 ,如图 4 . 表 3 网格进一步加密后 , 前四阶固有圆频率结果比较 运行时间Πs i 1 2 3 4 ω 757. 248 1 191. 341 1 967. 255 2 583. 461 0 ω 766. 309 1 219. 588 1 973. 668 2 599. 876 114. 08 s ω 766. 319 1 219. 664 1 979. 287 2 600. 835 58. 0 m ω(( ) ) : 用多重网格方法求得的细网格 8 ×4下的固有频率 radΠs,如图 3 . 0 ω(( ) ) : 用子空间迭代法直接求得的细网格 16 ×8下的固有频率 radΠs,如图 4. s ω(( ) ) :用多重网格方法求得的细网格 16 ×8下的固有频率 radΠs,如图 4 . m 从表 3 可以看出 ,两种不同的方法所得到的各阶固有频率仍十分接近 ,但多重网格方法的 运算时间要比传统有限元方法少了近一半. ()图 4 细网格 16 ×8 下面比较一下用线性插值和最小二乘拟合所得到的固有频率. 表 4 分别采用插值和拟合两项技术得到的前四阶固有频率结果比较 i 1 2 3 4 第一次加密后 :插值 757. 245 1 190. 615 1 954. 372 757. 248 1 191. 341 1 967. 255 2 583. 461 拟合 766. 315 1 219. 631 1 974. 156 2 596. 263 第二次加密后 :插值 766. 319 1 219. 664 1 979. 287 2 600. 835 拟合 由表 4 可以看到 ,通过最小二乘拟合所得到的结果比用线性插值所得到的结果略好 ,而 且 ,在实际计算中 ,通过最小二乘拟合所得到的初始向量进行的多重网格迭代次数要少 ,因 此 ,用最小二乘拟合来得到细网格的初始向量更理想一些. 需要说明的是 ,本节所举的实例是以均匀网格加密来进行求解的 ,这样做的目的一是为 了编程方便 ,二是容易对结果进行分析和比较 ,找出其中的规律 ,而对于非均匀的网格加密 , 本方法同样是可行的 ,因为通过插值或拟合技术总能求得结构中其它任意点处的近似特征 向量 ,而多重网格迭代过程并不会因此受到影响 ,只不过编程和分析起来稍麻烦一些 . 当然 , 插值过程并不一定是线性插值 ,Lagrange 插值或 Hermite 插值也是可以的 . 3 结论 本文通过插值和最小二乘拟合技术 ,对粗网格的特征向量进行处理 ,以得到细网格的近 似特征向量 ,以此近似向量作为多重网格迭代的初始迭代向量进行迭代来求解网格逐步加 密时的广义特征值问题 . 由于经过预处理后的特征向量是真实特征向量的一个好的近似解 , 因此 ,可以减少多重网格迭代的次数 ,以达到节省计算量的目的 ,同时也能保证解的精度和 收敛性. 在进行多重网格迭代的光滑步时 ,我们抛弃了传统多重网格法中的 Gauss2Seidel 迭 代或雅可比迭代 ,而选择了收敛梯度法 ,以提高求解大型稀疏线性方程组的收敛速度 . 在粗 网格的校正步中 ,导出了一个通过粗网格的残量来近似求解粗网格的特征向量的误差的方 程 . 数值算例的结果表明 ,在结构动力自适应有限元分析过程中 ,随着网格的逐步加密 ,自由 度数越来越多 ,本文所提出的多重网格方法将较为显著地减少运行时间 ,提高计算效率 . 因 第 4 期 刘寒冰等 : 一种有效的广义特征值分析方法 ?427 ? 参 考 文 献 Cullum J , Willoughby R A. Lanczos. Algorithms for Large Symmetric Eigenvalue Computations. Theory , volume 3 1 of Progress in Scientific Computing. Birkhauser , Boston , 1985 Grimes R G , Lewis J G , Simon H D. A shifted block Lanczos algorithm for solving sparse symmetric generalized 2 () eigenproblems. SIAM Journal of Matrix Analysis Application , 1994 , 151: 228,272 3 Sorensen D , Lehoucq R , Vu P , Yang C. ARPACK: an implementation of the Implicitly Restarted Arnoldi iteration that computes some of the eigenvalues and eigenvectors of a large sparse matrix. Available from ftp caam rice edu , directory , pubΠpeopleΠSorensenΠARPACK , 1995 4 Brandt A , McCormick S , Ruge J . Multigrid methods for differential eigenproblems. SIAM Journal of Science Statistics in Computation. 1983 ,4 :244,260 5 Bank R E , Chan T F. A multigrid continuation program for parameterized nonlinear elliptic systems. SIAM Journal of Science Statistics in Computation , 1986 ,7 : 540,559 6 Costiner S , Ta ’asan S. Simultaneous multigrid techniques for nonlinear eigenvalue problems : solutions of the ( ) nonlinear Schrodinger2Poisson eigenvalue problem in two and three dimensions. Physics Review E3, 1995 , 52 : 1181,1192 7 λMcCormick S. A mesh refinement method for Ax = B x . Mathematics Computation , 1981 , 36 :485,498 8 Cai Z , Mandel J , McCormick S. Multigrid methods for nearly singular linear equations and eigenvalue problems. SIAM. Journal Numerical Analysis , 1997 ,34 :178,200 9 Bank R E. Analysis of a multilevel inverse iteration procedure for eigenvalue problems. SIAM. Journal Numerical Analysis , 1982 ,19 : 886,898 10 D’Yakonov E G. Iteration methods in eigenvalue problems. Mathematical Notes , 1983 , 34 :945,953 11 D’Yakonov E G , Orekhov M Y. Minimization of the computational labor in determining the first eigenvalues of differential operators. Mathematical Notes , 1980 ,27 :382,391 12 Knyazev A V. Preconditioned eigensolvers. Electronics Transaction on Numerical Analysis , 1998 , 7 : 104, 123 13 ()electronic Neymeyr K. A geometric theory for preconditioned inverse iteration. I : Extreme of the Rayleigeh quotient . Linear 14 Algebra Application , 2001 ,322 :61,85 Neymeyr K. A geometric theory for preconditioned inverse iteration. II : Sharp convergence estimates. Linear 15 Algebra Application , 2001 ,322 : 87,104 16 Oliveira S. On the convergence rate of a preconditioned subspace eigensolver. Computing , 1999 ,63 : 219,231 哈克布思 W 著 ,林群译 . 多重网格方法 . 北京 :科学出版社 ,1988 AN EFFECTIVE AL GO RITHM FO R GENERAL IZED EIGENVAL UE PRO BL EMS Liu Hanbing Gong Guoqing Liu J ianshe ( )Traff ic and Transportation College , Jilin University , Changchun , 130025 An effective multi2grid method fit for the adaptive finite element analysis for Abstract generalized eigenvalue problems is proposed. The method utilizes sufficiently the results in the original meshes and takes interpolation method or least2square approximation to obtain new eigenvectors in the changed meshes. And then it accomplishes the solution of the generalized eigenvalue problems by multi2grid iteration procedure . In the smooth process we select Conjugate Gradient method other than the traditional Gauss2Seidel or J acobi iterations in order to improve the rate of convergence ; in the coarse grid correction we deduct a formula to solve approximately the errors of the eigenvectors. It combines the grid discretization with the numerical solution and establishes a fast reanalysis procedure after remeshing. The numerical results show that the method is computationally simple and feasible as compared with the traditional finite element method , and can be used as a particularly effective tool for the generalized eigenvalue problems by adaptive finite element method. Key words multi2grid iteration ,generalized eigenvalue problem ,adaptive FEM