Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
1/26
Organization (S): EDF/IMA/MN
Handbook of Référence
R5.01 booklet: Modal analysis
Document: R5.01.02
Algorithms of resolution for the problem
quadratic with the eigenvalues
Summary:
In this document, we fix the theoretical framework of the methods of search for eigenvalues of the problem
quadratic which is developed in the code of Aster mechanics.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
2/26
Contents
1 Introduction ............................................................................................................................................ 4
1.1 Position of the problem ....................................................................................................................... 4
1.2 Properties of the matrices ................................................................................................................... 4
1.3 Problem with the eigenvalues associated ........................................................................................... 4
1.4 Some traditional particular cases. ............................................................................................. 5
2 Reductions with a linear form ............................................................................................................ 6
2.1 Reductions with a linear form ...................................................................................................... 6
2.1.1 A particular choice for matrix Z: ±M ............................................................................. 7
2.1.2 Particular case of the positive matrix M definite ...................................................................... 7
2.1.3 Case of the symmetrical reduction for a singular matrix M ............................................. 8
2.2 Property of orthogonality of the clean vectors ............................................................................... 9
3 Method of determinant ...................................................................................................................... 10
3.1 General. ................................................................................................................................... 10
3.2 Method of Muller ......................................................................................................................... 10
3.2.1 Development of the method ............................................................................................ 10
3.2.2 Convergence of the method ................................................................................................ 12
3.2.3 Application of the method for search of eigenvalues ........................................ 12
3.2.3.1 Development ........................................................................................................ 12
3.2.3.2 Cost of the method in term of factorization ......................................................... 13
4 Methods of iteration reverses ................................................................................................................ 14
4.1 Method of opposite iteration suggested by Wilkinson ..................................................................... 14
4.2 Alternative developed by Jennings ................................................................................................ 14
4.3 The algorithm of iteration reverses of Aster .......................................................................................... 15
4.3.1 Implementation ..................................................................................................................... 15
4.3.2 Criterion of stop ....................................................................................................................... 15
5 Lanczos Method applied to the quadratic problem ....................................................................... 16
5.1 Choice of a problem to approach ................................................................................................. 16
5.1.1 Form reverses standard ....................................................................................................... 16
5.1.2 Strategies of shift ......................................................................................................... 16
5.2 Method of approximation .............................................................................................................. 17
5.2.1 Approximate problem and algorithm of Lanczos ..................................................................... 17
5.2.2 Choice of a pseudo scalar product ...................................................................................... 19
5.3 Application to quadratic Problème ............................................................................................ 19
5.3.1 Spectral operator ................................................................................................................ 19
5.3.2 Operator of pseudo scalar product .................................................................................. 20
5.3.3 Cost of the Lanczos phase .................................................................................................... 21
5.4 Implementation in Aster ............................................................................................................ 21
5.4.1 Parameters of the implementation ......................................................................................... 22
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
3/26
5.4.2 Under space of approximation ............................................................................................... 22
5.4.3 Strategy of reorthogonalisation ........................................................................................... 22
5.4.4 Implementation of the phase Lanczos ................................................................................... 22
5.4.5 Restoration of the approximations for the quadratic problem ........................................... 23
6 Bibliography ......................................................................................................................................... 24
Appendix 1 Interprétation of the complex eigenvalues ........................................................................ 25
Appendix 2 Réductions linear ................................................................................................................ 26
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
D
2
U
(
2n
T)
M
D
T
2
.
U
J
+
G
M
(
.
(
K
J
=
X
E
T
+
X
T
1
(
C
+
)
)
(
=
K
=
+
0
E
+
)
(
U
K
(
+
T
E
)
= F (T)
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
4/26
1 Introduction
1.1
Position of the problem
Dynamic analysis or the study of the stability of the balance of a mechanical structure led, in
tally of the linearized theory, to solve the matric differential equation of the second command
where
M
is the matrix of mass and inertia of the structure,
G
is the matrix induced by the gyroscopic effect (case of the revolving machines),
is a significant real parameter disk speed,
C
is the matrix of damping induced by dissipative forces.
K
is the matrix of rigidity of the structure,
E
is the matrix of viscous damping interns structure,
F
is the external force (which is null in the case of the search for balance).
1.2
Properties of the matrices
The matrices considered are with real coefficients.
Classically, one considers:
M
is symmetrical (semi) definite positive,
G
is antisymmetric,
C
is symmetrical,
K
is symmetrical not necessarily definite positive,
E
is antisymmetric.
Consequently, one realizes that the simultaneous presence of the matrix of damping and the matrix
of gyroscopic effect destroys the symmetry or the antisymetry of the term speed; pareillement
internal damping producing an antisymmetric matrix, destroyed the property of symmetry of
stamp rigidity.
In addition, the introduction of linear relations modifies the character of positivity of the matrices:
stamp K is indefinite (with positive or negative eigenvalues).
1.3
Problem with the eigenvalues associated
The sought solutions are form (separation of the variables of space and time).
U (T) = E t.x with IC and X ICN
What leads us to the quadratic problem with the eigenvalues according to:
the solution U (T) can be rewritten in the form:
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
D
U
dt
(
X
0
X
D
G
)
K
=
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
5/26
that one can put in matric form:
U (T) = [X] [E T] [K]
where
[X]
is the modal matrix
[N X 2n]
[E T] is a diagonal matrix
[2nx2n]
[K]
is a matrix unicolonne
[2nx1]
Proposal:
The modal matrix [X] perhaps used like stamps transformation to uncouple N
equations of the quadratic problem of origin, its 2n columns are not linearly
independent.
Note:
One determines [X], using the two following identities:
U (0) = [X] [K]
1.4
Some traditional particular cases.
These particular cases were lengthily studied (cf for example [MEI.67], [ROS.84])
M
G
C
K
E
values
vectors
clean
clean
rights
lefts
conservative
0
0
0
IR
IR
conservative
0
0
I IR
=
IC
gyroscopic
=
dissipative
0
0
IC
=
IC
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
Z
U
0
0
U
1
0
X
Z
U
-
U
With
-
K
- 1
M
C
X
M
C
=
0
K
.
B
Z
.
=
Z
=
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
6/26
2
Reductions with a linear form
One is interested in the possibilities of reduction of the quadratic problem in a generalized problem
equivalent.
The principles of reduction are applicable to unspecified matrices, and if we consider it
problem in (M, C, K) is only to simplify the talk.
2.1
Reductions with a linear form
There are several traditional methods to transform the quadratic problem into a problem
generalized with the eigenvalues.
We will develop the method which consists in introducing speed like auxiliary variable:
For that we introduce an additional equality:
Z (T) - Z (T) = 0
where Z is a matrix not identically null.
Our initial system can then be rewritten in the matric form in a space of dimension
double of initial space
And the problem with the eigenvalues associated is in its opposite form
By supposing K and Z regular, the problem with the eigenvalues generalized thus obtained can
formally to put itself in the standard form:
Notation: one poses y = X
The standard form associated the quadratic problem is independent of regular matrix Z
chosen.
Definition
One will call linear reduction of the quadratic problem, any generalized problem of which all them
clean elements ((X, y)) check:
,
(X) is clean solution of the quadratic problem,
·
,
and My = MX (condition of coupling)
·
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
M
y
0
1
0
y
0
T
1
T
QT
C
Q1
M
T
1
QT
K
Q1
X
-
Q
(
Q
.
,
With
M
C
-
X
C
T
Q
.
0
=
0
K
M
K
)
.
Q
.
(
)
=
+
T
)
= 0
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
7/26
Proposal:
If ((X, y) T) is a solution of the linearized equation and if Z is regular then (, X) is solution of
Pr,
oblème quadratic.
This result is immediate.
2.1.1 A particular choice for matrix Z: ±M
one takes for matrix Z the matrix - M of the initial system, consequently the linearized system (A.Z = B.z)
·
is written:
µ
If the matrices M, C and K are symmetrical and Si the matrix M is regular, this linear reduction
allows to build a symmetrical generalized problem.
one takes for matrix Z the matrix M of the initial system, consequently the linearized system (A.Z = B.z)
·
is written:
µ
If M and K are definite positive then B is also, but matrix A is not symmetrical.
Proposal:
If one chooses Z = ± M (matrix not identically null) and if (, (X, y) T), with X not no one, is solution
clean of the equation generalized then (, X) is solution propr
E of the quadratic problem.
2.1.2 Particular case of the positive matrix M definite
If the matrix of mass is definite positive then it a decomposition of Choleski admits
M = QT.Q
If one introduces the linear transformation U (T) = Q-1.q into the basic equation of the quadratic problem
after pre-multiplyhaving multiplied it by Q-T = (Q-1) T, one obtains:
and while defining v =
the preceding equation, leads us to the standard problem with the values
clean A. v = .v with
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
T
M
K
C
X
M
11
0
1
1
1
M
Q
M
R
With
0
B
T
y
M
=
0
Q
R
-
- M
M
R
R
T
R
22
=
M
Q
Q
=
Q
=
=
Q
=
M
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
8/26
2.1.3 Case of the symmetrical reduction for a singular matrix M
In this paragraph we suppose the matrices M, C and K symmetrical.
If M is singular and if the eigenvalue is semi simple then it Q exists an orthonormal matrix
such as.
0
Indeed, the eigenvalue being semi-simple, the core of M admits a base made up of vectors
clean. This base is C 0
omplétée and orthonormalized. The matrix Q admits for vectors columns them
vectors of this base.
Let us note that M11 is symmetrical (bus M is) and regular.
One introduces then a regularization of the matrix M represented by the matrix MR. and defined by:
where M22 is a regular symmetrical matrix (for example the matrix identity).
Property:
The generalized problem associated the matrices
is a symmetrical linear reduction of the quadratic problem.
Demonstration:
symmetry is ensured by construction.
· the “tildées” matrices check
·
thus
, that is to say still
However if
associated the eigenvalue of the generalized problem one is a clean vector has
MR. y = MX, therefore My = M
MX. The condition of coupling is established.
Let us suppose that under-vector X of the clean vector
that is to say no one. Then the equation
·
MR. y = MX led to y = 0 bus MR. is regular, which is absurd.
The condition X 0 is thus established.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
X
T
M
1
C
K
X
1
I
T
y
- M
R
- 1
M
0
X
I
,
M
J
-
y
.
I
,
y
I
I
+
J
X
X
C
J
has
I
J
X
I
+
B
ij
has
I
B
+
I
ij
K
I
ij
=
With
B
.
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
9/26
Note: Interpretation of the regularized matrix MR.:
This symmetrical linear reduction is distinguished from that presented previously by
substitution of M per MR. in matrix A.
If M is singuliére, this substitution forces the vector to have a component there
null in the core Mr. Ainsi the symmetrical linear reduction led to a generalized problem
posed in under space of
IR2n such as y ker M - {0}.
If K is regular, the preceding generalized problem admits a form reverses standard
equivalent.
2.2
Property of orthogonality of the clean vectors
It is immediate to show that the complex clean modes check the properties of orthogonality
exits of the following generalized problem if the matrices A and B are symmetrical:
If one develops the preceding expressions, by taking account of the linear reductions used one
obtains the expressions:
Note:
The first equality is independent of the regularity of M,
· the second equality is obtained directly if M is regular and is established if not in
·
using the basic change of the regularization of Mr.
The modes of the quadratic problem are thus not M, C or K orthogonal.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
B
0
Z
I
2
1
B
2
Z
I
)
1
F
(
B
0
Z
I
1
2
=
-
+
B
2
Z
I
)
=
1
=
-
2
1
2
F
+
(
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
10/26
3
Method of determinant
3.1 General.
The search for zero of the characteristic polynomial, a quadratic problem to the eigenvalues,
pose the problems inherent in the polynomials of the complex variable with complex values.
Not having relation of command in IC, usual methods, at two points, of dichotomy type
or secant are inapplicable.
We present here the “popular” method more for search of zero of polynomials of
variable complexes with complex values, the method of Muller.
3.2
Method of Muller.
The method suggested by Muller [MUL.56] is an iterative method using as curve
of interpolation a parabola with horizontal axis.
This method is relatively easy to implement but it lends itself badly to search
zeros of real functions with real roots, because it plunges the interpolation in the complex plan and this
even on the basis of actual values.
Its interest is related to the class of this method “the methods by curves of interpolation”, namely:
the safety of the method by dichotomy, since search is carried out in a ball “
·
reducing " gradually,
that only the calculation of the function is necessary (calculation of not derived as in
·
method of Newton),
convergence is connected with a quadratic convergence.
·
The most widespread method of the methods by curves of interpolation is the method at two points known as
secant.
3.2.1 Development of the method
Let us note F (Z) = a0zn + a1zn-1 +… + year the algebraic equation which one seeks the zeros, have are them
complexes and we suppose a0 not no one.
The quadratic formula of interpolation of Lagrange gives us:
Li (F (Z)) = b0z2 + b1z +b2
and we consider the curve which passes by the last three points (reiterated):
(zi, F (zi)), (zi-1, F (zi-1)), (zi-2, F (zi-2))
and thus the coefficients b0, b1 and b2 check:
and while posing
and
and
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
L
I
F
2
)
(
F (
Z
I
Z
I
)
2
Z
2
F
)
F
1
(
(
I
)
2
G
Z
F (
Z
I
-
Z
I
)
2
+
=
±
F (
)
I
2
-
+
)
1
)
I
F
1
4
2
+
(
(
+
F
(
F
F
(
(
+
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
11/26
one can rewrite the formula of interpolation of Lagrange
the new point is:zi+1 = Z
One can solve the quadratic equation in, we obtain then:
by taking the reverse of the traditional solution of a quadratic equation:
with
From i+1, one obtains:
hi+1 = i+1 hi,
zi+1 = zi + hi+1 which is one zero of the equation.
The sign of the denominator is then taken of such kind that it is of larger possible module and thus of
such kind that i+1 is of larger possible module, and finally zi+1 will be the root nearest to
zi
Note:
Muller [MUL.56] proposes a process of initialization while using:
“arbitrary” values z0 = - 1., z1 =1., z2 = 0.
and values
year + year-1+ year-2
for F (z0)
year - year-1+ year-2
for F (z1)
year
for F (z2)
This choice of values then results in considering:
L2 (F (Z)) = = year + year-1z +an-2z2
who is an approximation of F (X) in the vicinity of the origin.
The advantage of this process of starting is that it does not require any evaluation of the polynomial F (X) and
that it is thus fast.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
Z
I
F
K
(
H
Z
F
K
I
1
)
(
Z
1
2
)
Z
I
F
H
L
T
K
=
(
Z
I
)
-
<
Z
Z
I
=
1
-
1
,
=
2
K
+
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
12/26
3.2.2 Convergence of the method
Proposal: To consider that convergence is assured as soon as
for given is one
acceptable criterion.
Let us show initially that during the unfolding of the algorithm three the reiterated successive ones are distinct.
If it were not the case, by supposing that zi = zi-1 and that xi 1st are reiterated for which one has
convergence and thus
.
If it is supposed that Z
I = zi-2, one would then have I and i+1 identically null and thus one would have zi+1= zi
from where contradiction.
Then into constant that the difference between two reiterated can only decrease, one obtains the result
announced.
3.2.3 Application of the method for search of eigenvalues
3.2.3.1 Development
That is to say to determine the eigenvalues of the system (2M+ K+C) X = 0,
We thus seek the zeros of the polynomial Ca
ractér
istic that we note
F (Z) = dét (z2M + zK + C)
To calculate in sequence the zeros of the polynomial, we use a technique of deflation.
Consequently the polynomial considered is:
where
are K zeros already calculated.
The use of the algorithm is immediate, and we benefit from this adaptation to formulate it
slightly different way.
Let us note
the value of the function to be interpolated at the point zi at the time of the search from the kth zeros.
It is then practical, for the implementation, to reveal a few intermediate quantities:
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
2
F
I
F
I
T
2
2
1
I
G
0
F
1
±
2
Z
+
,
G
1
I
=
-
F
T
0
L
.
H
R
1
I
,
+
=
-
E
Im
Z
I
I
1
0
1
4
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
13/26
reiterated the li+1 sought being solution of:
with
One deduces some:
the sign of the denominator being taken so that li+1 (and thus hi+1) are of smaller module
possible
consequently
To converge towards the eigenvalues with positive imaginary part (and thus at positive frequency), one
takes
Notice on deflation:
When the eigenvalues appear per combined pairs, it is also necessary to eliminate
combined found eigenvalues and which physically corresponds to a negative frequency
3.2.3.2 Cost of the method in term of factorization
With each iteration of the method one makes:
Matrix algebra part (evaluation of the characteristic polynomial)
a linear combination of three matrices
O (N2)
· a LDLt factorization of the combination die
O (nb2), B width of
·
bandage
the product of the diagonal terms divided by deflation
O (N)
·
Method part:
operations of the method with properly spoken
O (1)
·
Moreover it is advisable to add, as an assumption of responsibility (since it is a method at three points)
twice the evaluation part of the characteristic polynomial.
Broadly this method costs (i+2) factorizations, for a solution calculated out of I iterations.
The cost of this method is such, that it must be to hold for small systems (case of development)
or when it is important to have a very high degree of accuracy on the sought frequency.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
0
2
y
X
S
M
C
y
2
M
1
K
1
K
-
M
- 1
K
-
µ
2
)
- 1
y
- 1
K
-
M
C
S
X
=
.
M
C
D (
)
- 1
.
D
D (
)
-
2
µ
-
C
=
+
-
X
S
M
M
C
(
1
M)
1
+
y
-
2
1
2
-
+
K
K
-
K
2
=
0
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
14/26
4
Methods of opposite iteration
The method of iteration opposite extends immediately to the quadratic problem by using its form
“linearized”.
The linearization of the quadratic problem not being single, there exist several alternatives.
4.1
Method of iteration opposite suggested by Wilkinson
Wilkinson [WIL.65] proposes to bring back the quadratic problem to the following standard problem
By supposing M invertible, and by posing y = X
Being given an approximation of the sought eigenvalue, one can define the iterative process
according to:
Matric equation, we deduce a system from equations to two unknown (ys+1, xs+1) and by
combination of this system we deduce the expression from ys+1 and xs+1.
4.2
Alternative developed by Jennings
When M and/or K are singular, one obtains a stable and equivalent quadratic equation in
introducing an auxiliary parameter
for a spectral shift of.
While replacing by
in the quadratic equation, one obtains the quadratic problem in:
µ
While noting
the dynamic matrix which is regular by construction, one
bring back to the standard problem such as it is proposed by Wilkinson while posing y = X
µ
This equation is stable in the direction where for strictly positive.
This process thus makes it possible “to regularize µ <
“the 1
orders of magnitude of the matrices M, C, and K through
stamp dynamic D ().
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
X
N
T
C
X
N
M
1
D (
0
)
0
2
T
K
1
C
MM
±
1
4
y
N
X
0
-
.
-
=
y
N
X
M
D (
0
)
X
N
+
0
2 X
T
1
-
M
+
M
X
N
1
C
.
-
.
M
.
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
15/26
4.3
The algorithm of iteration reverses of Aster
4.3.1 Implementation
This algorithm is available in Aster by operator MODE_ITER_INV.
That is to say a 0 approximate value of the sought eigenvalue, one builds the dynamic matrix D () that
one factorizes in LDLt form.
One initializes the iterative process by the vectors according to:
x0 = {(1, 0)}
y0 = 0 x0
the iterative process to obtain the nth one reiterated:
Standardization of xn-1 and yn-1 to avoid the overflows of capacity:
·
Resolution of:
·
Calculation of yn starting from xn:
·
Evaluation
of
N:
·
This diagram can be put in the matric form:
4.3.2 Criterion
of stop
We use a simple result [GOH. &al.86] relating to the polynomials of matrices of the form:
D () = M2 + C +K
That is to say O an eigenvalue of the operator D () and xo an associated clean vector not no one checking D (O) xo
= 0 then we have the equality:
The criterion of stop is done on the relative variation of 0
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
K
0
X
1
1
C
M
With
0
- M
R
y
- M
0
=
-
=
B
B
Z
=
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
16/26
5
Lanczos method applied to the quadratic problem
In this chapter we suppose the matrices M, C and K real symmetrical, so that the problem
quadratic associated to a linear form symmetrical Az = Bz can be tiny room
where
MR. is a regular matrix deduced from M which coincides with M if the latter is regular.
One seeks to develop the arithmetic method preserving real one overall in order to obtain
a problem reduces real.
5.1
Choice of a problem to be approached
5.1.1 Form reverses standard
One seeks an approximation of the couples (, X) of clean parts of the quadratic problem which
correspond to the close eigenvalues
of a shift = +i complexes given.
The approximation of Galerkin of Pr
oblème spectral of an operator
S in
an S
ous-space of Krylov km
= span (r0, Sr0., Sm-1 r0) makes it possible to approach the clean couples of elements of the operator S
corresponding to the eigenvalues of larger modules.
The passage to the standard form reverses after spectral shift of generalized problem
preceding, provides the spectral problem:
(A - B)-1 B Z = Z
This problem admits the same clean vectors as the problem generalized and of the eigenvalues
bound to those of the problem generalized by the relation:
Thus, an approximation of Galerkin of the spectral problem of the operator:
S = (A - B)-1 B =
in under space of Krylov km = span (r0, Sr0., Sm-1 r0) provides the approximation of the couples
clean parts of the quadratic problem sought.
5.1.2 Strategies of shift
The operator S preceding being complex, it calls in a natural way the use of the arithmetic one
complex. It is nevertheless possible to use the arithmetic real one to approach the couples
clean elements in which one is interested. It is enough to use a real operator of the same vectors
clean and whose eigenvalues of larger module correspond to those of the problem
quadratic closest to.
Technique of double DEC
alage spectral by and proposed by Francis within the framework of the method
QR makes it possible to build such an operator
noted to him also S:
S = [(A - B) (A B)] - 1 B = [AB-1A - 2 A + I I2 B] - 1 B
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
1
2
1
1
µ
-
1
2
2i
1
-
1
2i
-
=
-
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
17/26
The major disadvantage of this technique is the filling of the matrix
AB-1A - 2 A + I I2 B if B is not diagonal.
Parlett and Saad [PAR. SAA.87] propose an alternative which uses the real part or the imaginary part
of operator (A - B)-1B:
S+ = Re [(A - B)-1B]
S = Im [(A - B)-1B]
whose respectively noted eigenvalues + and - are related to the eigenvalues of the problem
µ
µ
quadratic by the relations:
+ and carry out the maximum of their module in the vicinity of and.
µ
µ
This approach preserves the arithmetic real one overall and avoids the disadvantage of the technique of
Francis, if the calculation of a vector S± v, where v are real, is carried out into arithmetic complex.
By noting that:
Re [(A - B)-1B] = [(A - B)-1 + (A B)-1] B
Im [(A - B)-1B] =
[(A - B)-1 - (A B)-1] B
This approach can be interpreted like a technique of double shift summons, by opposition
with the approach of the double produced shift suggested by J-C.F. Francis.
5.2 Method
of approximation
From now on S will indicate one of the real operators S+ or S, one of its eigenvalues and P
µ
stamp of a scalar pseudo-product (i.e. a symmetrical bilinear form not necessarily
defined positive).
5.2.1 Approximate problem and algorithm of Lanczos
When S is car-assistant for the pseudo scalar product induced by P, i.e.
(U, AV) P = (With, v) P U, v IR2n,
the method of Lanczos is used to generate a base of under space of Krylov km.
The spectral problem is then approached by projection P-orthogonal on km and the problem reduces thus
obtained is represented in the base of the vectors of Lanczos by a real matrix tridiagonale of command
Mr.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
J
1
2
1
2
\
1
Q
1
R
J
T
m
R
0
, R
J
P
+
\
m
=
,
R
0
P
1
R
J
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
18/26
The method of Lanczos extended to the scalar pseudo-products is defined by the formulas of
recurrence following cf [R5.01.01]:
ro = arbitrary
·
O =, qo =
·
0
0
1 = Sign ((ro, ro) P)
1 =
0
for J = 1, 2,…, m
·
J
=
(qj, Sqj) P
rj
=
Sqj - J J qj- j-1 J qj-1
j+1
=
sign ((rj, rj) P)
j+1 =
qj+1 =
If one notes Qm the matrix 2nxm vectors of Lanczos qj these formulas write in the form
matric real following:
Q T
m P Qm = Jm = diag (1,…, m)
S Qm - Qm Jm Tm = m+1 m+1 rm EM T with EM T = (,…, 1)
0
0
Tm the real matrix tridiagonale symmetrical m X m:
The product Jm Tm is a nonsymetric matrix tridiagonale as soon as Jm is not proportional to
identity.
The two preceding matric relations make it possible to write:
Q T
m P [S Qm - QmJmTm] = 0
The application of this relation to a couple ((m), S (m)) IC X ICm of clean elements of the operator
represented by the matrix J
µ
MTM gives:
Q T
(m)
m P [S Qm S (m) -
Qm S (m)] =
µ
0
who characterizes the couple ((m), Z (m) = Qm S (m)) IC X IC2n like approximation of Galerkin by
µ
projection P-orthogonal on km of a clean couple of element of the operator S.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
0
D
X
1
0
C
U
4
)
1
M
R
I
E
m
(
Z
With
B
M
R
0
y
B
- 1
M
0
M
R
C
- 1
M
U
4
)
-
M
R
- 1
=
-
- 1
M
- 1
+
M
(
I
X
m
=
-
D
-
R
E
(
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
19/26
5.2.2 Choice of a pseudo scalar product
The symmetry of matrices A and B ensures
that bilinear forms associated the matrices Re [(A - B)-1] - 1, Im [(A - B)-1] - 1 and B is
scalar pseudo-product,
that the operator left real S+ (respectively imaginary part S) is autoadjoint for
scalar pseudo-product induced by Re [(A - B)-1] - 1ou by B (respectively by Im [(A - B)-1] - 1 or by
B).
The pseudo-products scalar induced by Re [(A - B)-1] - 1 and Im [(A - B)-1] - 1 is extensions of
scalar product used in the alternative of Pipano-
Neuman of the algor
ithme Lanzcos.
If the matrix M is singular, the pseudo scalar product induced by B makes vectors
where y Ker M, of the quasi-null vectors. This disadvantage exists them also in the case of
pseudo-Pr
oduits scalar of Pipano-Neuman type (in particular if the matrix of the scalar pseudo-product
admits eigenvalues of null sum), but the occurrence of such an event is rare in
practical.
5.3
Application to quadratic Problème
The use of this approach requires, with the method of Lanczos, the calculation of the real vectors Sz and
Pz for Z IR2n.
5.3.1 Operator
spectral
That is to say the dynamic matrix D () associated the spectral shift = + I defined by:
2
D () =
M + C + K
If D () is regular then the operator (A complexes - B)-1B can be written in the form
The calculation of Sz for S = Re [(A - B)-1B] and S = Im [(A - B)-1B] can be carried out without destroying
structure digs matrices if the rear one
ithmetic complex be
T partially used in the algorithm:
Preparation into arithmetic complex
·
to form
D () = 2 M + C + K
- to factorize
D (
) under
form LDLT
-
Calculation of Sz
·
u1 = Cx u2 = MX u3 = My
in IR
- u4 = D () - 1 u1+ u2 +u3
in IC
-
according to the choice of the operator one obtains:
-
S+ Z = Re [(A - B)-1B] Z =
S Z = Im [(A - B)-1B] Z =
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
0
D
1
K
X
1
0
C
1
2
2
M
1
N
1
M
M
C
K
With
B
=
M
R
-
2
O
y
- 1
M
0
-
0
M
R
C
- M
- 1
B
B
ù
X
-
- 1
M
- 1
+
M
- M
C
+
- 2
M
,
IR
K
+
Z
y
E
D
K
Z
y
R
C
=
M
=
B
-
Z
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
20/26
5.3.2 Operator of pseudo scalar product
Choice P = B
·
This choice is valid for the operators obtained in real left approach or part
imaginary.
Calculations can be carried out without assembly of B and into arithmetic real.
Choice P = Re [(A - B)-1] - 1
·
This choice corresponds to the partly real approach of the operator S.
If the dynamic matrix D (), where is the real part of is regular then the real operator
(A - B)-1B is defined by:
,
and the scalar pseudo-product is written:
Re [(A - B)-1] - 1 = (A - B) + 2 B (A - B)-1B
Consequently the calculation of Pz can be carried out like that of Sz while using exclusively
the arithmetic real one.
This approach requires the use of an auxiliary real matrix to store factorized
the dynamic matrix D ().
Choice P = Im [(A - B)-1] - 1
·
This choice corresponds to the partly imaginary approach of the operator S.
Formally we have:
Im [(A - B)-1] - 1 = [(A - B) B-1 (A B) + 2 B B-1 B]
The matrix B is regular under the condition necessary and sufficient that M is it and
One obtains then:
Im [(A - B)-1] - 1
If M is singular, one can establish this equality by defining pseudo-opposite
of B
by
The calculation of Pz can then be carried out without assembling the matrix P explicitly and while using
exclusively the arithmetic real one.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
21/26
5.3.3 Cost of the Lanczos phase
The cost report corresponds to the allowance of additional vectors (3 realities and 1 complex or real)
and with the allowance of the dynamic matrices used for the calculation of the operators S and P.
The following table summarize these allowances
= 0
= IR *
= I I IR
= I IC
+
Approach
S
K IRnxn
D () IRnxn
D () ICnxn
D () ICnxn
part
P = B
-
-
-
-
real
P = Re (A - 1) - 1
-
-
K IRnxn
D () IRnxn
Approach
S
D () ICnxn
D () ICnxn
part
P = B
-
-
Imaginary
P = Im (A - 1) - 1
-
-
The cost in operation is divided into a fixed blow and a blow depend on the number of vector of Lanczos
to calculate.
The fixed cost corresponds to factorization LDLT of the additional matrices (carried out in
the arithmetic associated one) and is worth O (b2n) if B is the bandwidth common to the matrices M, C and K.
The calculation of a vector of Lanczos requires:
2 scalar products of vector of IR2n:
2
O (2n)
· 1 linear combination of 3 vectors of IR2n:
3
O (2n)
· the calculation of Sz:
·
3 products matrix-vector in IRn:
3
O (2bn)
1 linear combination of 3 vectors in ICn:
3
O (2n)
1 descent-increase in ICn:
O (2bn)
1 scalar product - vector in ICn:
O (N).
The calculation of Pz
·
For P = b:
3 products matrix-vector in IRn:
O (2bn)
For P = Re (A - 1) - 1:
5 products matrix-vector in IRn:
5
O (2bn)
2 combinations linear of 3 vectors of IRn:
6
O (N)
1 descent gone up in IRn:
O (2bn)
1 linear combinations of 5 vectors of IRn:
5
O (N)
For P = Im (A - 1) - 1:
6 products matrix-vector in IRn:
6
O (2bn)
2 linear combinations of 4 vectors of IRn:
8
O (N)
Broadly the cost of the Lanczos phase is out of O (b2n) + 10 m O (2bn).
5.4
Implementation in Aster
The matrices M and C are symmetrical semi-definite positive and the matrix K is symmetrical regular
indefinite. The pseudo scalar product retained corresponds to the extension of that proposed by Neuman and
Pipano [R5.01.01].
This algorithm is available in Aster by operator MODE_ITER_SIMULT.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
B
R
0
H
B
R
0
,
R
0
B
=
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
22/26
5.4.1 Parameters of the implementation
The spectral problem with real operator is parameterized by:
value of the spectral shift IC,
· the choice of Re (A - 1B)
or of Im (A - 1B)
·
The pseudo scalar product is then Re (A - 1) - 1 or Im (A - 1) - 1 and allowances and factorizations of
additional matrices are carried out has minimum in type and a number, this in agreement with the table
paragraph [§5.3.3].
5.4.2 Under space of approximation
The vector r0 IR2n generating under space of approximation breaks up into
where
RN; the choice selected consists in posing = 0 and drawing by chance the components from,
while imposing to him null components in ker Mr.
If the dimension of under space is not specified, it is calculated by the empirical formula:
m = 2 Min (Max (p+7), 2p, N)
where p is the number of couples of elements suitable to approach.
the dimension of the subspace of approximation are doubled because the couples of clean elements
complexes arise per combined pairs.
5.4.3 Strategy of reorthogonalisation
The use of arithmetic with finished precision, deteriorates the properties of orthogonality of the vectors of
Lanczos and with them the rate of convergence of the approximate couples.
The strategy of complete reorthogonalisation ensures the orthogonality of all the vectors of Lanczos, and
the algorithm then has a behavior close to that into arithmetic exact.
This strategy of reorthogonalisation forces to preserve the vectors P. qj j= 1,…, Mr.
The reorthogonalisation of the vectors is carried out by the process of Gram-Schmid modified cf.
[R5.01.01].
5.4.4 Implementation of the Lanczos phase
The selected implementation is that described in [R5.01.01] within the framework of the generalized problem.
El E is summarized by:
Inputs
:
·
the matrices P and S: i.e. the matrices M, C and K and factorized the LDLT of the matrices
-
dynamic additional.
m the number of vectors to be generated,
- r0 the vector generating under space of Krylov,
- precision of orthogonalization and the number maximum of authorized reorthogonalisation.
-
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
(
m
(
)
m
)
J
,
µ
X
N
S
X
J
S
=
O
Q
m
S
J
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
23/26
Outputs
:
·
vectors of Lanczos (q1, q2,…, qm),
- the diagonal (1, 2,…, m) and the on-diagonal (1, 2,…, m) of the matrix tridiagonale
-
Tm
the vector (1, 2,…, m) of the pseudo scalar products of the vectors of Lanczos.
-
Algorithm:
·
Generation of the first vector q1 and coefficients 1, 1 and 1
-
Loop generation of qj, J, J and J for J = 2,…, m
-
For J = 2,3,…, m to make
Calculation of the direction of qj
Standardization of qj, calculation of J and storage of Pqj
Réorthogonalisation so necessary compared to qj for I = 1,…, j-1
Reactualization of qj, J, J and J in the event of reorthogonalisation,
Calculation of J and J
5.4.5 Restoration of the approximations for the quadratic problem
Approximations
clean couples of the quadratic problem result from
clean couples
JmTm matrix by:
extraction of the “high” part of the vector.
·
choice
of, root of the flexible quadratic equation
with
who checks the condition of
·
coupling M [O N] Qm
=
M
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
24/26
6 Bibliography
[1]
[CUL.WIL.85] J.K.CULLUM & R.A. WILLOUGHBY, broad Lanczos algorithms for symetric
eigenvalue computations - VOL1 Theory, Birkhäuser, 1985.
[2]
[GOH.al.86.] I. GOHBERG, P. LANCASTER & L. RODMAN, Quadratic matrix polynomials
with has parameter, Advances in applied mathemetics 7, pp253-281, 1986.
[3]
[JEN.77] A. JENNINGS, Matrix computation for engineers and scientists, John Wisley & Sons,
1977.
[4]
[MEI.67] L. MEIROVITCH, Analytical methods in vibrations, The MacMillan Co. N.Y., 1967.
[5]
[MUL.56] D.E. MULLER, A method for solving algebraic equations using year automatic
computer, Maths. . Wash. Flight 10, p 208-215,1956.
[6]
[PAR. SAA.87] B.N. PARLETT & Y. SAAD, Complex shift and invert strategies for real
matrices, Linear will algebra and its N°88 applications/89, pp575-595, 1987.
[7]
[ROS. 84] Mr. ROSEAU, Vibrations of the mechanical systems, methods analytical and
applications, Masson, 1984.
[8]
[WIL.65] J.H. WILKINSON, The algebraic eigenvalue problem, Oxford University Press,
London 1965.
[9]
[R5.01.01] D. SELIGMANN - Algorithmes of resolution for the problem generalized with
eigenvalues. Note EDF HI-75/7815 - Documentation Aster [R5.01.01].
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
I
T
C
2
C
I
K
K
2
Re
I
I
*
T
m
D
K
,
M
I
=
m
±
2
I
I
±
1
I
I
C
I
=
I
I
+
+
1
I
-
E
I
I
=
-
I
T
=
2
2
R E
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
25/26
Appendix 1 Interprétation of the complex eigenvalues
In the case of a symmetrical damping and in internal absence of damping, relations
orthogonalities and owing to the fact that the clean elements appear per combined pairs, one has the relations
following:
If one notes
, one can then define:
and one can write the eigenvalue complexes in the following form:
for which one can give a physical interpretation of the eigenvalue
The imaginary part represents the oscillatory part of the solution
is the pulsation of the i-ème mode
the real term represents the dissipative character of the system
is the damping of the i-ème mode,
is the reduced damping of the i-éme mode.
Physical interpretation of the clean vectors:
The physical significance of the existence of a clean vector complexes, lies in the fact that if
·
structure vibrates on a clean mode, its various degrees of freedom do not vibrate with the same one
phase ones compared to the others.
The modal bellies and nodes do not correspond of the stationary points, but
·
move during the movement.
Note:
one finds the traditional formulation of the deadened systems with 1 degree of freedom
·
are real and are many quantities intrinsic with a mode (modal quantities) and
·
dependant on the standardization of the mode.
We point out that the modes of the quadratic problem do not diagonalisent the matrices M, K and
C.
Notice on the real term of the eigenvalue:
If the real part of the eigenvalue is negative, then the clean mode is a movement
·
deadened periodical of pulsation
If the part R
éelle of the eigenvalue is positive, then the clean mode is a movement
·
periodical of amplitude increasing and thus unstable.
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
2
M
0
1
C
C
X
M
0
0
N
M
1
K
K
X
1
0
M
C
0
N
.
- K
+
K
X
C
K
2n
- K
- 1
M
-
+
0
- 1
C
0
=
X
=
K
0
=
0
0
Code_Aster ®
Version
2.0
Titrate:
Algorithms of resolution for the quadratic problem
Date:
19/06/92
Author (S):
D. SELIGMANN, R. MICHEL
Key:
R5.01.02-A
Page:
26/26
Appendix 2 Réductions linear
form
problem
quadratic
generalized
(1)
generalized
(2)
standard
(1)
standard
(2)
Note: to obtain the forms standards it is necessary to suppose:
M and K regular for the form (1)
· M regular for the form (2)
·
Handbook of Référence
R5.01 booklet: Modal analysis
HI-75-7816/A
Outline document