(School of Information Science & Technology,Xiamen University Tan Kah Kee College,Zhangzhou 363105,China)
backward error; polynomial eigenvalue problem; vibration analysis; quadratic eigenvalue problem
DOI: 10.6043/j.issn.0438-0479.201705013
备注
给出了一类形如(λkAT+λlA)z=0(A为稀疏矩阵)的矩阵方程的多项式特征值问题向后误差分析.并通过高速列车的振动分析中的一类二次特征值问题(λ2AT+λQ+A)z=0的例子,应用该方法讨论此类二次特征值问题的向后误差.
This article will give backward errors to polynomial eigenvalue problems like(λkAT+λlA)z=0),where A is a sparse matrix.For example,in studying the vibration of fast trains,we encounter a palindromic quadratic eigenvalue problem(QEP)λ2AT+λQ+A)z=0.We will discuss the backward error for this problem with the new method in this article.
引言
本文中考虑如下多项式特征值的向后误差:
(λkAT+λlA)z,(1)
其中,k,l∈Z,A是一个n×n复矩阵,它只在右上角有一个m×m块为非零块且m≤n/2,λ和z为特征值和特征向量.
本文中引用了
参考文献[1]中多项式特征值问题向后误差的定义和标记,其重要性在于研究算法的稳定性和质量,它与向前误差、条件数之间有如下关系:向前误差≤向后误差×条件数.可以看出向后误差与条件数会直接影响到数值结果的误差估计.
扰动理论和向后误差分析有广泛应用,如:线性系统[2]、最小二乘问题、一般特征值问题和广义特征值问题[3-4],以及物理学中的过阻尼物理系统[5-6].近年来有很多文献对向后误差的概念进行了阐述[1,3-4],Rigal等[2]给出了线性系统的向后误差分析,Fraysse等[3-4]将向后误差分析推广到了多项式特征值问题,Duffin[6]介绍了一些特殊矩阵的向后误差,Stewart等[7]也提到了一些特殊矩阵的向后误差.
本文中考虑的多项式特征值(λ4AT+λlA)z=0来自于二次特征值问题(λ2AT+A)z=0,本文中把它推广到(λkAT+λlA)z=0,在后面的应用举例中,将以铁轨在高速列车通过时的振动问题的有限元模型中产生的二次特征值问题为例,使用本文中的方法求出该问题的向后误差.
二次特征值问题在许多领域有丰富的应用[8-13],如汽车制动系统中的有限元系统[11],地震工程[9],保守结构系统和非保守结构系统分析[12-13],以及最小二乘问题[10]等.但是二次以上特征值问题的向后误差分析并没有很好地得到解决,很多问题难以求出向后误差.
本文中的例子是铁轨在高速列车通过时的振动问题的有限元模型中产生的二次特征值问题[14-16].该问题对文献[17]的对称多项式特征值问题以及文献[18-22]有很大的推动作用.Chu等[18]为该问题建立了有限元模型,Guo等[23]对该问题提出一种高精度的保持特征值结构的doubling算法,Lu等[24-25]对该算法进行改进,减少了计算量和计算时间.袁飞等[26]使用回文特征值问题向后误差的分析方法给出了该问题在文献[22]的定义下的一种向后误差.由于
参考文献[1]定义的多项式特征值问题向后误差更为精准和实用,本文中将会用新的方法,给出该二次特征值问题在
参考文献[1]定义下的向后误差.
1 相关定义定理
1.1 一些基本定义考虑如下非线性特征值问题:
P(λ)z=0,(2)
其中,矩阵P(λ)的元素为含λ的多项式,形式为
P(λ)=λsAs+λs-1As-1+…+A0,
其中,Al∈Cn×n,l=0,1,…,s,称P为λ-矩阵[1].
在本文中,ΔAl表示Al的误差.为方便表示,令
ΔP(λ)=λsΔAs+λs-1ΔAs-1+…+ΔA0.
对于复数λ,令
sign(λ)={(-overλ)/λ, λ≠0,
0, λ=0,
本文中使用的2-范数与文献[1]相同,为‖x‖2=(x*x)1/2,‖A‖2=max{‖Ax‖2:‖x‖2=1}.
1.2 多项式特征值问题向后误差的定义本文中引用了文献[1]中多项式特征值问题向后误差的定义和标记.
对于方程(2)的一对近似解((～overx),(～overλ)),矩阵El(l=0:s)为任意给定的矩阵,定义向后误差η((～overx),(～overλ))如下:
η((～overx),(～overλ)):=min{:(P((～overλ))+ΔP((～overλ)))(～overx)=0,
‖ΔAl‖2≤‖El‖2,l=0,1,…,s}.(3)
由此η((～overx),(～overλ))可由文献[23]中的定理求出.
定理1[23] 向后误差η((～overx),(～overλ))可由如下公式求出:
η((～overx),(～overλ))=(‖r‖2)/((～overα)‖(～overx)2‖),(4)
其中r=P((～overλ))(～overx),(～overα)=∑sl=0|(～overλ)l|‖El‖2.
从定理1可以看出,求多项式特征值问题向后误差的关键在于求出每个‖ΔAl‖2对应的上界‖El‖2.
作为对比,给出文献[26]中定义的二次特征值问题的向后误差.
对于P(λ)=λ2AT+λQ+A,定义A0=A,A1=Q,A2=AT.假设(ξ,ν)是P(λ)的一个近似特征值对,则定义高速列车震动分析问题的向后误差ΔF:
ΔF=min(∑1l=0‖ΔAl‖2)1/2,
其中,矩阵A0,A1满足(∑2l=0(Al+ρlΔAl)ξl)ν=0,常量ρ2-l=ρl≥0,通常取自Al的某些范数,等式右边的范数可以为2-范数或者F-范数.
2 向后误差的求法
在本节中,将讨论多项式特征值问题(1)的向后误差:
(λkAT+λlA)z=0,
其中k,l∈Z,A是一个n×n复矩阵,它只在右上角有一个m×m块为非零块且m≤n/2,λ和z为特征值和特征向量.
在1.2的讨论中可以知道,求该问题向后误差的关键在于求出‖ΔA‖2的上界.由本文中对2-范数的定义可以得到
‖A‖2≤‖A‖F.
于是,只需要求出‖ΔA‖F的下界,即可得出‖ΔA‖2的上界.本文中利用如下定理来求式(1)的向后误差.
定理2 对于多项式特征值问题:
(λkAT+λlA)z=0,
其中k,l∈Z,A是一个n×n复矩阵,它只在右上角有一个m×m块为非零块且m≤n/2,λ和z为特征值和特征向量,该问题对应于一对近似解((～overλ),(～overz))的向后误差可以表示为
‖ΔA‖2≤21/2/2‖Z+a‖2.
其中
Z=[(～overz)T3 0 … 0
0(～overz)T3 … 0
0 0 …(～overz)T3 0m×m2
0m×m2 (～overz)T1 0 … 0
0(～overz)T1 … 0
0 0 …(～overz)T1],
α=[α3
α1],{A12(～overz)3=-a3,
AT12(～overz)1=-a1,
Z+为Z的Moore-Penrose逆矩阵[27].
证明 首先把式(1)写成分块矩阵的形式:
[λk(0 0
AT12 0)+λl(0 A12
0 0)]z=0,(5)
其中A12是m×m块为非零块且m≤n/2.把z写为如下形式:
z=(z1 z2 z3)T,
其中z1和z3是m×1向量,z2是剩余部分.
于是式(5)可以表示为
{λlA12z3=0,
λkAT12z1=0.(6)
对于方程(5)的一对近似解((～overλ),(～overz)),A12的误差为ΔA12,由向后误差的定义式(3)可以得到:
[(～overλ)k(0
AT12+ΔAT12 0)+
(～overλ)l(0 A12+ΔA12
0 0)](～overz)=0,(7)
把它写成类似式(6)的形式:
{(～overλ)l(A12+ΔA12)(～overz)3=0,
(～overλ)k(AT12+ΔAT12)(～overz)1=0.(8)
设{A12(～overz)3=-a3,
AT12(～overz)1=-a1,则可以得到关于ΔA12的方程组:
{ΔA12(～overz)3=a3,(9)
ΔAT12(～overz)1=a1.(10)
只需求出‖ΔA12‖F的下界即可得到‖ΔA‖F的下界.令
ΔA12={aij},i,j=1,2,…,m,
ai:表示ΔA12的第i行; a:j表示ΔA12的第j列,则式(9),(10)可分别写成如下形式:
[(～overz)T3 0 … 0
0(～overz)T3 … 0
0 0 …(～overz)T3][aT1:
aT2:
aTm:]=a3,(11)
[(～overz)T1 0 … 0
0(～overz)T1 … 0
0 0 …(～overz)T1][aT:1
aT:2
a:m]=a1,(12)
其中,0都表示1×m的一行0,那么整个系数矩阵(11),(12)均是m×m2阶的.
把式(11)和(12)整合起来写成如下方程:
Zx=a,(13)
其中
Z=[(～overz)T3 0 … 0
0(～overz)T3 … 0
0 0 …(～overz)T3 0m×m2
0m×m2 (～overz)T1 0 … 0
0(～overz)T1 … 0
0 0 …(～overz)T1],
x=[a1: a2: …am: aT:1 aT:2 … aT:m]T,
a=[a3
a1],{A12(～overz)3=-a3,
AT12(～overz)1=-a1.
令x*=Z+a,Z+,为Z的Moore-Penrose逆矩阵[27].那么min=ΔA12=2F=1/2=x*=22=1/2=Z+a=22.则得到了如式(3)中定义的多项式特征值问题向后误差:
=ΔA=2≤(min=ΔA=2F)1/2=
(min=ΔA12=2F)1/2=21/2/2=Z+a=2.
证毕.
3 高速列车振动分析问题的向后误差
在高速列车的振动分析[23-25]中,P(λ)=λ2AT+λQ+A,对应特征方程为:
(λ2AT+λQ+A)z=0.(14)
其中A和Q为n×n复矩阵,且具有如下特殊结构:
Q=Kt+iωDt-ω2Mt∈Cn×n,
A=Kc+iωDc-ω2Mc∈Cn×n,
其中,i为虚数单位,ω>0为外力的频率,Kt、Dt、Mt、Kc、Dc、Mc都是m×m的分块实矩阵,每个块有k×k个元素
.
A和Q都是m×m的分块矩阵,每个块有k×k个元素,即n=m×k; 此外,Q是块三对角阵,A只有位于(1,m)位置的一个块为非零块A13.
本文中仅讨论矩阵A的向后误差.下面将用第2节中的方法讨论此二次特征值问题在第1节定义下的向后误差.
按照第1节中的定义,对于特征方程(14)的一对近似解((～overx),(～overλ)),矩阵El为任意给定的矩阵,定义向后误差η((～overx),(～overλ))如下:
η((～overx),(～overλ)):=min{:((～overλ)2AT+A+(～overλ)2ΔAT+ΔA)(～overx)=
0,=ΔA=2≤=El=2}.(15)
该问题等价于本文中第2节中研究的特征多项式(1)在k=2,l=0时的向后误差问题.故可直接用第2节中的方法解决这个问题.
由第2节的定理1,可以得到问题(15)的向后误差为
=ΔA=2≤21/2/2=Z+a=2.
其中,
Z=[(～overx)T3 0 … 0
0(～overx)T3 … 0
0 0 …(～overx)T3 0k×k2
0k×k2 (～overx)T1 0 … 0
0(～overx)T1 … 0
0 0 …(～overx)T1],
a=[a3
a1],{A12(～overx)3=-a3,
AT12(～overx)1=-a1.
其中Z+为Z的Moore-Penrose逆矩阵[27].
4 数值实验
下面给出数值结果,本节使用的矩阵与文献[23-25]给出的模型相同.用于数值实验的3个矩阵的大小分别为(k,m)=(159,11),(303,19),(705,51).系数矩阵A和Q的形式由本文引言部分给出.振动频率ω的取值为1 000.
对于每组系数矩阵A和Q对应的二次特征值问题P(λ)z=0,本文中首先使用文献[5]中的算法求出它的全部近似特征值和特征向量,然后随机从中取出一个近似特征值对(λ,z).使用文献[26]中定理3方法求出该问题的向后误差表示为ΔF,用本文中的方法求出的向后误差表示为ΔA,结果如表1所示,可见ΔA的误差更小.
表1 数值实验结果
Tab.1 Numerical experiment results(k,m)Λ ΔF =ΔA=2(159,11)0.87-0.12i 7.8e-08 7.1e-09(303,19)0.77-0.18i 1.2e-07 9.8e-08(705,51)0.60-0.27i 3.6e-09 1.3e-09
5 结 论
本文中介绍了多项式特征值问题向后误差相关的概念,并且使用文献[1]中多项式特征值问题向后误差的定义,求出了特征值问题(λkAT+λlA)z=0的向后误差.铁轨在高速列车通过时的振动问题的有限元模型中产生的二次特征值问题(λ2AT+λQ+A)z=0属于回文二次特征值问题,所以
参考文献[26]中处理回文二次特征值问题的方法来分析该二次特征值问题的向后误差.由于文献[1]定义的多项式特征值问题向后误差更为精准和实用,本文中使用新的方法,给出该二次特征值问题在
参考文献[1]定义下的向后误差.在本文中最后,对于
参考文献[23-25]给出的3个铁轨在高速列车通过时的振动问题的有限元模型,给出了一些数值实例,并且与
参考文献[26]给出的向后误差进行对比.
- [1] TISSEUR F.Backward error and condition of polynomial eigenvalue problems[J].Linear Algebra and Its Applications,2000,309:339-361.
- [2] RIGAL J L,GACHES J.On the compatibility of a given solution with the data of a linear system[J].Journal of the ACM,1967,14(3):543-548.
- [3] FRAYSSE V,TOUMAZOU V.A note on the normwise perturbation theory for the regular generalized eigenpro-blem[J].Numerical Linear Algebra with Applications,1998,5(1):1-10.
- [4] HIGHAM D J,HIGHAM N J.Structured backward error and condition of generalized eigenvalue problems[J].SIAM Journal on Matrix Analysis and Applications,1998,20(2):493-512.
- [5] LANCASTER P,SNEDDON I N,STARK M,et al.Lambda-matrices and vibrating systems[M].Oxford:Pergamon Press,1966:187-190.
- [6] DUFFIN R J.A minimax theory for overdamped networks[J].J Rat Mech Anal,1954,4(2):221-233.
- [7] STEWART G W,SUN J.Matrix perturbation theory[M].London:Academic Press,1990:1-94.
- [8] MACKEY D S,MACKEY N,MEHL C,et al.Numerical methods for palindromic eigenvalue problems:computing the anti-triangular schur form[J].Numerical Linear Algebra with Applications,2009,16(1):63-86.
- [9] GANDER W,GOLUB G H,VON MATT U.A constrained eigenvalue problem[J].Linear Algebra and Its Applications,1989,114/115:815-839.
- [10] CLOUGH R W,MOJTAHEDI S.Earthquake response analysis considering non-proportional damping[J].Earthquake Engineering & Structural Dynamics,1976,4(5):489-496.
- [11] KOMZSIK L.Implicit computational solution of genera-lized quadratic eigenvalue problems[J].Finite Elements in Analysis & Design,2001,37(10):799-810.
- [12] SMITH H A,SINGH R K,SORENSEN D C.Formulation and solution of the non-linear,damped eigenvalue problem for skeletal systems[J].International Journal for Numerical Methods in Engineering,1995,38(18):3071-3085.
- [13] ZHENG Z C,REN G X,WANG W J.A reduction method for large scale unsymmetric eigenvalue problems in structural dynamics[J].J Sound and Vibration,1997,199(2):253-268.
- [14] KOMZSIK L.Implicit computational solution of genera-lized quadratic eigenvalue problems[M].Manuscript:the MacNeal-Schwendler Corporation,1998:56-72.
- [15] SMITH H A,SINGH R K,SORENSEN D C.Formulation and solution of the non-linear,damped eigenvalue problem for skeletal systems[J].International Journal for Numerical Methods in Engineering,1995,38(18):3071-3085
- [16] IPSEN I C F.Accurate eigenvalues for fast trains[J].SIAM News,2004,37(9):1-2.
- [17] MACKEY D S,MACKEY N,MEHL C,et al.Structured polynomial eigenvalue problems:good vibrations from good linearizations[J].SIAM Journal on Matrix Analysis and Applications,2006,28(4):1029-1051.
- [18] CHU E K W,HWANG T M,LIN W W,et al.Vibration of fast trains,palindromic eigenvalue problems and structure-preserving doubling algorithms[J].Journal of Computational and Applied Mathematics,2008,219(1):237-252.
- [19] HUANG T M,LIN W W,QIAN J.Structure-preserving algorithms for palindromic quadratic eigenvalue problems arising from vibration of fast trains[J].SIAM Journal on Matrix Analysis and Applications,2008,30(4):1566-1592.
- [20] KRESSNER D,SCHRÖDER C,WATKINS D S.Implicit QR algorithms for palindromic and even eigenvalue problems[J].Numerical Algorithms,2009,51(2):209-238.
- [21] LANCASTER P,PRELLS U,RODMAN L.Canonical structures for palindromic matrix polynomials[J].Oper Matrices,2007,1(4):469-489.
- [22] LI R C,LIN W W,WANG C S.Structured backward error for palindromic polynomial eigenvalue problems[J].Numerische Mathematik,2010,116(1):95-122.
- [23] GUO C H,LIN W W.Solving a structured quadratic eigenvalue problem by a structure-preserving doubling algorithm[J].SIAM Journal on Matrix Analysis and Applications,2010,31(5):2784-2801.
- [24] LU L,YUAN F,LI R C.A new look at the doubling algorithm for a structured palindromic quadratic eigenvalue problem[J].Numerical Linear Algebra with Applications,2015,22:393-409.
- [25] LU L,YUAN F,LI R C.An improved structure-preserving doubling algorithm for a structured palindromic quadratic eigenvalue problem[R].Arlington,UT:University of Texas at Arlington,2014.
- [26] 袁飞,卢琳璋,李仁仓.一类二次特征值问题的向后误差分析[J].厦门大学学报(自然科学版),2016,55(1):97-102.
- [27] PENROSE R.A generalized inverse for matrices[J].Mathematical Proceedings of the Cambridge Philosophical Society,1955,51(3):406-413.
