Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
1/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Organization (S):
EDF/MTI/MN
Manual of Reference
R5.01 booklet: Modal analysis
Document: R5.01.01
Algorithm of resolution for the problem
generalized
Summary
In this document, we present the algorithms of resolution for the generalized modal problems which
are established in Code_Aster via the operators
MODE_ITER_INV
and
MODE_ITER_SIMULT
:
·
method of the iterations opposite,
·
method of Lanczos,
·
the method WILL GO (known as of “Sorensen”),
·
method of Bathe & Wilson.
One gives to the reader the properties and the limitations, theoretical and practical, of the modal methods approached
while connecting these considerations, which can sometimes appear a little “éthérées”, to a precise parameter setting of
operators. For each method, one recapitulates in the form of tables the aforementioned parameter setting with his values by
defect and of the references to the paragraphs of the document.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
2/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Contents
1 Introduction - Description of the document ................................................................................................. 4
2 Context ................................................................................................................................................. 6
2.1 Problems ................................................................................................................................. 6
2.2 Taking into account of the conditions limit ........................................................................................... 7
2.3 Properties of the matrices ................................................................................................................... 8
2.4 Properties of the clean modes ........................................................................................................ 9
2.5 Estimate of the real spectrum ............................................................................................................. 10
2.6 Establishment of the test of Sturm ....................................................................................................... 14
2.7 Spectral transformation ............................................................................................................... 16
2.8 Modal calculation .................................................................................................................................. 17
2.9 Establishment in Code_Aster .................................................................................................. 19
3 Method of the powers opposite (
MODE_ITER_INV
) ....................................................................... 22
3.1 Introduction .................................................................................................................................... 22
3.2 Localization and separation of the eigenvalues ............................................................................. 23
3.2.1 Method of bisection ......................................................................................................... 23
3.2.2 Method of the secant ......................................................................................................... 24
3.3 Method of the powers opposite ................................................................................................ 25
3.3.1 Principle ................................................................................................................................ 25
3.3.2 Method of iteration of the quotient of Rayleigh ........................................................................ 27
3.3.3 Establishment in Code_Aster ......................................................................................... 28
3.3.4 Display in the file message ....................................................................................... 29
3.3.5 Summary of the parameter setting ............................................................................................... 30
4 Method of subspace (
MODE_ITER_SIMULT
) .............................................................................. 31
4.1 Introduction .................................................................................................................................... 31
4.2 Analyze of Rayleigh-Ritz ............................................................................................................... 31
4.3 Choice of the space of projection ..................................................................................................... 33
4.4 Choice of the spectral shift ........................................................................................................... 35
5 Method of Lanczos (
METHOD = “TRI_DIAG”
) ............................................................................. 36
5.1 Introduction .................................................................................................................................... 36
5.2 Theoretical algorithm of Lanczos .................................................................................................. 36
5.2.1 Principle ................................................................................................................................ 36
5.2.2 Estimates of errors and convergences ........................................................................... 38
5.3 Algorithm of Lanczos practices .................................................................................................... 40
5.3.1 Problem of orthogonality ..................................................................................................... 40
5.3.2 Capture multiplicities ...................................................................................................... 41
5.3.3 Phenomenon of Lanczos ....................................................................................................... 41
5.4 Complementary processing ....................................................................................................... 42
5.4.1 Detection of spaces invariants ............................................................................................. 42
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
3/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5.4.2 Strategies of restartings ................................................................................................. 43
5.5 Establishment in Code_Aster ................................................................................................... 44
5.5.1 Alternative of Newmann & Pipano ........................................................................................... 44
5.5.2 Parameter setting ......................................................................................................................... 46
5.5.3 Warning on the quality of the modes ............................................................................... 47
5.5.4 Perimeter of use ............................................................................................................ 47
5.5.5 Display in the file message ........................................................................................ 48
5.5.6 Summary of the parameter setting ................................................................................................ 49
6 Algorithm WILL GO (
METHOD = “SORENSEN”
) ....................................................................................... 50
6.1 Introduction .................................................................................................................................... 50
6.2 Algorithm of Arnoldi ....................................................................................................................... 50
6.2.1 Principle ................................................................................................................................. 50
6.2.2 Estimates of errors and convergence ............................................................................. 52
6.3 The stakes ...................................................................................................................................... 53
6.4 Algorithm “Implicit Restarted Arnoldi” (WILL GO) .................................................................................. 54
6.5 Establishment in Code_Aster ................................................................................................... 56
6.5.1 ARPACK ............................................................................................................................... 56
6.5.2 Adaptations of the algorithm of Sorensen ............................................................................. 56
6.5.3 Parameter setting ......................................................................................................................... 57
6.5.4 Display in the file message ........................................................................................ 58
6.5.5 Summary of the parameter setting ................................................................................................ 59
7 Method of Bathe and Wilson (
METHOD = `JACOBI
) ........................................................................ 60
7.1 Principle .......................................................................................................................................... 60
7.2 Tests of convergence .................................................................................................................... 60
7.3 Establishment in Code_Aster ................................................................................................... 60
7.3.1 Dimension of the subspace .................................................................................................. 60
7.3.2 Choice of the initial vectors ................................................................................................... 61
7.3.3 Parameters in Code_Aster ........................................................................................... 62
8 Conclusion - Synthesis .......................................................................................................................... 63
9 Bibliography ......................................................................................................................................... 65
Appendix 1 General information on algorithm QR .............................................................................................. 67
Appendix 2 Orthogonalization of Gram-Schmidt ....................................................................................... 72
Appendix 3 Method of Jacobi .................................................................................................................. 75
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
4/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
1
Introduction - Description of the document
A majority of study concerning the dynamic behavior of structures are carried out in
carrying out a transitory analysis on modal basis. To exhume these modes of vibrations, one
string of algorithms were developed since about fifty years. In order to face
continual increase in the size of the problems and with the degradation of conditionings of
the discretized operators, only most effective and most robust, in practice, were built-in
in the two modal operators of Code_Aster.
The optimal perimeters of use of these operators can be dissociated. When it is about
to determine some eigenvalues (typically a half-dozen) or to refine some
estimates, the operator
MODE_ITER_INV
is completely indicated. It gathers algorithms
heuristic and those of type powers (cf [§3]).
On the other hand, to capture a significant part of the spectrum, one A resorts to
MODE_ITER_SIMULT
.
This last federates the methods known as of “subspace” (Lanczos [§4], [§5], IRAM [§6], Bathe & Wilson
[§7]) which projects the operator of work in order to obtain an approximated spectrum of more reduced size (treated
then by a total method of type
QR
or Jacobi).
Until now, these algorithms stumbled regularly on the same shelves: detection
correct of multiple modes, modes of rigid body and a general way, processing of
packed spectrum. All this led to the appearance of “phantom” modes sometimes badly easy to detect
(modes corresponding to multiplicities missed and being able to generate correct residues with the direction
of Aster, and this, more especially as the criteria of the post-modal checks were sometimes permissive
(residue in 10
2
instead of the 10
6
current), decontaminated even insufficient (test of Sturm limited to the values
clean positive) which causes distortions of results downstream from calculation, during projections on
base modal.
To be been free from the recurring problems to this type of approach, one thus proposed to enrich
MODE_ITER_SIMULT
(starting from V5) from the algorithm WILL GO (“Implicit Restarted Arnoldi” [§6]). This
alternative of Arnoldi, initiated by D.C. Sorensen in 1992, makes real great strides for the resolution of large
modal systems on parallel supercomputers.
It tries to bring an elegant remedy for the recurring numerical problems raised by the others
approaches. In short, IRAM gets a better total robustness In V5, it made it possible to balance
all of software faults related to the modal problems generalized (in V5, it allowed
to balance all software faults related to the modal problems generalized) while improving
complexities calculation and memory for a fixed precision, and, it allows a real control on
quality of the modes via a suitable parameter setting. Its use (method by defect in
MODE_ITER_SIMULT
) is thus with advising in all the cases of figures.
This document is articulated around the following parts:
·
initially, one recalls the context of modal calculation in Code_Aster to
through matrices and limiting conditions used, particular properties of modes
clean exhumed and of the difficulties of estimate of the spectrum. One recapitulates also what it is necessary
to know with minima about the spectral transformations, the families of modal algorithms and
on the broad outline of the course of a modal calculation in the code,
·
in the second time, one describes the method of the powers opposite and his phases
preliminaries of localizations of eigenvalues (method of bisection and the secant)
installations in
MODE_ITER_INV
,
·
the third part treats the general framework of the methods of subspace installation in
MODE_ITER_SIMULT
. One details in particular the implications of the choice of the shift and the type of
subspace of projection in the parameter setting of the operator,
·
the three following chapters show, in the order, the three methods of this operator:
Lanczos, IRAM and Bathe & Wilson,
·
in the appendices one approaches more in details of the processes “of second level” to which make
call three preceding methods. They are the algorithms
QR
and Jacobi who are known as
general practitioners because they allow to capture all the spectrum of an operator. The first,
the algorithm
QR
, is fundamental because it intervenes in the majority of the methods. It is
the algorithm of reference which offers a very great robustness but calculation complexities and
memory prohibitory. One approaches also the various algorithms of orthonormalisation put in
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
5/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
place in the code. Indeed, one will cease hammering, throughout this document, only
their quality and their robustness are crucial for the good unfolding of the algorithms.
This last version of the document was almost entirely rewritten while taking as a starting point the the indices
written precedents, respectively by D. Seligmann [R5.01.01] index A and B. Quinnez [R5.01.01]
index B, this in order to try to approach the code while clarifying more in details them
characteristics of the methods and the subjacent numerical phenomena. A particular effort was
brought to put in prospect the choices led in Code_Aster compared to search,
passed and current, like clarifying the general philosophy of a modal calculation.
One gives to the reader the properties and the limitations, theoretical and practical, of the modal methods
approached while connecting these considerations, which can sometimes appear a little “éthérées”, to one
precise parameter setting of the operators. For each method, one recapitulates in the form of tables the aforementioned
parameter setting with its default values and of the references to the paragraphs of the document.
At the time of the first passages, one strongly engages the user has to modify only the parameters
main noted in fat in these tables and with reading the traces of the file message. Others,
concerning more the mysteries of the algorithms, were initialized empirically with values
standards.
One tried constantly to bind different the items approached and to limit to the bare minimum the recourse
with long mathematical demonstrations. In any event, the many references which enamel
the text must make it possible to seek precise information.
The object of this document is not to detail all the aspects approached, of the complete works having
already fulfilled this mission. One will quote in particular F. Chatelin [bib3], G.H. Golub [bib6], P. Lascaux
[bib11], B.N. Parlett [bib18] and Y. Saad [bib32], and one recommend the synthesis particularly
brought up to date and exhaustive made by J.L. Vaudescal [bib23].
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
6/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
2 Context
2.1 Problems
We consider the problem generalized with the eigenvalues:
To find (
,
U
) such as
With U
B U
=
,
U
0, where
With
and
B
are symmetrical matrices with coefficients
realities. This type of problem corresponds, in mechanics, in particular with:
·
The study of the free vibrations of a not deadened and nonrevolving structure. For this
structure, one seeks the smallest eigenvalues or those which are in one
interval given to know if an exiting force can create a resonance. In this case,
stamp
With
is the matrix of rigidity, noted
K
, (possibly increased matrix of
noted geometrical rigidity
K
G
, if the structure is prestressed) and
B
is the matrix of mass
or of noted inertia
M
. The eigenvalues obtained are the squares of the associated pulsations
at the sought frequencies.
The system to be solved can be written:
(
)
K K
U
M U
+
=
G
2
where
=
2 F
is the pulsation,
F
the Eigen frequency and
U
the vector of associated clean displacement.
·
The search for linear mode of buckling. Within the framework of the linearized theory, in
supposing a priori that the phenomena of stability are suitably described by
system of equations obtained by supposing the linear dependence of displacement by report/ratio
at the level of critical load, the search of the mode of buckling
U
associated this level of
critical load
, brings back itself to a problem generalized to the eigenvalues of the form:
(
)
K
K
U 0
+
=
G
with
K
stamp rigidity and
K
G
stamp geometrical rigidity.
Note:
·
this type of generalized modal problem is treated in Code_Aster by two operators:
MODE_ITER_INV
and
MODE_ITER_SIMULT
. Each one having its perimeter of application, its
functionalities and its limitations,
·
in V5, the user can specify the class of membership of his calculation by initializing it
key word
TYPE_RESU
with '
DYNAMICS
“(default value) or with”
MODE_FLAMB
'. Display of
results will then be formatted by taking account of this specificity. In the first case one
will speak about frequencies whereas in the second, one will speak about critical load,
·
in the presence of depreciation and gyroscopic effects, the study of dynamic stability
of a structure leads to the resolution of a modal problem of a nature higher, known as quadratic:
(
)
K
C
M U 0
+
-
=
I
2
. It is solved by the two modal operators and is the object
of a specific note [R05.01.02].
Now that links between the mechanics of the structures and the resolution of modal problems
generalized were recalled, we will be interested in the processing of the conditions limit in
code and with their angles of attack on the matrices of mass and rigidity.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
7/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
2.2
Taking into account of the limiting conditions
There are two ways, during the construction of the matrices of rigidity and mass, to take into account
the boundary conditions (this description in dynamic term of problem is extrapolated easily with
buckling):
·
double dualisation, by using ddls of Lagrange [R3.03.01], makes it possible to check
C U =
0 (CL for Linear Limiting Condition),
with
C
real matrix of size
p X N
(
K
and
M
are of command
N
). Matrices of rigidity and of
mass dualized have the form then
~
~
K
K
C
C
C
Id
Id
C
Id
Id
M
M 0 0
0
0 0
0
0 0
=
-
-
=
T
T
with
and
strictly positive realities which will be initialized to 1 without that involving a loss
of general information.
The dimension of the problem was increased by
2p
, because with
N
ddls known as “physics”, one has
added ddls of Lagrange. There are two ddls of Lagrange by linear relation affected to
p
limiting conditions.
·
The setting has zero of
p
lines and columns of the matrices of rigidity and mass. This is not
valid that for blockings of ddls. One cannot take into linear account of relation
and one will speak about kinematic blocking (CLB for Limiting Condition of Blocking). Matrices
of rigidity and mass become:
~
~
K
K
0
0
Id
M
M 0
0
0
=
=
The dimension of the problem remains unchanged but it is however necessary to withdraw the participations of
ddls locked with the components of the initial matrices (
K
is obtained from
K
in
eliminating the lines and the columns from the ddls which are locked; idem for
M
).
When limiting conditions are imposed, the number of eigenvalues (with all their multiplicities)
really implied in physics (with the flat of modeling close (mesh, frequency of
rupture…) phenomenon is thus lower than the size
N
transformed problem:
·
N
N
p
with p
p
ddl active
-
= -
=
3
2
2
'
'
(double dualisation),
·
N
N p
ddl active
-
= -
(kinematic blocking).
Framed below the display dedicated to these parameters in the file message shows.
------------------------------------------------------------------------
THE NUMBER OF DDL
TOTAL EAST: 220
N
LAGRANGE EAST: 58
p'
= 2
p
A NUMBER OF ACTIVE DDL EAST: 133
N
ddl_actifs
-----------------------------------------
Example 1: Management of the ddls
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
8/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
In addition, in the modal calculation algorithms, one must ensure oneself of the membership of the solutions with
acceptable space. One brings back oneself there via auxiliary processing. Thus when blockings are used
kinematics (CLB), it is necessary in the various algorithms and for each iteration, to use a vector “of
positioning "
U
bloq
, defined by:
·
if
ième ddl is not locked
()
U
bloq
I
=
1
,
·
if not
()
U
bloq
I
=
0
,
()
()
() (
)
U
U
U
U
1
0
1
1
I
I
I
I
N
bloq
=
=
.
…
If one uses the method of double dualisation, one needs a vector of positioning of the ddls for
lagrange
U
lagr
defined like
U
bloq
. It is only used at the time of the choice of the random initial vector.
So that this vector
U
0
check the conditions limit (CL) one operates in the following way:
()
()
() (
)
U
U
U
K U
U
U
1
0
2
1
2
1
I
I
I
I
N
lagr
=
=
=
.
…
~
Thereafter, to simplify the notations, we will not make the distinction between the initial matrices and theirs
hanging dualized (noted with a tilde) that if necessary. Very often, they will be indicated by
With
and
B
in order to approach the usual modal notation without being attached to such or such class of
problems.
2.3
Properties of the matrices
As we wrote previously the matrices considered are symmetrical and with
real coefficients. According to the cases of figure indexed in the table below, they can be
defined positive (noted
> 0
), semi-definite positive (
0
), indefinite (
0
0
or
) even singular
(S).
Structure
free
Lagranges *
Buckling
Fluid
structure
K
0
and S
< 0 or > 0
> 0
0
and S
M
(resp.
K
G
)
> 0
0
and S
0
0
or
> 0
Table 2.3-a: Properties of the matrices of the generalized problem
* This column relates to of course the properties of the matrices dualized made up to leave
initial matrices [R3.03.01]
Columns of this mutually being excluded table, in practice, a problem of using buckling
Lagranges doubles to modelize some of its limiting conditions, sees its dualized matrices
(
~
~
K
K
and
G
) to become potentially indefinite.
Note:
This range of properties must be taken into account at the time of the choice of the couple operator of work -
(pseudo) produced scalar. This framework can thus reinforce, with effectiveness and transparency,
robustness and the perimeter of the modal calculation algorithm in all the encountered cases of figure
by Code_Aster.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
9/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
The following paragraphs go enable us to measure the angle of attack of these properties on the spectrum
of generalized problem.
2.4
Properties of the clean modes
Let us recall first of all that if the matrix of the standard modal problem
With U
U
=
is real
symmetrical, then its own elements are real; The clean elements of a matrix are its
clean values and its vectors. In addition,
With
being normal, its own vectors are orthogonal.
In the case of the generalized problem
With U
B U
=
, this condition is not sufficient. Thus,
let us consider the following generalized problem:
1 1
1 0
1
0
0
1
1
2
1
2
=
-
U
U
U
U
its own modes are
(
)
±
±
±
±
=
±
=
+
-
1
2 1
3
1
1
1
2
I
and U
.
If one adds the assumption “one of the matrices
With
or
B
is definite positive ", then the problem
generalized has its real solutions. There is even the characterization (sufficient condition) more precise
following.
Theorem 1
Are
With
and
B
two real symmetrical matrices. If there exists
and
such as
With
B
+
that is to say definite positive, then the generalized problem has its real own elements.
Proof:
This result is obtained immediately by multiplying the problem generalized by
and by carrying out one
spectral shift
. The problem then is obtained
(
) (
)
With
B U
B U
+
=
+
. Like
(
)
With
B
+
is definite positive, it admits a single decomposition of Cholesky in the form
DC
T
with
C
regular matrix.
The problem is written then
C B C
Z
Z
-
-
=
1
T
µ
with
Z C U
=
T
and
µ
=
+
1
, which allows
to conclude, because the matrix
C B C
-
-
1
T
is symmetrical.
Notice
This characterization is not necessary, thus the generalized problem associated the matrices
(
)
With
=
- -
diag 1 2 1
,
and
(
)
B
=
-
diag
2 11
,
a spectrum admits real all while not answering
condition of definite positivity.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
10/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Proposal 2
If matrices
With
and
B
are real and symmetrical, the clean vectors of the generalized problem are
With
and
B
- orthogonal, which means that they check the relations
U B U
U With U
iT
J
ij
J
iT
J
J ij
J
has
has
=
=
where
has
J
is a scalar depend on the standard of the J
ième
clean vector,
ij
is the symbol of Kronecker
and
U
J
is the clean vector associated the eigenvalue
J
.
Proof:
Immediate for distinct eigenvalues, by writing them
With
and
B
- scalar product between two
couples
(I, J)
and
(J, I)
, then by using the symmetry of the matrices (cf [bib9]).
Note:
·
it is shown that them
With
and
B
- orthogonalities of the clean vectors are a consequence of
hermiticity of the matrices. They are clearly a generalization of the properties of the problem
square standard (even normal),
·
orthogonality compared to the matrices does not mean especially that the clean vectors are
orthogonal for the conventional euclidian norm. The aforementioned can be only the fruit of
particular symmetries (cf TP n°1 [bib25]),
·
this property simplifies calculations of modal recombinations (
DYNA_TRAN_MODAL
[R5.06.01]), when one handles matrices of generalized rigidity and mass which are
diagonals. Quantities
K
J
=
J
has
J
and
m
J
=
has
J
are called, respectively, modal rigidity
and modal J masses
ième
mode.
Knowing that the modes are real, we now will worry we about their estimate.
2.5
Estimate of the real spectrum
Because of membership of the spectrum to the real axis, the problems of counting of eigenvalues go
to be largely simplified. One does not have to treat régionnement complex plan and one can
to be based on the corollary of the law of inertia of following Sylvester.
Corollary 3
Are
With
and
B
two symmetrical real matrices,
B
being of more definite positive. The number of
eigenvalues, strictly lower than
, of the generalized problem
With U
B U
=
is then equal to
a number of strictly negative diagonal coefficients of the matrix
D
such as
(
)
With
B
LDL
-
=
T
.
Proof:
Cf paragraph n°1 of the article of Y. Haugazeau [bib7].
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
11/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Remarks
·
this corollary extends to the square matrices and results from the properties of the series Sturm
and of the principle of inclusion (one notes
()
()
()
()
(
)
p
dét
N
N
N
=
-
With
B
the characteristic polynomial
generalized problem
() ()
() () ()
With
U
B
U
N
N
N
N
N
=
obtained by removing them
N
last
lines and columns of the matrices
With
and
B
. The continuation of polynomials
()
()
p
N
N
constitute a continuation of
Sturm because roots of I
ème
polynomial frame those of (i+1)
ème
. This property
of interlacing of the spectrum of symmetrical real submatrices “principle is called
of inclusion ") [bib7], [bib9],
·
the possible multiple eigenvalues are taken into account their multiplicity,
·
thereafter, one will call modal position of
and it will be noted
pm
(
), this number of
strictly negative diagonal coefficients.
This corollary thus makes it possible to easily determine the number of eigenvalues contained in one
interval
[]
µ
,
and the modal position of these eigenvalues in the spectrum
.
It is enough to carry out two
decompositions
LDL
T
, that of
(
)
With
B
-
and that of
(
)
With
B
- µ
and to enter the difference
number of strictly negative terms between the two diagonal matrices. In the jargon of the code,
one indicates (improperly besides!) this test under the term of “test of Sturm”.
It must however be extended to the quite particular shapes of the matrices met in
Code_Aster, and in particular, it is necessary to be able to take into account buckling with or without Lagrange.
With this intention, one widened the criterion with a matrix
B
unspecified and one generalized it while taking in
count the dualized matrices.
Let us recall that one defines usually the signature of a matrix (and that of the quadratic form
associated) like the triplet of natural entireties
(R, S, T)
where
R
indicate the number of eigenvalues
> 0
,
S
a number of null eigenvalues and
T
those
< 0
. The latter entirety is thus what one notes
pm
(
)
in our case of figure.
Property 4
Are
With
and
B
two symmetrical real matrices (of command
N
) related to the generalized modal problem
(S):
With U
B U
=
. Let us note
~
With
and
~
B
, their matrices associated resulting from the double dualisation with
p
Lagranges allowing to check
C U =
0 (CL), with
C
real matrix of size
p X N
.
Then,
, the signature of
(
)
~
~
With
B
-
is written, while noting
[]
card
,
B has
the number of values
clean of the generalized problem included in the interval
[]
B has
,
:
if
B
is indefinite and
With
is definite positive
then
{}
card
0
0
=
and
if
<
0
:
]
[
[
[
{}
()
]]
R
S
p
T
p
=
-
+
+
= +
=
=
+
card
,
card
,
card
pm
card
,
0
0
éq 2.5-1
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
12/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
if
=
0
:
]
[
()
R
S
p
T
p
=
- +
=
=
=
card
,
pm
éq 2.5-2
if
>
0
:
]
]
]
[
{}
()
[[
R
S
p
T
p
=
- +
+
=
+
=
=
+
card
,
card
,
card
pm
card
,
0
0
éq 2.5-3
if
B
is definite positive and if
With
is semi-definite positive
then
]
[
card
,
- =
0
0
and
if
=
0
:
]
[
{}
()
R
S
p
T
p
=
+
= +
=
=
card
,
card
pm
0
0
éq 2.5-4
if
>
0
:
]
[
{}
()
[[
R
S
p
T
p
=
+
= +
=
=
+
card
,
card
pm
card
,
0
éq 2.5-5
Proof:
To impose boundary conditions linear, a technique of double dualisation is used
[R3.03.01] which leads to the transformed generalized system
(
)
()
~
~ ~
~
With
B U
With
C
C
C
Id
Id
C
Id
Id
B 0 0
0 0 0
0 0 0
U
v
W
0
-
=
-
-
-
=
T
T
S
One will note thereafter
v
ij
the null vector column of size
p (i= 1. .p; j= 2,3)
except with the index
I
, for
which it is worth 1, and
0
N
, the null vector column of size
N
. Now let us consider the three families of
vectors following:
·
N p
-
vectors
V
U
v
v
I
I
I
I
1
1
1
1
=
who are the clean vectors of the system
(~)
S
. Clearly,
according to proposal 2, they are
~
With
and
~
B
- orthogonal (according to the properties of
With
and
B
) and
one notes their eigenvalues
I
.
·
p
independent vectors
V
0
v
v
I
N
I
I
2
2
2
=
,
·
p
independent vectors
V
0
v
v
I
N
I
I
3
3
3
=
-
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
13/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
These
N p
+
vectors form a base B of space
E
=
=
=
~
/
U
U
v
W
Cu 0
solutions
acceptable of the dualized problem (as a free family of a space of finished size).
Let us consider, for a real number
given, the quadratic form associated the matrix
(
)
~
~
With
B
-
()
(
)
~
~
~ ~, ~, ~
U
With
B U U
U
=
-
E
While breaking up on the basis B generated by the preceding vectors
~
,
,
,
U
V
V
V
=
+
+
= -
=
=
has
has
has
I
I
I
N p
I
I
I
p
I
I
I
p
1
1
1
2
2
1
3
3
1
,
one obtains
()
()
(
)
(
)
()
()
~
,
,
U
U B U
=
-
+
-
= -
=
has
has
I
I
I
I
I
N p
I
I
p
T
1 2
1
1
1
3 2
1
4
(one used in particular the properties of orthogonality of the families of vectors
V
ij
and the relation
(
) (
)
U
C v
C U
v
I
T I
I
I
1
2
1
2
0
,
,
=
=
). From where, while noting
L
I
the linear form which associates
~
U
its
I
ème
coordinated in the base B
()
()
(
)
(
)
()
()
()
~
L ~
L
~
L
~
,
,
,
U
U
U B U
U
U
=
-
+
+
-
= -
- +
=
+
=
I
I
I
I
I
N p
N p I
I
p
N I
I
p
T
2
1
1
1
2
1
2
1
0
4
Clearly them
()
L
,
I I N p
= +
1
are linearly independent from where an immediate reading of the terms of
the signature according to the sign of the factors. Moreover, it is noticed that this decomposition comprises
obligatorily
p
zeros and
p
minus signs. The relations of the property result all then
naturally, by using the principle of inertia of Sylvester (who ensures us that this decomposition
is invariant), the relations of
With
and
B
- orthogonality of property 2, properties of the table
[Table 2.3-a] and by noticing that:
~ ~ ~
~ ~ ~
U B U
U B U
U With U
U With U
iT
J
I
T
J
iT
J
J
I
T
J
J
J
=
=
=
for
0
.
The restriction
0
case of figure (2) comes owing to the fact that if
B
then the spectrum is definite positive
problem (S) is positive and thus the shift with a strictly negative value does not have any interest (one
find within the framework of application of corollary 3).
According to the table [Table 2.3-a] this property applies to the matrices handled by
Code_Aster and one can build the following corollary.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
14/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Corollary 5 (theoretical wide Sturm)
In the matric configurations of Code_Aster, the number of clean modes of the problem
generalized (S) whose clean vector checks the linear conditions limiting (CL) and of which the value
clean is contained in the interval
] [
µ
,
is
If
µ
0
] [
()
()
card
,
pm
pm
µ
µ
=
-
,
if not
] [
()
()
card
,
pm
pm
µ
µ
=
+
-
2 p
.
If no dualisation of Lagranges is required, one poses
p
= 0.
Proof:
It is immediate by combining the results of the table with those of property 4, and, while noticing
that the generalized modal configurations of Code_Aster can be only of two types:
·
buckling
With
is definite positive and
B
is indefinite, the spectrum can be negative and one uses
relations of the type (1),
·
dynamics
With
is semi-definite positive and
B
is definite positive, the spectrum is positive and one
can use the relations of the type indifferently (1) or (2).
Of course, if no dualisation of Lagranges is required, the same reasoning applies in
posing
p
= 0,
~
, ~
~
With
WITH B
B
U U
=
=
=
and
.
To finish, the principle of inertia of Sylvester ensures us that the signature of factorization
(
)
~
~
~ ~ ~
With
B
LDL
-
=
T
is identical to that of the shiftée matrix. In spite of this transformation, them
modal positions thus remain many perennial information.
However this corollary is tributary of arithmetic exact, in practice it is necessary to adapt it and control
its application.
2.6
Establishment of the test of Sturm
This test is confronted into arithmetic finished with two concomitant problems:
·
the factorization of the shiftée matrix when
is very close to an eigenvalue,
·
the calculation of the strictly negative terms of
~
D
to evaluate the modal position
()
pm
.
The first point requires, au préalable, the determination of a criterion of membership of the shift with
spectrum of the problem. In Code_Aster that Ci is founded on the loss of decimals at the time of
factorization of the matrix
(
)
~
~
~ ~ ~
With
B
LDL
-
=
T
. More precisely,
is regarded as being one
eigenvalue, if during factorization one loses more
NPREC_SOLVEUR
decimals (in any rigor it
is only one test of bad numerical conditioning). It is then necessary to modify the value of
in
shifting this shift of
PREC_SHIFT
% following the algorithm:
(
)
For I = 1,
If loss of more than
decimals
then
1
,
If not
exit;
Fine loops.
NMAX_ ITER_ SHIFT
NPREC_ SOLVEUR
PREC_ SHIFT
+
Algorithm 1: Shift of the shift
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
15/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
If at the end of
NMAX_ITER_SHIFT
attempts, the matrix is still not numerically invertible,
one continues all the same calculation with the last shiftée value. Framed below shows
trace of such a shift in the file message.
ONE MODIFIES THE VALUE OF SHIFT OF: 5.00000E+00POURCENT
THE VALUE OF SHIFT BECOMES: 5.96840E+04
ONE MODIFIES THE VALUE OF SHIFT OF: 5.00000E+00POURCENT
THE VALUE OF SHIFT BECOMES: 5.66998E+04
VALEUR_MIN IN FREQUENCY EAST: 2.00000E-01
VALEUR_MAX IN FREQUENCY EAST: 5.78813E+01
THE VALUE OF SHIFT IN FREQUENCY EAST:
3.78975E+01
Example 2: Shift of the shift
In algorithm 1 if
= -
F
F
corig
corig
then
. This parameter
F
corig
corresponds to a value
threshold in lower part of which it is considered that one has a numerically null eigenvalue (in
conventional mechanics, that corresponds to a mode of rigid body, cf [§5.5.4]). The imposition
= -
F
corig
thus allows to dissociate these modes of the remainder of the spectrum and to avoid instabilities
numerical during the test of Sturm. This threshold is skeletal with the key word
SEUIL_FREQ
.
Note:
·
the preceding key words are accessible starting from the key word factor
CALC_FREQ
from both
modal operators,
·
Y. Haugazeau [bib7] proposes more generally to deal with these problems of pivots
numerically “small” by building a matrix, unitairement similar (i.e.
similar via a unit matrix of passage (orthogonal in reality) to the shiftée matrix, of which
factorization would have less instability.
By taking into account these elements and knowing that, numerically, the calculation of the pivots
strictly negative of the diagonal matrix includes in fact also the elements (theoretically) null of
the signature, one can rewrite the preceding corollary.
Corollary 5bis (numerical wide Sturm)
According to the assumptions of corollary 5, there is the numerical criterion of accountancy of the modes according to:
If
µ >
0
[]
()
()
card
,
pm
pm
µ
µ
=
-
,
if
µ <
0
[]
()
()
card
,
pm
pm
µ
µ
=
+
-
4 p
Proof (heuristic):
The fact is applied that numerically the operator of factorization provides the “modal position
numerical "
()
pm
= +
S T
with property 4 and corollary 5. In addition, establishment of the criterion
neither the nullity of the product allows, nor the estimate of
{}
card
.
Note:
This criterion is used in preprocessings in
MODE_ITER_INV
(options '
ADJUST
“or”
SEPARATED
') and
in
MODE_ITER_SIMULT
(option '
BANDAGE
'), in postprocessings to check the validity of the number
modes (
MODE_ITER_SIMULT
) and in auxiliary controls (
IMPR_STURM
).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
16/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Now that we are able to enter the spectrum of the generalized problem, it remains with
to determine! The generic algorithms being intended for the standard problems, it is necessary to transform
our initial problem.
2.7 Transformation
spectral
These techniques make it possible to answer objective triple:
·
to exhume a standard modal problem,
·
to direct the search of the spectrum,
·
to separate the eigenvalues.
The spectral calculation algorithms converging of as much better than the spectrum (of work) which they treat
is separated, these techniques can be regarded as prepacking of the problem of
departure. They make it possible to return the separation of certain modes much more important than
those of other modes, and to improve their convergence thus.
Most widespread of these transformations is the technique known as of “shift and invert” which consists with
to work with the operator
With
such as:
(
)
With U
B U
With
B
B
With
U
U
=
-
= -
-
µ
1
1
! “
# #
$
# #
!“$
Appear 2.7-a: Effect of the “shift and invert” on the separation of the eigenvalues.
The figure [Figure 2.7-a] shows that this separation and this orientation of the spectrum of work in
µ
are
had with the particular properties of the hyperbolic function. In addition, it is observed that only them
eigenvalues are affected by the transformation. At the end of the modal process it is thus enough to
to pass by again in the plan of
by a suitable change of variable.
Note:
·
the variable
is usually indicated by the term of “shift” or spectral shift,
·
the matrix of work
With
must of course be invertible, that can become one besides of
motivations of this shift (cf [§5.5.1]).
µ
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
17/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
For memory, let us note that with a shift complexes several scenarios are possible:
·
to work completely into arithmetic complex,
·
into arithmetic real,
·
mixer two steps by isolating the real and imaginary contributions from
With
()
~
Re
~
With
With
µ
=
=
- + -
1
2
1
1
,
()
~
~
Im
~
~
With
With
µ
=
=
- - -
1
2
1
1
I
.
Each one of its approaches has its advantages and its disadvantages. For the quadratic problems,
it is for the moment the second solution which was adopted in the code.
Remarks
·
this choice of the operator of work is indissociable among that of (pseudo) - scalar product. It
allows to direct itself towards such or such algorithm and can thus influence the robustness of calculation,
·
another class of spectral transformation, with double shifts, allows to select them
eigenvalues located on the right vertical axis. It is about the rational transformation of
Cayley:
(
) (
)
With
B
With
B
With
U
U
-
-
= --
-
µ
1
1
2
2
1
!
“
# # # #
$
# # # #
!“$
.
We now will enumerate the various families of methods allowing to solve it
standard modal problem.
2.8 Calculation
modal
The methods of calculation modal can gather in (at least) four families:
·
Algorithms of the type
QR
(cf [Appendix 1]) which was presented per H. Rutishauser (1958)
and formalized jointly by J.C. Francis and V.N. Kublanovskaya (1961). It is one
fundamental algorithm often implied in the other methods.
Perimeter of use: Calculation of all the spectrum.
Advantages: Good convergence, robustness, calculation of the form of Schur (any matrix
complex
With
the decomposition of Schur admits
Q WITH Q T
*
=
, with
Q
and
T
matrices
respectively, unit and triangular higher. In the real case,
T
is only diagonal by
block 1x1 or 2x2 (real form of Schur). Columns of the matrix unit, known as of Schur,
vectors of Schur are called).
Disadvantages: Prohibitory memory complexities and calculation, sensitivities to the “balancing (them
techniques of balancing (English balancing) consist in transforming reversibly
the operator of work so that its matric standard decreases. Thus its handling will be
less sensitive to the effects of round (cf [§10.3]))“.
Alternatives: With implicit shift or clarifies, simple or double…
·
The algorithms of the powers type which were historically developed the first for
to solve generic modal problems. These are basic algorithms of which them
others are an improvement. They are implied in the operator
MODE_ITER_INV
(cf [§3]).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
18/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Perimeter of use: Calculation of the extreme values of the spectrum.
Advantages: Simplicity, very good estimate of the clean vector in some iterations.
Disadvantages: Convergence can become problematic, bad capture of
multiplicities, of the clusters… (together of eigenvalues “very close”).
Alternatives: Powers opposite, (Bi) iteration of the quotient of Rayleigh (cf [§3.3.2])…
·
methods of subspace which consist in projecting the operator of work on a space
H
such that the spectrum of the projected operator is a good approximation of the part of
initial spectrum which one seeks
.
These algorithms are the hard core of the operator
MODE_ITER_SIMULT
(cf [§4]).
Perimeter of use: Calculation of part of the spectrum.
Advantages: Reduction of the size of the problem and memory complexities and calculation,
require only the calculation of a product matrix-vector and not the knowledge of
all the matrix.
Disadvantages: Use the many pre one and postprocessings, convergence can become
problems, more or less easily capture the multiplicities and the clusters according to
alternatives.
Alternatives: Iterations of subspace, Bathe and Wilson (1971), Lanczos (1950), Arnoldi
(1951), Davidson (1975), Sorensen (1992)…
·
other approaches more or less empirical and are specialized. They are often
connected to other problems: search for roots of polynomials, functions
unspecified… One can thus quote the method of bisection used in preprocessings in
MODE_ITER_INV
, but also that of Jacobi, Laguerre etc…
Note:
Many parallels can be led between these families (the method
QR
is not thus
that a method of iterations of subspaces applied to entire space), but they
also lead to processes similar to those developed for other problems:
·
optimization; The method of the quotient of Rayleigh is with the method of the powers opposite,
what the method of Newton is for a method of descent conventional,
·
resolution of linear systems; The method of the combined gradient is a method of
subspace for the positive definite symmetrical systems,
·
search for roots of polynomials: the method of the powers is a method of
Bernoulli applied to the matrix “companion” of the polynomial (cf [bib11] pp502/503).
The following paragraph will synthesize what precedes in the total flow chart by resolution of one
generalized modal problem of Code_Aster.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
19/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
2.9
Establishment in Code_Aster
The aforementioned can break up into four phases:
1)
The first operation consists in determining the shift like certain parameters of the problem.
That is carried out in a more or less transparent way according to the option of calculation chosen by
the user:
·
MODE_ITER_INV
+ option: '
SEPARATED
“or”
AJUSTE'
the shift is determined by the first
phase of the algorithm and the number of clean modes sought by frequential tapes
(provided by the criterion of Sturm) is limited by
NMAX_FREQ
,
·
MODE_ITER_INV
+ option: '
NEAR
'
the shift is fixed by the user and the number of modes
clean is equal to the number of shifts,
·
MODE_ITER_SIMULT
+ option:
“PLUS_PETITE”
shift is null and the number of modes is
parameterized by
NMAX_FREQ
,
·
MODE_ITER_SIMULT
+ option: '
BANDAGE
'
shift is equal in the middle of the tape fixed by
the user and the number of modes are determined by the criterion of Sturm,
·
MODE_ITER_SIMULT
+ option: '
CENTER
'
the shift is fixed by the user and the number of
clean modes is parameterized by
NMAX_FREQ
.
2)
In one second phase of preprocessings, one factorizes partially the matrix of work
(
)
With
With
B
B
=
-
-
1
while being interested only in its reversed part. The action of this operator on an unspecified vector
U
,
noted
U
2
, will build itself thus more effectively
(
)
U
B U
With
B U
U
U
1
2
1
2
=
-
=
This factorization undergoes the same risks besides as the criterion of Sturm when the shift is
near to an eigenvalue. One carries out the same shifts then according to algorithm 1.
Notice
In V5 when
NMAX_ITER_SHIFT
calculation is reached stops in fatal error
(only for this preprocessing), this problem of factorization which can give place to
numerical instabilities.
3)
Modal calculation itself is carried out: the standard problem is solved, then one returns to
initial problem. In
MODE_ITER_INV
two alternatives are available: method of the iterations
opposite (option: “DIRECT”) and its acceleration by quotient of Rayleigh (“RAYLEIGH”). For it
who is of
MODE_ITER_SIMULT
, it allows the use of three distinct methods: method of
Bathe and Wilson (“JACOBI”), that of Lanczos (
'TRI_DIAG
') and finally, that of Sorensen
(“SORENSEN”).
Each one of these methods have internal tests of stops. Without counting that methods
of projection employ auxiliary modal methods:Jacobi (cf [Appendix 3]) for the first
and
QR QL
/
(cf [Appendix 1]) for the others. They require also tests of stops.
The user often has access to these parameters, although it is warmly recommended, at least
initially, to preserve their default values.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
20/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
4)
This last part gathers the tests of total stops which check the good unfolding of
calculation. They are of two types:
·
(One normalizes recalls that
U
U
=
max
I
I
,
With
With
=
max
I
ij
J
,
U
U
2
2
1
2
=
I
I
and
(
)
()
With
AA
With
With
2
1
2
=
=
max
*
max
I
I
I
I
if
square
residue of the initial problem
U
U
U
With U
B U
With U
With U
B U
>
-
<
-
<
If S
then
If not
End if.
EUIL FREQ
THRESHOLD
THRESHOLD
_
?
,
?
,
2
2
2
Algorithm 2: Test of the standard of the residue
This sequence is parameterized by the key words
THRESHOLD
and
SEUIL_FREQ
, pertaining
respectively with the key word factor
VERI_MODE
and
CALC_FREQ
. In addition, it
postprocessing is activated by initialization with “YES” (default value) of
STOP_ERREUR
in
the key word factor
VERI_MODE
. When this rule is activated and non-observance, calculation
arrète, if not the error is just announced by an alarm. One would know of course only too
to recommend not to decontaminate this parameter preferential treatment!
·
Counting of the eigenvalues
This postprocessing is set up in
MODE_ITER_SIMULT
and it makes it possible to check that it
a number of eigenvalues contained in a tape test
[
]
1
2
,
is equal to the number
detected by the algorithm. The inclusion of the initial tape (in option
“TAPE”
it is about
values indicated by the user, if not they are the extreme values of the calculated spectrum)
[
]
I
F
,
in this tape test is led in order to detect possible problems of
clusters or of multiplicities at the initial boundaries (cf TP n°1 [bib25]).
Appear 2.9-a: Counting of the eigenvalues
I
i+
1
F
2
F
Initial tape
Bandage tested
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
21/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
While noting
I
+
and
F
-
, respectively, largest and the smallest eigenvalue not requested by
the user and including the initial tape (cf [Figure 2.9-a]), one has the algorithm of construction of
bandage following test:
(
)
()
(
)
()
(
)
(
)
If
exist resp.
and if
resp.
then
=
2
resp.
=
2
,
If not
=
1 - sign
resp.
=
1 + sign
End if.
i+
F
1
I
I
F
F
I
F
I
I
I
F
F
F
I
F
PREC SHIFT
+
-
+
-
-
<
-
+
+
_
,
1
2
2
PREC_ SHIFT
PREC_ SHIFT
Algorithm 3: Construction of the tape test
This sequence is parameterized by key word PREC_SHIFT of the key word factor VERI_MODE.
Framed below the trace of the error messages shows when these post-checkings
start.
CHECKING A POSTERIORI OF THE MODES
<E> <MODE_ITER_SIMULT> <VERIFICATION OF THE MODES> FOR THE CONCEPT >MOD_4< IT
MODE NUMBER 2 OF CRITICAL LOAD 8.7243E+07 A a STANDARD Of ERROR
OF 6.3919E-01
<E> <MODE_ITER_SIMULT> <VERIFICATION OF THE MODES> FOR THE CONCEPT >MOD_4< IN
The INTERVAL (- 1.0608E+08,1.8130E+08) IT THERE A THEORETICALLY 3 LOAD (S)
CRITICAL (S) AND ONE OF A CALCULEE (S) 4.
Example 3: Post-checkings
The activation of the second postprocessing is subordinated to the initialization of
STURM
with
“YES”
in
key word factor
VERI_MODE
. When this rule is activated and not respected, the continuation of the events
is controlled by
STOP_ERREUR
(if '
OUI'
calculation stops, if not the error is just announced by one
alarm). It would not be known of course that too much to recommend not to decontaminate these parameters 'passes
droit'!
Now that the context of the generalized modal problems of Code_Aster was brushed, us
let us be interested more particularly in the methods of the power type and their establishment
in the operator
MODE_ITER_INV
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
22/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
3
Method of the powers opposite (
MODE_ITER_INV
)
3.1 Introduction
To calculate several eigenvalues of the generalized problem, the general method is not used
iterations opposite just as it is.
One can, for example, to combine it with a technique deflation (technique of consistent filtering with
to transform the operator of work so that it has the same eigenvalues except some
preset modes which see themselves affecting a zero value) so as to filter automatically
spectral information updated more to find it with the following iteration. With deflation by
restriction (of other types of deflation exist, such as the method of Ducan-Collar which uses for
to filter spectral information the first line of the matrix and the clean vector. These techniques
vectorial spread of course per blocks) of Wielandt one builds the operator repeatedly of
work (in the symmetrical case), while noting
(
µ
K
,
U
K
)
mode to be filtered,
With
With
U U
K
K
K
K
K
T
+
=
-
1
µ
.
This strategy, which does not apprehend the multiplicities, must be supplemented by a criterion of Sturm.
In addition, the fact of working into arithmetic finished and of not building the operator indeed
With
K
with each iteration, constrained with mâtiner this step of process of orthogonalization
powerful.
All these complications resulted in choosing another channel, which breaks up into two parts:
1) localization of the eigenvalues (determination of an approximate value of each value
clean contained in an interval given by a technique of bisection, refined or not,
by a method of the secant),
2) improvement of these estimates and the calculation of their own vectors associated by one
method of iterations opposite.
The search for a value approached for each eigenvalue considered is selected in
key word factor
CALC_FREQ
via
OPTION
:
·
if
OPTION
=
“SEPARATE”
, in each interval of frequencies defined by the key word
FREQ
,
an approximate value of each eigenvalue contained in this interval is calculated in
using the method of dichotomy (cf [§3.2.1]),
·
if
OPTION
= '
ADJUST
', one carries out the same operations first of all as previously and
then, on the basis of these approximations, one refines the result by the method of the secant
(cf [§3.2.2]).
For these two options, one calculates at the same time the modal position of each value
clean what makes it possible to detect the multiple modes. Either only they are retained
NMAX_FREQ
the lowest frequencies contained in the maximum interval specified by the user, is
one calculates all the values of this interval (if
NMAX_FREQ
= 0).
·
if
OPTION
=
“NEAR”
, frequencies given by the key word
FREQ
, are considered
like the approximate values of the sought eigenvalues.
Note:
·
of course, as one already specified (cf [§2.8]), this operator is to be used only for
to determine or refine some eigenvalues. For a wider search it is necessary
to use the operator
MODE_ITER_SIMULT
,
·
with the option '
NEAR
'one cannot calculate multiple modes.
·
it is an expensive algorithm because it calls much upon the test of Sturm and thus with its
associated factorizations.
·
the terminals of the intervals of search are provided by
FREQ
or
CHAR_CRIT
according to
the initialization of
TYPE_RESU
.
We now will detail the various algorithms (and their parameter settings) which are put in
place in the first part of the process.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
23/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
3.2
Localization and separation of the eigenvalues
3.2.1 Method of bisection
As one already saw previously, the corollary 5bis of the law of Inertia of Sylvester (cf [§2.5])
allows to determine the number of eigenvalues contained in an interval given while carrying out
two decompositions
LDL
T
. This criterion can thus result in refining the interval until not including
that an eigenvalue. This piloting being set up, one passes from one iteration to the other by using it
principle of the dichotomy.
On a starting interval
[
]
has
B
,
, one thus operates in the following way, knowing
()
pm
has
and
()
pm
B
:
(
)
()
()
()
()
(
)
[]
[]
(
)
[]
[]
Calculation of
Calculation of pm
If pm
pm
resp. pm
then
to start again on
resp.
If not
to start again on
and on
End if.
*
*
*
*
*
*
*
.
.
,
,
,
,
,
,
=
+
=
1
2
has
B
has
B
B
has
has
B
Algorithm 4: Method of Bisection
The process is stopped if one cut out more
NMAX_ITER_SEPARE
time the starting interval, or if
for a given interval, one has
has
B
SEPARATE PREC
-
*
_
(in this case one does not refine any more
seek in this interval).
One obtains finally a list of frequencies. In each interval defined by the arguments of this
list, is an eigenvalue having a certain multiplicity. Like approximation of this value
clean, the medium of the interval is taken.
Note:
·
one could have used as criterion the change of sign of the characteristic polynomial, but
in addition to the fact that it is very expensive to evaluate, it does not make it possible, just as it is, to detect them
multiplicities,
·
the initialization of the process can be carried out in an empirical way according to the needs for
the user. To include part of the spectrum one can also use the régionnements
plan complexes theorems of Gerschgörin-Hadamard (on
WITH A
,
T
…). Accordingly,
the method of bisection can prove more effective than one
QR
in the presence of cluster. Its
convergence, although linear, is indeed raised by ½ whereas that of
QR
can tend
towards 1 [bib11].
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
24/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
has
B
*
1 era division
has
B
*
p (
)
B
has
*
B
p (
)
B
has
has
*
*
*
B
B
p (
)
B
has
has
*
*
B
has
B
*
*
has
the 2nd division
the 3rd division
the 4th division
has
Interval of will tra
Appear 3.2.1-a: Method of bisection
These approximated eigenvalues can be improved by a unidimensional minimization
seeking to cancel the characteristic polynomial
()
(
)
p
dét
=
-
With
B
. But only the calculation of one
value of this determinant requires NR! operations.
3.2.2 Method of the secant
The method of the secant is a simplification of the method of Newton-Raphson
.
With the stage
unspecified
K
, knowing a value
K
and by approximating the nonlinear function
()
p
by its
tangent in this point, one determines the next value
k+1
as being the intersection of this line
with the axis of
, and so on, according to the iterative diagram
() () ()
K
K
K
K
K
K
K
p
p
p
+
-
-
=
-
-
-
1
1
1
Appear 3.2.2-a: Method of the secant
The tangent being approximated by a difference finished in order not to have to calculate of derived from
()
p
, only the estimate of the polynomial is necessary.
k-1
k+1
K
p
(
)
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
25/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
It is considered that one reached convergence when
K
K
K
PREC ADJUSTS
+
-
<
1
_
and, in addition, one
limits itself to
NMAX_ITER_AJUSTE
iterations if this criterion is not reached.
Note:
This method has a quadratic convergence when it is close to the solution, in the case
opposite, it can diverge. From where interest to combine the method of bisection with this
approach.
We now will detail the algorithm of the powers opposite (coupled to an acceleration of
Rayleigh) constituting the second part of the process.
3.3
Method of the powers opposite
3.3.1 Principle
To determine the eigenvalue of the generalized problem
With U
B U
=
nearest in module to
, one applies the method of the powers to the operator
(
)
With
B
B
-
-
1
. In fact, one only builds
the factorized matrix (this notation should not be confused that symbolizing the “shift and invert”
(
)
With
With
B
B
=
-
-
1
)
(
)
With
With
B
=
-
and that amounts dealing with the generalized problem
()
With
Drunk
U
µ
-
= -
1
1
!“$
. Method of the powers opposite converging proritairement towards
eigenvalues of stronger module (in
µ
), they thus will be captured
closest to the shift.
The principle is as follows, knowing an estimate
eigenvalue sought and on the basis of one
standardized initial vector
y
0
, one builds a clean vector series approximate
()
y
K K
by
recurring formula
To make
Fine loops.
K
K
K
K
K
K
K
K
=
=
=
=
-
1
1
,…
~
,
~,
~
,
With y
B y
y
y
y
µ
µ
Algorithm 5: Method of the powers opposite
With
1
,
2
… first eigenvalues (of clean vector
U
I
) closest in module to
,
it is shown that one has a linear convergence of
y
U
K
1
and a quadratic convergence of
() ()
(
)
µ
K
K
K
I
I
I
N
-
-
=
-
1
1
1
1
1
1
and
~
.
y
y
(the factor of convergence of these continuations is of
the command of
1
1
-
-
min
I
I
).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
26/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Note:
·
this result is not acquired (with the factor of convergence near) only when the eigenvalue
dominant (here that nearest in module to the shift) is single. In the contrary case,
results remain however possible, even in the complex framework, but the analysis
rigorous of the convergence of this algorithm is still incomplete [bib6], [bib23],
·
in theory, these results require that the initial vector is not orthogonal with
clean subspace on the left required. In practice, the round-off errors avoid it
problem (cf [bib11] pp500-509),
·
even if the estimate of the eigenvalue is coarse, the algorithm provides one quickly very
good estimate of the clean vector,
·
the major disadvantage of this method is that it is necessary to carry out a factorization of
With
for
each eigenvalue to calculate.
In general one works in euclidian norm or infinite standard, but to facilitate calculations
post-modal one seeks here with
B
- to standardize the clean vectors (when
B
is indefinite, one works
with the associated pseudo-standard). The basic algorithm can be rewritten while posing
Z
B y
0
0
=
()
(
)
For K 1… to make
Fine loops.
=
=
=
=
=
=
-
-
With y
Z
Z
B y
y
y
Z
y
Z
y
Z
Z
Z
y
Z
K
K
K
K
K
K
T
K
K
T
K
K
K
T
K
K
K
K
K
T
K
sign
1
1
1
2
,
~
,
.
. ~,
. ~,
~
. ~
,
Algorithm 6: Method of the powers opposite with
B
- pseudo-standard
Let us note that
()
y
K
-
1
1
. This presentation avoids the products matrix-vector by the matrix
B
during the calculation of the scalar products and already the coefficient of Rayleigh of the following paragraph precedes.
It is observed that the factor of convergence is all the more small as the spectral shift
is close
sought eigenvalue and thus that
With
B
-
is close to the singularity. That is not in fact not
prejudicial with the process because the error made by solving the system is “mainly” in
direction generated by the clean vector which is the sought direction. That means that at the time of
resolution of
With y
Z
K
K
=
-
1
, one does not find the solution exact
y
K
but that round-off errors
lead to solution close to the form
~
y
y
W
K
K
=
+
. The aforementioned is proportional to the solution
exact, but as standardization is arbitrary, all happens correctly [bib18].
Note:
This bad conditioning, far from having an unfavourable effect, improves even convergence of
the algorithm.
This algorithm is thus used to improve the clean vector associated with the value approximated with
phase 1. To refine this estimate of the eigenvalue, a quotient of Rayleigh is introduced.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
27/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
3.3.2 Method of iteration of the quotient of Rayleigh
Let us recall that the quotient of Rayleigh applied to the generalized problem is defined by the number
R (X)
, with
X
a vector not no one of
N
, such as:
()
R X
X With X
X B X
=
T
T
This quotient has the remarkable property of stationnarity in the vicinity of any clean vector and
to reach a extremum (local) which is the corresponding eigenvalue: for each
X
fixed
, R (X)
minimize
(
)
-
With
B X
2
.
What we can translate by “if
X
is an approximation of a clean vector of the system
With X
B X
=
, then
R (X)
is an approximation of the eigenvalue associated the vector
X
“, and
reciprocally, we saw that if one had a good estimate of an eigenvalue,
method of the iterations opposite made it possible to obtain a good estimate of the clean vector
corresponding.
From where the idea to combine these two properties by considering the algorithm of iteration reverses with
spectral shift for which one revalues, with each iteration, the eigenvalue via the quotient of
Rayleigh. One then obtains the algorithm known as of iterations of the quotient of Rayleigh (in
B
- pseudo-standard)
()
(
)
()
()
(
)
For K 1… to make
Fine loops.
=
-
=
=
=
+
=
=
-
-
-
-
With
y
B y
Z
Z
B y
y
y
Z
y
Z
y
y
Z
Z
Z
y
Z
R
,
~
,
R
.
. ~
R
,
. ~,
~
. ~
,
K
K
K
K
K
K
K
T
K
K
T
K
K
K
K
T
K
K
K
K
K
T
K
sign
1
1
1
1
1
2
Algorithm 7: Method of iterations of the quotient of Rayleigh with
B
- pseudo-standard
For the standard modal problem, one can show [bib18] that the convergence of this algorithm is
cubic if the operator of work is normal (a matrix
With
is known as normal if
AA *
=
WITH * A
.
It is thus the case of the operators square, antihermitiens or unit) (a fortiori in the case
square) and at best quadratic in the other cases.
If one used this method without modifying it, it would be necessary for each iteration of the process of improvement
of each eigenvalue, to carry out a factorization
LDL
T
, which would be very expensive. From where the idea of
not to even carry out (this strategy of the coupling of methods to the complementary characteristics
antagonists is often used in numerical analysis. For example, in optimization, the method of
Levenberg-Marquadt couples one steepest-descent and a Newton) this shift of Rayleigh, that if one
is in a vicinity (arbitrary concept to define) solution.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
28/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Notice
Within a more general framework, certain authors introduced an algorithm known as “of Bi-iterations of
quotient of Rayleigh ". Based on the stationnarity of the Bi-quotient
()
R
,
B
T
T
X y
y With X
y B X
=
in the vicinity
clean vectors on the right and on the left, it provides (in nonsquare) the two types of vectors
clean. Its cost is however crippling because it requires twice more factorizations
(cf B.N. Parlett, 1969 [bib18]).
3.3.3 Establishment in Code_Aster
This spectral shift is activated only if the key word
OPTION
key word factor
CALC_MODE
is initialized
in “RAYLEIGH”. By defect, one has '
DIRECT
'and the shift is then conventional (total correction rather
that progressive of the eigenvalue). The algorithm set up in the code cuts out as follows
(in
B
- pseudo standard):
·
Initialization of the eigenvalue starting from the estimate of the first phase:
.
·
Calculation of an initial vector
y
0
random checking the limiting conditions.
·
B
- orthonormalisation of
y
0
compared to the modes previously calculated (if it is a mode
multiple according to the first phase) by Gram-Schmidt Modified (since the development
of this operator, more effective and more robust methods of orthogonalization are
spread, such as Gram-Schmidt Modified Iterative (IGSM) used in
MODE_ITER_SIMULT
(cf [Appendix 3] and B.N.Parlett [bib18])) (GSM).
·
Calculation of
(
)
0
0
0
=
sign
T
y By
.
·
For
k=1,
NMAX_ITER
to make
To solve
(
)
With
B y
B y
-
=
-
-
K
K
K
1
1
.
B
- orthonormalisation (possible) of
y
K
.
Calculation of the correction of the eigenvalue
C
K
T
K
K
T
K
=
-
y B y
y B y
1
.
If
y B y
K
T
K
PREC
-
-
1
1
then
=
+
C
,
exit;
If not
If OPTION = “RAYLEIGH” and if
C
01
.
then
=
+
C
;
End if.
End if.
Fine loops.
Algorithm 8:
MODE_ITER_INV
The standard of maximum error acceptable
PREC
and the maximum number of authorized iterations
NMAX_ITER
are arguments of the key word factor
CALC_MODE
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
29/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Note:
·
the clean vector being
B
- normalized, one considers to have reached convergence when the value
absolute of the unsuitable scalar product is close to the unit,
·
to avoid taking an initial vector
B
- orthogonal with the sought eigenvalue one uses
a method of hard copy random for the components of this vector,
·
in addition to be able to determine multiple or close modes, one is used
B
- orthogonalization with the already calculated modes,
·
the option “of acceleration” of the algorithm by quotient of Rayleigh being expensive, it is not
used with each iteration that if one is in the vicinity of the required eigenvalue.
3.3.4 Display in the file message
In the file message, the results are displayed in the form of table
Position
modal
Frequency Deadened
sow
Dichotomy
Method
Secant
Method
Opposite
Numbers
iterations
Numbers
iterations
Precision Numbers
iterations
Precision
*
*
0.
*
*
*
*
*
Table 3.3.4-a: Traces
MODE_ITER_INV
in the file message
------------------------------------------------------------------------
THE NUMBER OF DDL
TOTAL EAST: 192
LAGRANGE EAST: 84
THE NUMBER OF ACTIVE DDL EAST: 66
------------------------------------------------------------------------
CHECKING OF THE FREQUENCY SPECTRUM (CONTINUATION OF STURM)
THE NUMBER OF FREQUENCIES IN THE TAPE (1.00000E-02, 6.00000E-02) IS 4
------------------------------------------------------------------------
MODAL CALCULATION:
METHOD Of ITERATION OPPOSITE
OPPOSITE SECANT DICHOTOMY
NUMBER FREQUENCY (HZ) AMORT NB_ITER/NB_ITER/PRECISION/NB_ITER/PRECISION
4 1.97346E-02 0.00000E+00 4 6 2.97494E-07 4 1.22676E-07
5 2.40228E-02 0.00000E+00 4 5 4.21560E-05 3 4.49567E-09
6 4.40920E-02 0.00000E+00 3 5 2.19970E-05 3 2.62910E-09
7 5.23415E-02 0.00000E+00 3 5 2.34907E-07 5 1.32212E-07
------------------------------------------------------------------------
CHECKING A POSTERIORI OF THE MODES
------------------------------------------------------------------------
Example 4:
MODE_ITER_INV
With the option
“NEAR”
, the columns “Dichotomy” and “Secant” do not appear, while with
the option
“SEPARATE”
, only the “Secant” column disappears. The last column precision gathers
intermediate data and does not represent, as for the other modal operator
MODE_ITER_SIMULT
, the standard of error of the residue. It is an artefact which will be brought to disappear.
Let us recapitulate the parameter setting of the operator now
MODE_ITER_INV
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
30/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
3.3.5 Summary of the parameter setting
Key word factor
Key word
Default value
References
“DYNAMIC” TYPE_RESU
“DYNAMIC”
[§2.1]
“MODE_FLAMB”
[§2.1]
CALC_FREQ
FREQ
[§3.1]
CHAR_CRIT
[§3.1]
OPTION “SEPARATE”
“ADJUSTS”
[§2.1]
“ADJUSTS”
[§3.1]
“NEAR”
[§3.1]
NMAX_FREQ
0
[§3.1]
NMAX_ITER_SEPARE
30
[§3.2.1]
PREC_SEPARE
1.E-04
[§3.2.1]
NMAX_ITER_AJUSTE
15
[§3.2.2]
PREC_AJUSTE
1.E-04
[§3.2.2]
NPREC_SOLVEUR
8
[§2.6]
NMAX_ITER_SHIFT
5
[§2.6]
PREC_SHIFT
0.05
[§2.6]
SEUIL_FREQ
1.E-02
[§2.9]
CALC_MODE
OPTION “DIRECT”
“DIRECT”
[§3.3.3]
“RAYLEIGH”
[§3.3.3]
PREC
1.E-05
[§3.3.3]
NMAX_ITER
30
[§3.3.3]
VERI_MODE
STOP_ERREUR “YES”
“YES”
[§2.9]
“NOT”
THRESHOLD
1.E-02
[§2.9]
Summary table 3.3.5-a: of the parameter setting of
MODE_ITER_INV
Note:
·
in V5, the user can specify the class of membership of his calculation by initializing it
key word
TYPE_RESU
. According to this value, one informs the vector
FREQ
or
CHAR_CRIT
,
·
one finds all the “tripaille” of parameters related to the preprocessings of the test of Sturm
(
NPREC_SOLVEUR
,
NMAX_ITER_SHIFT
,
PREC_SHIFT
) and with postprocessings of checking
(
SEUIL_FREQ
,
VERI_MODE
),
·
at the time of the first passages, it is strongly advised to modify only them
main parameters noted in fat. The others relate to more the mysteries of
the algorithm and they were initialized empirically with values standards.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
31/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
4
Method of subspace (
MODE_ITER_SIMULT
)
4.1 Introduction
If one only wishes to calculate the p clean elements of a generalized problem of command
N
where
N
p
<<
(for example, p smaller eigenvalues or all eigenvalues included/understood in
a given interval), one saw that it was preferable to have resorts to methods of subspace. They
are based on the analysis of Rayleigh-Ritz which consists in effectively projecting the problem considered
on under space of dimension
m
(
p < m < N
) and to seek certain clean elements of it
problem projected (much easier to treat) using robust algorithms (
QR
or
QL
for
Lanczos and IRAM, Jacobi for the method of Bathe and Wilson).
The criteria of effectiveness of the aforesaid projection are:
·
small size of the space of projection (directly related to calculation complexities and memory),
·
facility of its construction,
·
the robustness of orthogonal projection (in the nonsquare case, of oblique projections
were developed by F. Chatelin [bib2],
·
the setting in the canonical shape of the projected matrix,
·
the good approximation of the part of the initial spectrum sought by that of the projected operator,
·
and, of course, the minimization of calculation complexities and memory and that of the effects of flares
(those are especially related to the problems of orthogonality raised by the 3
ème
not).
One will first of all present the analysis of Rayleigh-Ritz before detailing three methods resulting from this
analyze: method of Lanczos, that known as of Sorensen (WILL GO) and that of Bathe and Wilson.
4.2
Analyze of Rayleigh-Ritz
Let us consider the standard modal problem of command
N
,
With U
U
=
, and the subspace
H
of
N
generated by a orthonormée base
(
)
Q
Q
Q
1
2
,
,…,
m
. The latter constitutes the matrix
orthogonal
Q
m
allowing to define the operator of projection
P
Q Q
m
m
m
T
=
. Method of
Galerkin used by this analysis consists in solving the following problem
()
(
)
()
To find
such as
With
to find such as
~, ~
~
~~
~
,
~, ~
*
U
P With
U
0
U Q X
X
Q AQ
B
X
X
×
-
=
=
×
=
H
m
m
m
m
m
m
! “
# $
#
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
32/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
The search of the elements of Ritz
()
~, ~
U
who thus interest cost to us to find the modes
clean of the matrix of Rayleigh
B
Q AQ
m
m
T
m
=
. The eigenvalues remain unchanged in
two formulations, on the other hand knowing a clean vector
X
of
B
m
one must go up with space all
entirety via the transformation
~
U Q X
=
m
. The step can thus be summarized in the form:
·
Choice of a space
H
and of a base
(H
1
, H
2
,…, H
m
)
represented by
H.
·
Orthonormalisation of the base
·
Stamp of Rayleigh
·
Clean elements of
B
m
:
()
~,
X
·
Elements of Ritz:
(
)
~, ~
U Q X
=
m
,
·
Test of error with the residue:
R
With
U
=
-
~
~~
.
Algorithm 9: Procedure of Rayleigh-Ritz
Note:
·
the elements of Ritz are in fact the clean modes of the matrix (of command NR) of
the approximation of Galerkin
C
m
=
P
m
WITH P
m
,
·
calculation complexity of this process, the command at worst of
O (m
2
(4n+m))
(much less
with a good space, cf Lanczos and IRAM), are without common measurement with that of a good
QR
(
O (N
2
)
) or that of an iteration of the quotient of Rayleigh encapsulated in a process
heuristics (the cost of factorizations is prevalent on that of the products matrix-vector
algorithm itself) such as that developed in
MODE_ITER_INV
(>> pn
3
). But
it is quite clear that these methods are more complementary than concurrent.
To consider the error made by using the clean elements of the matrix
B
m
there is the result
according to.
Theorem 6
That is to say
()
, U
a clean element of the matrix diagonalisable (it is the case of the normal matrices and of
nondefective matrices (matrices of which all the eigenvalues have even arithmetic multiplicity and
geometrical)
With
and,
B
m
the matrix of Rayleigh associated, then the aforementioned admits an eigenvalue
~
checking
(
)
-
-
~
I P
U
P U
m
m
2
2
2
where
is a number dependant on
,
With
and
P
m
.
,
U
H
P
m
U
~ ~
, U
(I - P
m
) U
X
H
=
Q
m
R
0
N
m
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
33/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
In the square case the numerator of this increase is squared. The associated clean vector,
~
U
,
check as for him
()
(
)
sin, ~
U U
I P
U
P U
-
m
m
2
2
.
Proof:
Cf [bib11] pp531-537.
One shows that for
m
enough large,
can be raised by a number which depends only on
With
and of
min
I
J
I
J
-
. The estimate of the second member
(
)
I P
U
-
m
2
depends on the choice of the subspace.
Note:
The number
min
I
J
I
J
-
(and its alternatives) is omnipresent in the analyzes of errors, of
convergence or of conditioning spectral. One sees, once more, the importance of
separation of spectrum of work on the good numerical behavior of the algorithms.
Of course, to take account of spectral information already obtained, even to improve it or filter it,
one reiterates this procedure of Rayleigh-Ritz by modifying the subspace of projection. One will connect
then the construction of this new space (recurring of the precedent) and this process of projection. Two
types of spaces (each one attached to distinct methods) are most often used.
4.3
Choice of the space of projection
Two choices of space of projection
H
are most often set up:
·
The first, that of the method of Bathe and Wilson (cf [§7]
METHOD = “JACOBI”
),
consist in choosing under space
H
of dimension
m
then to build successively:
H
WITH H
H
WITH H
WITH H
H
WITH H
WITH H
1
2
1
2
1
=
=
=
=
=
-
…
I
I
I
This method, which holds at the same time of generalization in form block of the method of
powers and of the truncation of the algorithm
QR
, led to an impoverishment of space
of work:
dim
dim
H
H
I
I
-
1
. Not to always find the same clean mode
dominating it is necessary to insert in the process a reorthogonalisation (with all the problems
of calculation complexity and flares that that implies).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
34/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
·
The second, that of the method of Lanczos (cf [§5],
METHOD = “TRI_DIAG”
) and of IRAM
(cf [§6],
METHOD = “SORENSEN”
) consists starting from an initial vector
y
, to build
continuation of spaces of Krylov (
H
I
)
I = 1, m
where
H
I
is the space generated by the family of vectors
(
)
y Ay
With
y
,
,…,
I
-
1
. This last is called the space of Krylov of
With
and of command
I
initiate by
y
:
()
(
)
H
K
Vect
I
I
I
=
-
With y
y Ay
With
y
,
,
,…,
1
.
They check the following property:
dim
dim
H
H
I
I
+
1
. Contrary to the preceding choice, one
thus a certain enrichment of the workspace has. In addition, it will be seen that they
allow to project in an optimal way within the meaning of the criteria defined previously.
Historically, the first type of space was very much used in mechanics of the structures. But for
to decrease calculation complexity related on the size of the subspaces and the orthogonalizations, one prefers to them
from now on methods of the type Lanczos/Arnoldi.
Note:
One meets an impressive variety of algorithms using a space of the first type. They
are called “'iterations of subspace”, “iterations orthogonal” or “iterations
simultaneous ". Beyond the various terms, it should especially be retained that these adaptations
the strategies of restartings (technique concern consisting in starting again an algorithm
modal with an initial vector comprising spectral information already. Typically, one chooses
a linear combination of the clean vectors or already exhumed vectors of Schur
(cf [§5.4.2]), techniques of acceleration (technical consisting in replacing the matrix of work
With
by a matric polynomial
P (A)
who with the characteristic to be of great amplitude in the areas
spectral of interest (cf [§6.4]) and of factorization implemented to enrich the workspace.
There are even versions based on the powers opposite making it possible to calculate them
p
more
small modes (cf [bib11] pp538-45, [bib23] pp49-54).
For these methods of projection, by supposing that initial space
H
is not too poor according to
p
dominant directions, the factor of convergence of
I
ème
mode When they are arranged
classically, i.e. by order descending of module (when they are classically arranged,
i.e. by order descending of module) at the end of
K
iterations is written:
O
m
I
K
+
1
.
It expresses two phenomena clearly:
·
the priority convergence of the dominant modes,
·
improvement of these convergences (and thus of their standards of error for a number
iterations fixed) when the size of the subspace increases.
To transform the modal problem generalized into a standard problem one A resorts to
spectral transformations. They also make it possible to direct the search of the spectrum and to improve
convergence (cf [§2.7]).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
35/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
4.4
Choice of the spectral shift
To calculate the smallest eigenvalues close to a given frequency or all the values
clean included/understood in a given frequency band, one carries out a spectral shift of
generalized problem. That is to say
the value of shift, one transforms the initial problem
With U
B U
=
in
a shifted standard problem.
WITH B U
U
µ
-
=
1
with
With
With
B
= -
and
µ
= -
This transformation spectral, known as of “simple shift”, is used in the method of Bathe and Wilson.
As the process detects the smallest eigenvalues in
µ
, those gradually are captured
who are closest to
(in module).
One carrying out the reverse transformation (“shift and invert”) it occurs
With U
U
µ
=
with
(
)
With
With
B
B
=
-
-
1
and
µ
= -
1
It is the problem dealt with with Lanczos and IRAM which makes it possible to calculate the greatest eigenvalues
µ
thus the eigenvalues closest to the shift
(in module).
In the control
MODE_ITER_SIMULT
, there are three ways of choosing this shift (for buckling,
transposition at the levels of critical loads is immediate):
·
=
0
, one seeks the smallest eigenvalues of the starting problem. This corresponds to
OPTION = “PLUS_PETITE”
under the key word factor
CALC_FREQ
.
·
=
0
with
(
)
0
0
2
2
=
F
, one seeks the frequencies close to the frequency
FREQ
=
F
0
.
(
OPTION = “CENTER”
).
·
=
+
0
1
2
with
(
)
0
0
2
2
=
F
and
(
)
1
1
2
2
=
F
, all the frequencies are sought
included/understood in the tape
[
]
F
F
0
1
,
(
OPTION = “TAPE”
with
FREQ
= {
F
0
, F
1
}).
The number of frequencies to be calculated is given in general by the user using
NMAX_FREQ
under
the key word factor
CALC_FREQ
. But for the option
“TAPE”
, it is automatically given in
using the corollary 5bis (cf [§2.6]).
Note:
It is reminded the meeting that the factorization of the matrix of work, just like the tests of Sturm of the option
“TAPE”
are dependant on numerical instabilities when the shift is close to an eigenvalue.
Palliative, skeletal processing by the user, were established (cf [§2.6], [§2.9]).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
36/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5
Method of Lanczos (
METHOD = “TRI_DIAG”
)
5.1 Introduction
This algorithm, developped at the point by Lanczos [bib10] in 1950, was hardly used until the medium of
the Seventies. First of all very simple and effective, it is the seat nevertheless of great instabilities
numerical being able to cause the capture of phantom modes! Those are dependant mainly on
problems of hold of orthogonality enters the vectors generating the space of Krylov.
comprehension of this type of behavior vis-a-vis arithmetic finished computers was long with
to exhume.
Since, many palliative solutions were born and an abundant literature covers the subject.
Let us quote for example the work of J.K. Cullum [bib3] and Al which is entirely devoted to him, that of
B.N. Parlett [bib18] and the brought up to date and exhaustive synthesis of J.L. Vaudescal [bib23] (pp55-78).
In the continuation of this paragraph, we first of all will delay we on the theoretical framework of
operation of the algorithm. Then, before detailing the alternative installation in Code_Aster,
we will stick to the realistic framework of arithmetic finished.
5.2
Theoretical algorithm of Lanczos
5.2.1 Principle
Its perimeter of application covers the couples operator with work (pseudo) produced scalar ensuring
hermiticity of
With
. It makes it possible to build a matrix of Rayleigh repeatedly
B
m
of size
flexible, tridiagonale and square (with a true scalar product, if not it loses this
last property). This particular form largely facilitates the calculation of the modes of Ritz: with
QR
one loses a command busy magnitude then, when one seeks
p
clean modes while projecting on one
subspace of size
m
, of
O (pm
2
)
with
O (20pm)
[bib6].
The algorithm consists in gradually building a family of vectors of Lanczos
Q
1
,
Q
2
,…,
Q
m
while projecting orthogonally, with the iteration K, the vector
With
Q
K
on the two preceding vectors
Q
K
and
Q
k-1
. The new vector becomes
Q
k+1
and thus, gradually, one ensures structurally
the orthogonality of this family of vectors. The iterative process is summarized thus.
(
)
Q
0
Q
Q
Z A Q
Q
Z Q
v Z
Q
v
Q
v
0
0
1
1
0
1
0
=
=
=
,
.
/
.
,
,
,
;
Calculation of
For K = 1, m to make
=
-
,
=,
= -
,
=
So then
=
If not
Deflation
End if.
Fine loops.
K
k-1
K - 1
K
K
K
K
K
K
K +1
K
Algorithm 10: Theoretical Lanczos
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
37/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
While noting,
E
m
m
ième
vector of the canonical base, the vector residue of the factorization of
Lanczos is written
R
Q
E
m
m
m
m
T
=
+
1
.
Appear 5.2.1-a: Factorization of Lanczos
The matrix of Rayleigh takes the form then (in complex with a square scalar product)
B
m
m
m
m
=
-
-
1
1
1
2
1
1
0
0
0
0
0
0
…
…
…
.
By construction, this algorithm ensures us of the following results (into arithmetic exact):
·
(
Q
K
)
k=1, m
constitute an orthonormal family,
·
they generate the space of Krylov of command
m
initiate by
Q
1
,
(
)
(
)
K
Vect
m
m
With Q
Q With Q
With
Q
,
,
,…,
1
1
1
1 1
=
-
,
·
they make it possible to build a matrix
B
m
tridiagonale, square and of flexible size,
·
the spectrum of
B
m
does not admit that the m dominant simple modes on which
Q
1
one has
component.
Note:
·
the algorithm thus does not allow, in theory, the capture of multiple modes. That can
to be explained by noticing that any matrix of irreducible Hessenberg (a matrix
With
is
said of Hessenberg higher (resp. lower) if
With
I, J
= 0
for
i>j+1
(resp.
j>i+1
). It is
of more irreducible if
With
I
I
I
+
1
0
,
) (thus a fortiori a matrix tridiagonale) does not admit that
simple modes,
·
the cost of an iteration is weak, about that of a product matrix-vector
With
Q
K
, that is to say
for
m
first iterations a calculation complexity of
O (Nm (c+5))
(with C the number
means of nonnull terms on the lines of the matrix of work). In addition, complexity
necessary memory is weak bus one does not need to build
With
in entirety, one has right need
to know its action on a vector. This characteristic is very interesting for
to solve problems of big sizes [bib23],
·
in general the sphere of activity directs the choice of a vector of initial Lanczos, the aforementioned must by
example to belong to a space of acceptable solutions (checking stresses,
conditions limit…) like with the image unit of the operator. This last point is important
because it makes it possible more quickly to enrich modal search while not being limited to the core,
·
the algorithm of Lanczos is only one means of approaching the subspaces of Krylov. Its
generalization with the nonsymmetrical configurations is called algorithm of Arnoldi (cf [§6.2]).
With
Q
m
=
N
N
m
Q
m
B
m
+
m
m
0
R
m
=
m
Q
m+1
E
MT
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
38/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5.2.2 Estimates of errors and convergences
Because of particular form of
B
m
, the extraction of the required elements of Ritz is simplified, and in
more, the following result allows us quickly (via a product of two realities!) to estimate theirs
qualities.
Property 7
The euclidian norm of the residue of the element of Ritz
(
)
~, ~
U Q X
=
m
is equal to
(
)
R
With
I U
E X
2
2
=
-
=
~
~
m
m
T
Proof:
Commonplace by taking the euclidian norm of the factorization of Lanczos
WITH Q X
Q B X
Q
E X
m
m
m
m
m
m
T
=
+
+
1
where
Q
m+1
is normalized with the unit.
Note:
A more general result, independent of any method of resolution, us ensures that for
each clean mode
()
~,
X
of
B
m
, there is a clean mode
()
, U
of
With
such as:
()
-
=
-
~
sin, ~
min
~
~
~
R
U U
R
U With U
2
2
2
with
I
I
T
Therefore, more than minimization of the residue, the criterion of stop of the methods of the Lanczos type/Arnoldi
m
<
allows to approach the original spectrum as well as possible. These increases are not accessible
that in the square case.
In the general case (and in particular nonnormal) they are more difficult even impossible with
to build without additional information (level of nonnormality, “defectivity”.).
The estimate of the residue is not then any more the good criterion to estimate the quality of the approximated modes.
Let us interest being maintained in the quality of convergence of the algorithm. Let us particularize theorem 6
for this algorithm. By limiting the standard of error of
(
)
I P
U
-
m
2
increases are obtained
following:
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
39/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Theorem 8
That is to say
(
)
I
I
, U
I
ème
clean mode dominating of
With
diagonalisable, there exists a mode of Ritz then
()
~, ~
I
I
U
such as:
()
(
)
()
I
I
J I
J
N
J
I
m I
I
I
I
J I
J
N
J
I
m I
I
T
T
-
-
-
-
-
<
-
-
<
-
-
~
sin
, ~
2
2
2
1
U U
with
I
= 1 + 2
-
-
I
i+1
1+i
N
,
T
semi
(
X
)
the polynomial of Tchebycheff of degree
semi
,
the constant of
theorem 6 and
(
)
=
P U
Q U
m I
I
2
1
tan
,
.
Proof:
Cf [bib11] pp554-557.
These increases indicate (cf [bib23] pp59-63) that:
·
if the initial vector does not have any contribution along the clean vectors, one cannot capture
the aforementioned modes (
+
),
·
one
has
firstly convergence of the extreme modes and as much better than the spectrum
is separate (these methods are less sensitive to the clusters than the methods of the types
reiterated powers), without packing of eigenvalues (famous “the clusters”),
·
the error decrease when
m
increase (like the reverse of the polynomial
T
semi
(
X
)
, therefore like
the reverse of exponential),
·
the estimate on the eigenvalues is better than that their associated own vectors.
Note:
The convergence of the method was studied by P. Kaniel then improved (and corrected) by
S.C. Paige; One will find a synthesis of these studies in the paper of Y. Saad [bib20].
Let us look at now how behaves the algorithm into arithmetic finished. We will see that it
is the seat of surprising phenomena.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
40/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5.3
Practical algorithm of Lanczos
5.3.1 Problem
of orthogonality
At the time of implementation the effective of this algorithm, one realizes that very quickly the vectors of
Lanczos lose their orthogonalities. The matrix of Rayleigh is not then any more the matrix projected of
the initial operator and the captured spectrum are sullied more and more with errors. A long time this phenomenon
inescapable was allotted wrongly to the effects of rounded the arithmetic one finished. In 1980, D.C. Paige showed
that the loss of orthogonality was especially ascribable with the convergence of a clean mode. Thus, as of
that one captures a mode dominating, one disturbs the layout of the preceding vectors of Lanczos.
following result expresses this paradox clearly.
Property 9 (Analysis of Paige)
By taking again the notations of algorithm 10 and while noting
the precision machine, with the iteration
k+1
the orthogonality of the new vector of Lanczos is expressed in the form:
(
)
Q
Q
v Q
With
K
T
I
T I
K
I K
+
=
+
1
2
Proof:
Cf [bib26].
The problem occurs in fact during the standardization of the new vector of lanczos
Q
k+1.
When it
mode converges, according to property 7, that comes from the smallness of the coefficient
K
and thus in spite of
theoretical orthogonality (
(
)
v Q
T I
I
K
=
0
), effective orthogonality is far from being acquired
(
(
)
Q
Q
K
T
I
I
K
+
>>
1
).
The digital processing of these problems has been the search object many for thirty years
and of many palliative strategies were developed. The choice of such or such method depends
type of required spectral information and average data processing available, the synthesis of
Besides JL.Vaudescal proposes a very good overflight of these alternatives (cf [bib23] pp66-78):
·
Algorithm of Lanczos without reorthogonalisation, developed by J.K. Cullum and
R.A. Willoughby [bib3] who consists with expurger the calculated spectrum of his phantom modes of
studying interlacings of the eigenvalues of the matrix of work and one of its
submatrices. This alternative admits a weak overcost calculation and memory to determine them
eigenvalues, but it complexes (sometimes largely) the search of the vectors
clean.
·
Algorithms of Lanczos with total reorthogonalisation (B.N. Parlett) or selective
(J. Scott) which consists with each pitch with réorthogonaliser the last vector of Lanczos obtained
compared to all the already calculated vectors or simply compared to the vectors of Ritz
converged (this alternative thus makes it possible to control the loss of orthogonality dynamically
acceptable). These alternatives are much more expensive in complexity calculation (them
automatic réorthogonalisations) and memories (storage of the preceding vectors) but
they prove more effective and more robust.
In the alternative of Newmann-Pipano [bib14] (
METHOD = “TRI_DIAG”
) of Code_Aster, it is
strategy of reorthogonalisation supplements which was selected: the vector
Q
k+1
is thus not completely
calculated as in the theoretical algorithm.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
41/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Note:
·
these problems of loss of orthogonality occur with all the more of acuity that the size of
subspace increases. It is an additional argument for the installation of one
iterative process limiting this size (cf [§5.4.2]),
·
these strategies are declined vectoriellement or per blocks, they are established in
simple or iterative algorithms. The explicit latter being able to profit from restarts or
implicit, controlled or automatic. Alternative IRAM profits thus from an implicit restart
calculated automatically,
·
the central point of these algorithms is consisted the method of orthogonalization put in
place. All the process is dependant on its success and its robustness. With
orthogonal of type Housholder or Givens, very expensive transformations but very
robust, one prefers from now on the algorithms of the type Gram-Schmidt Itératifs (IGSM) which are
a better compromise between effectiveness and complexity calculation (cf [Appendix 2]). It is it besides
choice which was made for the alternatives of Lanczos/Arnoldi installation in the code,
·
the alternative of selective reorthogonalisation compared to the converged modes returns to
to carry out an implicit deflation by restriction on the operator of work not to have with
to recompute (cf [§3.1]).
5.3.2 Capture multiplicities
It was seen that in theory the algorithm of Lanczos did not allow, some is his strategy of
reorthogonalisation, the theoretical capture of multiple modes. In practice, for once, effects
from flare come to the rescue and “powder” with small components along almost
all clean vectors. One can thus capture multiplicities, however they can be
erroneous and require a complementary postprocessing of checking.
Note:
In fact, only a version per blocks (G. Golub & R. Underwood, 1979) can allow us one
correct detection of the multiplicities, at least as long as the size of the blocks is sufficient. This
version was wide (Mr. Sadkane [bib22], 1993) with the algorithm of Arnoldi but it would be a pity
to deprive itself of the alternative of Sorensen (IRAM) who is more effective and more robust.
5.3.3 Phenomenon of Lanczos
The success of this algorithm rested at the beginning on what is called the “phenomenon of Lanczos”.
This conjecture predicts that for a sufficiently large size of subspace (
m >> N
) one is
able to detect all the spectrum (with load to thereafter distinguish the “good grain from the ryegrass”) of
the operator of work. Taking into account weak the pre-necessary report of “basic” Lanczos
(grosso-modo a matrix tridiagonale and some vectors), this is particularly interesting for
to treat hollow systems of very big sizes (other algorithms of reiterated the powers type or
QR
they require it knowledge of all the matrix of work).
We now will detail (a little) some elements whose interest largely exceeds the framework
of this algorithm.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
42/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5.4 Processing
complementary
5.4.1 Detection of spaces invariants
In algorithm 10, the nullity of the coefficient
l-1
prevent the standardization of the new vector and
an additional processing called of deflation requires. The subspace of calculated Krylov is then under
space invariant and its spectrum is thus exact. In other words, them
l-1
first eigenvalues
exhumed by the algorithm of support (
QR
or
QL
) are those of
With
.
To continue, it is then necessary to choose a new random vector observing the limiting conditions,
the orthogonaliser with all the already calculated vectors of Lanczos and to build a second family of
vectors of Lanczos that one orthogonalise compared to the first family constituting space:
(
)
(
)
H
K
Vect
L
L
L
=
=
-
With Q
Q With Q
With
Q
,
,
,…,
1
1
1
1 1
.
The matrix tridiagonale obtained with the following form:
B
m
L
L
L
L
m
m
m
=
+
-
-
1
1
1
2
1
1
1
0
0
…
…
…
…
…
…
0000
0000
The intermediate factorization of Lanczos (to the command
L
) is written:
WITH Q
Q B
L
L
L
=
Values
(
K
,
K
)
k=1, L
are obtained starting from the first random vector and the values
(
K
,
K
)
k=l+1, m
are obtained starting from a second random vector. To detect the nullity of the term
extra-diagonal the numerical criterion is used
L
L
PREC LANCZOS
-
1
_
then
L
-
=
1
0
,
where
PREC_LANCZOS
is initialized under the key word factor
CALC_FREQ
.
Note:
·
a more robust criterion retained by large mathematical bookstores EISPACK [bib26],
LAPACK [bib26] and ARPACK [bib27] use the preceding diagonal term, a parameter
flexible
C
and the precision machine
,
(
)
L
L
L
C
-
-
<
+
1
1
it is this criterion which at summer appointed for IRAM (but the user cannot modify the parameter
C
,
the aforementioned is fixed at a standard value by the code), various messages preventing besides
the user of the detection of a space invariant.
·
Obtaining such a space is completely sympathetic nerve since it ensures us of very
good quality of the first modes. Moreover, one reduces the size of the problem by dividing it into
two parts: one solved, the other to solve. Many techniques seek with
to reproduce artificially this phenomenon (cf [§3.1], [§5.3.1]).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
43/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5.4.2 Strategies of restartings
To limit the endemic problems of orthogonalization, to limit to them pre-necessary report and
to take dynamically into account spectral information already obtained, of the strategies of
restartings were coupled with the algorithm of Lanczos. It loses its “simple” character and becomes
“iterative”. One reiterates the “restarts” until satisfying the desired criteria of convergence. The algorithm
does not undergo the convergence of the modes imposed by theorem 8 and it can support intéractivement
that of certain modes.
In theory however, plus the size of the subspace is large, better is convergence. One
compromise is thus to find between the size of the subspace and the frequency of the restartings.
Various vectors of restartings can be used, being generally written like a sum
balanced
p
required clean modes (Y. Saad 1980)
Q
Q X
1
1
=
=
I
m I
I
p
.
Indeed, if one starts again with an initial vector pertaining to the subspace invariant generated by
sought modes, one is then sure to obtain (with a good algorithm of orthogonalization) a standard
residual about the precision machine and of the almost exact clean modes. In addition to the decision
to start the restarting, all the problems lie in the search for these weights
I
. Us
will see that with IRAM, the restarts set up via implicit polynomial filters solve
elegantly and automatically this question.
Note:
·
this philosophy of the restarts was initiated by W. Karush (1951),
·
restartings based on the vectors of Schur associated with the spectral decomposition
sought seem more stable (J. Scott 1993),
·
criteria of effectiveness were required to decide restart appropriateness (B. Vital
1990),
·
instead of a simple linear combination of clean vectors, one can determine a vector
of restarting via polynomials (Tchebycheff in reality and Faber in complex) allowing
to center modal search in such or such area. One speaks then about polynomial accelerations
explicit [bib22].
The following paragraph will show some of the concepts presented hitherto to clarify
alternative installation in the code.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
44/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5.5
Establishment in Code_Aster
5.5.1 Alternative of Newmann & Pipano
This alternative developed by Mr. NEWMAN & A. PIPANO [bib3], [bib14], in 1977, is a method of
Simple Lanczos, into arithmetic real, with total reorthogonalisation (via a GSM). It uses it
crossbred “shift and invert” conventional of the scalar pseudo-product introduced by the shifted matrix (
With
-
B)
.
(
)
()
(
)
()
()
With
B
B
With
U
U
X y
y With
B X
X X
X
X X
-
= -
=
-
=
=
-
µ
1
1
1
2
!
“
# #
$
# #
!“$
with
,
,
,
,
,
.
T
sign
This choice returns the couple operator of work - symmetrical scalar product and it are adapted to the matrices
particular of the mechanics of the structures. One can thus deal with problems:
·
of free structure and fluid structure (
With
can be singular),
·
of buckling (
B
is indefinite).
The scalar pseudo-product introduced by the shifted matrix is thus regular (
answers this
prerogative cf [§2.6], [§2.9]) and it makes it possible to seek vectors of Lanczos (
With
-
B)
- orthogonal.
The total gonalisation is carried out only if it proves to be necessary and this criterion is skeletal.
Note:
·
B
- orthogonality and a fortiori that with the Euclidean direction cannot be but the fruit any more of
particular configurations. This does not modify the properties of orthogonalities of the modes
clean (cf proposal 2),
·
as one already announced [§1.7] there is a whole zoology of couples operator - product
scalar, the aforementioned being only one possibility among others. Thus, the libraries [bib29]
conventional propose to deal with the problems of buckling in “buckling mode” via the same one
“shift and invert” and the scalar pseudo-product introduced by
With
. But because of introduction
quasi-systematic of Lagranges, this matrix becomes indefinite even singular, which
disturb the process largely. The same causes produce the same effects when
for a calculation of dynamics, one uses it
B
- scalar product.
The price to be paid for this gain in robustness is the loss of possible symmetry of the matrix of Rayleigh
(
)
B
Q
Q
m
m
m
m
I
I I
I
I
I
sign
=
=
=
-
-
+
+
1
1
1
2
1
1
1
1
0
0
0
0
0
0
…
…
…
,
with
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
45/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
After balancebeing balanced and putbeing put in the form of Hessenberg higher,
B
m
is diagonalisée by one
method
QL
implicit if it remains in spite of very symmetrical or by a method
QR
if not
(cf [Appendix 1]). The difference in cost calculation between these solveurs remains negligible vis-a-vis the costs of
réorthogonalisations. The diagram of established Lanczos becomes then:
(
)
()
()
(
)
()
()
(
)
(
)
(
)
(
)
(
)
Random hard copy of
For K = 1, m to make
=
-
,
=
= -
,
Réothogonalisation from report/ratio with
K
k-1
k-1
k-1
K
K
K
K
Q
Q
With
B
Q
U
Q
With Q
U
Q With
B Q
Q
Q Q A
B Q
Q
0
Z A Q
Q
Q With
B Z
v Z
Q
v
Q
init
init
lagr
I
bloq
I
T
T
K
T
I
I
I
I
sign
.
~
.
.
%
~
.
.
%
%
% %
%
,
,
.
,
1
1
1
1
1
1
1
1
1
1 1
1
1
2
0
0
0
0
0
=
-
=
=
-
=
-
=
=
=
-
-
-
()
(
)
(
)
(
)
(
)
I
K
K
K +1
K
IGSM
=
So then
=
If not
Deflation
End if.
Fine loops.
I
K
T
K
T
K
sign
=
-
=
-
1
1
2
0
,
,
,
,
,
;
v A
B v
v A
B v
Q
v
Algorithm 11: Alternative of Newman-Pipano
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
46/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
The reorthogonalisation is carried out thanks to an alternative “house” of the Method of Gram-Schmidt
Iterative (IGSM) according to the following process:
For
I = 1, K
(
)
has
K
T
I
=
-
+
Q
With
B Q
1
,
If
(
)
(
)
Q
With
B Q
K
T
I
PREC ORTHO
+
-
1
_
then
For
J = 1,
NMAX_ITER_ORTHO
(
)
(
)
(
)
X Q
Q
With
B Q Q
X
With
B Q
=
-
-
=
-
+
+
K
K
T
I
I
T
I
B
1
1
,
If
(
)
B
PREC ORTHO
_
then
Q
X
K
+
=
1
,
I
I
+
1
, exit;
If not
If
(
)
B has
then
Q
X
K
+
=
1
,
B has
=
;
If not
Failure, emission of a message of alarm,
I
I
+
1
, exit;
End if.
End if.
Fine loops in
J
.
End if.
Fine loops in
I
.
Algorithm 12: Procedure of reorthogonalisation of
“TRI_DIAG”
Note:
After some tests, it seems that this alternative specific to the code is less effective than
the IGSM of the type Kahan-Parlett [bib18] selected for the IRAM (cf [Appendix 2]).
5.5.2 Parameter setting
To be able to activate this method, the key word should be initialized
METHOD
with
“TRI_DIAG”
. In addition,
the size of the subspace of projection is determined, either by the user, or empirically to leave
formula:
(
)
(
)
m
p p
N
ddl active
=
+
-
min max
,
,
4
7
where
p
is the number of eigenvalues to calculate,
N
ddl-credits
is the number of active degrees of freedom of the system (cf [§2.2])
The user can always impose to him even dimension by indicating it with the key word
DIM_SOUS_ESPACE
key word factor
CALC_FREQ
.
Parameters of the total reorthogonalisation,
PREC_ORTHO
and
NMAX_ITER_ORTHO
, the number
maximum of iterations
QR,
NMAX_ITER_QR
, and the criterion of deflation (cf [§5.4.1]),
PREC_LANCZOS
, are
accessible by the user under the key word factor
CALC_FREQ
. When this phenomenon of deflation
product, a specific message specifies its row
L
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
47/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5.5.3 Warning on the quality of the modes
Let us recall that this alternative is not iterative. From the number of desired frequencies, one
estimate the dimension of the subspace of calculation and one hopes that the clean modes of the problem
projected standard will be a good approximation of those of the generalized problem. One checks has
posteriori the validity of the results (cf [§2.9]) but one does not have any active control on this precision.
If they are not satisfactory, one has as only only resorts to carry out a new test in
increasing the dimension of the subspace.
Note:
·
the method WILL GO that we will approach proposes a flexible and dynamic control of
precision of the results, it is an iterative method,
·
in V5, in parallel of the modal position and eigenvalue (in frequency if one is in
dynamics), one displays from now on the standard of the residue of the initial problem calculated in
postprocessing (after having standardized the clean vectors ad infinitum (in order to exhume
clean vectors “neutral” with respect to the string of standardizations which they can undergo at the time
calculations of the modal parameters of the structure (factors of participation…)).
5.5.4 Perimeter
of use
Attention, this method is usable one to the only deal with quadratic problem. In the event of error
at the time of the initialization of the key word
METHOD
, calculation stops and a trace is left in the files of
exit.
In V5, an algebraic method calculating the null eigenvalues of the modal problem
generalized was introduced. Physically those correspond to the movements with energy of
null deformation of a free structure (cf [bib25] TP n°2). Except the numerical difficulties of
handling of quasi-null quantities, because of their multiplicities, their capture correct was until
often problematic present for Lanczos established in Code_Aster: phantom modes
appeared corresponding to missed multiplicities!
This algorithm of detection of the modes of rigid body, which one activates while initializing
OPTION
with
“MODE_RIGIDE”
(default value `
WITHOUT
`), intervenes in preprocessing of modal calculation properly
known as. It is based on the analysis of the matrix of rigidity and breaks up into three phases:
·
detection of the null pivots of this matrix,
·
blocking of these pivots,
·
resolution of a linear system whose are solutions the associated clean vectors.
During the process of Lanczos it is enough consequently to orthogonaliser, progressively with their
determination, basic vectors with the latter.
However, the introduction of the method WILL GO reduced the interest (if it is not as comparison) of one
such option (which is not free in calculation complexity since it requires inversions of systems).
What good was it to deprive itself of such an algorithm which is by no means affected by the presence of these modes
a little particular multiples!
This method remains nevertheless useful as solution of help in the event of failure of IRAM (that
can occur in the presence of clusters in the vicinity of the terminals or for operators
pathological (badly conditioned, nonnormal, defective…) who should however be rather rare
(that often reveals a badly posed problem). One can also plan to encapsulate it in one
auxiliary operator (like
IMPR_STURM
).
Note:
The quadratic processing of problem is incompatible with the options
“MODE_FLAMB”
and
“MODE_RIGIDE”
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
48/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5.5.5 Display in the file message
The example below resulting from the list of case tests of the code (sdll112a) recapitulates the whole of the traces
managed by the algorithm. The iteration count
QR
effective (or
QL
) can be only identical for
all eigenvalues. It is an artefact of information which will be brought to disappear.
------------------------------------------------------------------------
THE NUMBER OF DDL
TOTAL EAST: 86
LAGRANGE EAST: 20
THE NUMBER OF ACTIVE DDL EAST: 56
------------------------------------------------------------------------
The SELECTED OPTION EAST:
PLUS_PETITE
THE VALUE OF SHIFT IN FREQUENCY EAST: 0.00000E+00
------------------------------------------------------------------------
INFORMATION ON CALCULATION REQUIRES:
A NUMBER OF REQUESTS MODES: 10
THE DIMENSION OF REDUCED SPACE EAST: 0
IT IS LOWER THAN THE NUMBER OF MODES, ONE TAKES IT EQUALIZES A 40
------------------------------------------------------------------------
FREQUENCIES CALCULEES INF. AND SUP. ARE:
FREQ_INF: 1.54569E+01
FREQ_SUP: 1.01614E+02
THE FIRST HIGHER FREQUENCY NOT SELECTED EAST: 1.29375E+02
------------------------------------------------------------------------
MODAL CALCULATION: METHOD Of SIMULTANEOUS ITERATION
METHOD OF LANCZOS
NUMBER FREQUENCY (HZ) STANDARD Of ERROR ITER_QR
1 1.54569E+01 1.29798E-12 4
2 1.54569E+01 7.15318E-13 4
3 3.35823E+01 3.98618E-13 4
4 3.35823E+01 4.25490E-12 4
5 4.73076E+01 4.26463E-12 4
6 4.73076E+01 1.43391E-12 4
7 5.45850E+01 9.52006E-12 4
8 8.80156E+01 3.45489E-13 4
9 1.01614E+02 6.13949E-12 4
10 1.01614E+02 3.18663E-12 4
------------------------------------------------------------------------
CHECKING A POSTERIORI OF THE MODES
IN the INTERVAL (1.54182E+01, 1.16326E+02)
IT THERE A WELL 10 FREQUENCY (S)
------------------------------------------------------------------------
Example 5:
MODE_ITER_SIMULT
with
“TRI_DIAG”
Now let us recapitulate the parameter setting available of the operator
MODE_ITER_SIMULT
with this
option
METHOD = “TRI_DIAG”
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
49/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
5.5.6 Summary of the parameter setting
Key word factor
Key word
Default value
References
“DYNAMIC” TYPE_RESU
“DYNAMIC”
[§2.1]
“MODE_FLAMB”
[§2.1]
METHOD “TRI_DIAG”
“SORENSEN”
[§5.5.3]
OPTION “MODE_RIGIDE”
“WITHOUT”
[§5.5.4]
“WITHOUT”
[§5.5.4]
CALC_FREQ
FREQ
[§4.4]
CHAR_CRIT
[§4. 4]
OPTION “PLUS_PETITE”
`PLUS_PETITE ''
[§4.4]
“TAPE”
[§4.4]
“CENTER”
[§4.4]
NMAX_FREQ
10
[§4.4]
DIM_SOUS_ESPACE
Calculated
[§5.5.2]
PREC_ORTHO
1.E-12
[§5.5.1], [§5.5.2]
NMAX_ITER_ORTHO
1.E-04
[§5.5.1], [§5.5.2]
PREC_LANCZOS
1.E-04
[§5.5.1], [§5.5.2]
NMAX_ITER_QR
15
[§5.5.1], [§5.5.2]
NPREC_SOLVEUR
8
[§2.6]
NMAX_ITER_SHIFT
5
[§2.6]
PREC_SHIFT
0.05
[§2.6]
SEUIL_FREQ
1.E-02
[§2.9]
VERI_MODE
STOP_ERREUR “YES”
“YES”
[§2.9]
“NOT”
[§2.9]
PREC_SHIFT
5.E-03
[§2.9]
THRESHOLD
1.E-06
[§2.9]
STURM “YES”
“YES”
[§2.9]
“NOT”
[§2.9]
Summary table 5.5.6-a: of the parameter setting of
MODE_ITER_SIMULT
with
“TRI_DIAG”
Note:
·
in V5, the user can specify the class of membership of his calculation by initializing it
key word
TYPE_RESU
. According to this value, one informs the vector
FREQ
or
CHAR_CRIT
,
·
one finds all the “tripaille” of parameters related to the preprocessings of the test of Sturm
(
NPREC_SOLVEUR
,
NMAX_ITER_SHIFT
,
PREC_SHIFT
) and with postprocessings of checking
(
SEUIL_FREQ
,
VERI_MODE
),
·
at the time of the first passages, it is strongly advised to modify only them
main parameters noted in fat. The others relate to more the mysteries of
the algorithm and they were initialized empirically with values standards,
·
in particular, to improve quality of a mode, the only flexible parameter is
dimension of the subspace,
DIM_SOUS_ESPACE
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
50/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
6
Algorithm WILL GO (
METHOD = “SORENSEN”
)
6.1 Introduction
We saw that one of the crucial problems of the method of Lanczos is the loss of orthogonality
inescapable of its basic vectors (cf [§5.3.1]). A generalization of this algorithm to the case not
square imagined by W.E. Arnoldi [bib31] in 1951 makes it possible to solve partially this
problems. It was given to the style of the day by Y. Saad [bib32] in 1980 and one abundant
literature covers the subject. Put aside the article founder and papers of Y. Saad and Mr. Sadkane [bib22],
one recommends the brought up to date and exhaustive synthesis of J.L. Vaudescal [bib23] (pp79-112).
This method being at the base of algorithm IRAM known as “of Sorensen” (IRAM for `Implicit
Restarted Arnoldi Method'), we first of all will detail his operation, its
behaviors and its limitations. Thereafter, we will determine the stakes to which IRAM must answer
(and it does it in the majority of the standard cases!) and we will consider its theoretical mysteries
and numerical. We will conclude by the summary from its effective parameter setting in Code_Aster and
for an example of file message.
6.2 Algorithm
of Arnoldi
6.2.1 Principle
Its perimeter of application covers all the couples operator with work (pseudo) produced
scalar. During this opening is the filling of the matrix of Rayleigh which becomes
form Hessenberg higher. It is not very prejudicial, because one can thus apply to him
directly the algorithm
QR
(the first stage of a good
QR
, except balancing, consists in reducing
the matrix of work in the form of Hessenberg. That makes it possible to gain a command magnitude at the time of
the resolution itself (cf [Appendix 1]), from where a gain about
O (10m
3
/3)
.
The algorithm is very similar to that of Lanczos, it consists in building a family gradually
of vector of Arnoldi
Q
1
,
Q
2
,…,
Q
m
while projecting orthogonally, with the iteration
K
, the vector
With
Q
K
on
K
vectors precedent. The new vector becomes
Q
k+1
and thus, gradually, one ensures
the orthogonality of this family of vectors.
With the difference of Lanczos, the orthogonality of the new vector compared to all the precedents is
thus explicitly and not assured implicitly. The aforementioned is managed by the algorithm of Gram-Schmidt
Modified (GSM) (cf appendix 2) which proves sufficiently robust in the majority of the cases.
While noting,
E
m
m
ième
vector of the canonical base, the vector residue of the factorization of Arnoldi
is written
R
Q
E
m
m
m
m
m
T
B
=
+
+
1
1
,
.
Appear 6.2.1-a: Factorization of Arnoldi
N
m
With
Q
m
=
N
Q
m
B
m
+
m
m
0
R
m
=
B
m+1, m
Q
m+1
E
MT
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
51/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
The iterative process is summarized as follows:
()
Calculation of
For K = 1, m to make
=
,
For L = 1, K to make
B =,
B
Fine loops;
B
,
If B
then
= B
If not
Deflation
End if.
Fine loops.
K
lk
L
lk
L
K 1, K
K 1, K
K +1
K 1, K
Q
Q
Z A Q
Z Q
Z Z
Q
Z
Q
v
1
1
1
0
/
.
,
,
,
;
=
= -
=
+
+
+
(GSM)
Algorithm 13: Theoretical Arnoldi
The matrix of Rayleigh is written then:
B
m
m
m
m
m
m m
mm
B
B
B
B
B
B
B
B
B
=
-
-
11
12
1
21
22
2
1
1
0
0
0
…
…
…
…
,
,
.
Except the shape of this matrix and a less acuity with the problems of orthogonality, this algorithm
us ensures the same theoretical and numerical properties that Lanczos. However if one
let grow indefinitely the size of the subspace until convergence of the eigenvalues
wished, the effects of flare despite everything will take again the top and one goes to the front of large troubles.
From where need, as for Lanczos, to make this process iterative via restartings.
Note:
·
this method can be seen like an implicit alternative of the algorithm of Lanczos with
total reorthogonalisation (cf [§5.3.1]),
·
the implicit orthogonalization of the algorithm can be led by more expensive algorithms
but more robust such as
QR
or IGSM. This is strongly necessary when the operator of
work presents a too strong defect of normality,
·
because of structure of
B
m
, memory complexity is requested than with Lanczos, by
against complexity calculation remains of the same order of magnitude
O (Nm (c+3+m)
)
(with
C
numbers average nonnull terms on the lines of the matrix of work),
·
to improve this first point, Y. Saad [bib32] showed that the structure of Hessenberg
higher can see to cancel its extreme on-diagonals if one carries out only partially
reorthogonalisation,
·
the vectorial algorithms having a tendency natural to miss by the multiplicities, one prefers
often to use a version blocks (Mr. Sadkane [bib22], 1993). But the size of those influence
on the quality of the results, for this reason one prefers the vectorial versions to them or
blocks of IRAM,
·
the choice of the vector of initial Arnoldi is carried out same manner as for Lanczos.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
52/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
6.2.2 Estimates of errors and convergence
With regard to the evaluation of the quality of approximation of the clean modes obtained, there is one
as simple and effective criterion as for Lanczos.
Property 10
The euclidian norm of the residue of the element of Ritz
(
)
~, ~
U Q X
=
m
is equal to
(
)
R
With
I U
E X
2
2
1
=
-
=
+
~
~
,
B
m
m
m
T
Proof:
Commonplace by taking the euclidian norm of the factorization of Arnoldi
WITH Q X Q B X
Q
E X
m
m
m
m
m
m
m
T
B
=
+
+
+
1
1
,
and like
Q
m+1
is normalized with the unit.
Always by focusing us on the standard of the additional projector
(
)
I P
U
-
m
2
one can
to generalize the theorem of convergence 8 with the nonsquare case.
Theorem 11
That is to say
(
)
1
1
, U
the first (conventional arrangement, by order descending of module) clean mode
dominating of
With
diagonalisable and is
Q
U
1
1
=
=
K
K
K
N
the initial vector of Arnoldi broken up on
base clean vectors, it exists a mode of Ritz then
()
~, ~
1
1
U
such as:
1
1
2
1 1
-
~
m
with
=
P U
m I 2
,
the constant of theorem 6,
1
1
1
1
=
=
K
K
K
N
,
and
1
1
2
1
2
1
1
m
K
K
J
K
K J
m
J
m
=
-
-
=
+
=
+
-
,
. This result is declined in the same way on the other modes.
Proof:
By taking again to the demonstration of Y. Saad ([bib33] pp209-210) and the result of theorem 6.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
53/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
These increases very different from those obtained with Lanczos guide the same ones however
phenomena:
·
if the initial vector does not have any contribution along the sought clean vectors, one cannot
to capture (
I
+
),
·
one
has
firstly convergence of the peripheral modes of the spectrum, and this, of as much
better than it is separate (property of
im
),
·
the decrease of the error is proportional to the increase in
m
(property of
im
).
Note:
·
when an eigenvalue is badly conditioned
I
+
then it is necessary to increase
m
so that
im
decreases,
·
similar results were exhumed in the case of a defective operator (cf Zia 94 [bib23]).
Extremely from this lesson, we now will recapitulate the stakes to which must answer
IRAM.
6.3
stakes
The algorithm tent will bring an elegant remedy for the recurring problems raised by the others
approaches:
·
minimization of the space of projection: it proposes with minima
m > p+1
instead of
m= 4p
of Lanczos,
·
optimal management of the overcosts of orthogonalization establishing a compromise enters the size of
subspace and the frequency of the restartings,
·
transparent, dynamic and effective management of these restarts,
·
taking into automatic account of spectral information,
·
fixing
pre-necessary report and of the quality of the results.
There is thus more question to be posed concerning the strategy of reorthogonalisation, the frequency of
restarts, their establishments, criteria of detection of possible phantom modes… “super-IRAM”
charge of all!
In short, it gets:
·
one
better
total robustness,
·
calculations complexities
O (4nm (m-p))
and memories
O (2np+p
2
)
improved (especially by
report/ratio in simple Lanczos such that of Newmann & Pipano) for a fixed precision,
·
a more rigorous capture of the multiplicities, clusters and rigid modes of body
(thus less parasitic modes!).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
54/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Note:
·
on this last point, only a version per blocks of IRAM or a version incorporating of the “purging
and lock " (techniques of capture and filtering, cf D. Sorensen & R.B. Lehoucq, 1997) can
us to guarantee a correct detection of the spectrum of a standard operator (i.e not too badly
conditioned),
·
it seems, that in practice in Code_Aster, the report/ratio in calculation complexity between IRAM and
Lanczos/Bathe & Wilson are with minima of command 2 in favor of the first. With sizes of
reasonable problems (a few tens of thousands of ddls and
p= O (100)
) the aforementioned can
to go up up to 10 (without the encapsulation of
MACRO_MODE_MECA
). In certain cases
semi-industrial, it made it possible to unroll a search for spectrum which had failed with
Jacobi and who were inaccessible with Lanczos (taking into account the time limits),
·
a class of algorithm known as of “Jacobi-Davidson” (cf R.B. Morgan 1990) seems even more
promising to treat pathogenic cases. It uses an algorithm of the Lanczos type that it
packaged via a method of Rayleigh.
In the following paragraph we will clarify the operation of IRAM.
6.4
Algorithm “Implicit Restarted Arnoldi” (WILL GO)
This algorithm was initiated by D.C. Sorensen [bib30] in 1992 and makes real great strides for the resolution
great modal systems on parallel supercomputers. Its framework of application is very with
general fact. It deals with as well the real problems as complex, square or not. It is summarized in
a succession of factorizations of Arnoldi whose results control automatically
static and implicit restartings, via polynomial filters modelized by
QR
implicit.
First of all it carries out a factorization of Arnoldi of command
m= p+q
(in theory
Q = 2
is enough, in
practical
q= p
is preferable. It is this last default value besides which was retained
(cf [§6.5.3]) of the matrix of work. Then once this preprocessing carried out it reiterates a process of
filtering of the part of the undesirable spectrum (numerically and méthodologiquement, it is easier
to exclude to incorporate).
It starts by determining the spectrum of the matrix of Rayleigh (via the indétronable
QR
) and it in
dissociate the nondesired part (while referring to the criteria fixed by the user) which it uses then
like shift to set up a series of
Q
QR
implicit with simple shift (cf [Appendix 1]).
factorization is written then:
WITH Q
Q B
R Q
m
m
m
m
+
+
+
=
+
where
Q
Q Q
m
m
+
=
,
B
Q B Q
m
T m
+
=
and
Q Q Q
Q
=
1 2
.
Q
the unit matrix associated
QR
. Afterwards
to have updated the matrices, one truncates them until the command
p
WITH Q
Q B
R
p
p
p
p
+
+
+
+
=
+
and thus, at the price of
Q
new iterations of Arnoldi, one can find a factorization of Arnoldi of command
m
who is viable. All the subtlety of the process rests on this last sequence. Let us note:
()
(
)
With
With
Id
=
-
+
=
~
p I
I
Q
1
,
the matric polynomial of command
Q
generated by the operator of work. In fact, them
QR
implicit have acts
on
p
first lines of the matrix of Arnoldi, Rayleigh and the residue, so that
factorization complementary to command
Q
produces the same effect as a factorization of command
m
initiated by the vector
()
~
Q
With
Q
1
1
=
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
55/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
The initial vector was implicitly modified (it does not have construction and effective application there of
polynomial) so that it generates the desired modes preferentially and this, by withdrawing them
components considered to be “unsuitable”. As one already pointed out (cf [§5.4.2]) this type of restart
allows to decrease the residue by reducing the components of the initial vector following the modes
undesirable. Into arithmetic exact, one would have immediately
R
m
= 0
and the process would stop there!
Iteration after iteration, one thus improves quality of the modes sought while resting on
auxiliary modes. Of course the aforementioned is estimated at each stage and conditions the opportunity of
to simulate another restart or not. The process can be summarized (very macroscopically!) like
follows:
Factorization of Arnoldi of command m:
.
For K 1,
to make
To calculate
with implicit shifts,
Update
and
Truncation of these matrices has the command p,
,
Quality estimate of the p modes.
Fine loops.
1
2
p
p+1
p +q
AQ
Q B
R
QR
Q
B
R
AQ
Q B
R
AQ
Q B
R
m
m m
m
m
m
m
p
p p
p
m
m m
m
NMAX_ITER_SOREN
=
+
=
=
+
=
+
~, ~… ~, ~… ~,
,
,
Preserved for
améliorat ion
Used like
shifts
! “
#
$
#
! “
#
#
$
# #
Algorithm 14: Method WILL GO (known as of Sorensen)
The maximum number of iterations is controlled by the key word
NMAX_ITER_SOREN
key word factor
CALC_FREQ
.
It should be noticed that the Arnoldi algorithm is prudently supplemented by a reorthogonalisation
total (started that if that proves to be necessary). This overcost is all the more acceptable as it is
established via algorithm IGSM of Kahan-Parlett (cf [Appendix 2]) which is particularly effective.
All this makes it possible to ensure us of the good behavior of projection with respect to the initial spectrum.
In addition, the evaluation of the quality of the modes is not carried out simply by building them
p
residues
(
)
R
With
I U
I
I
I
2
2
=
-
~
~
via property 10. One already mentioned that in the case not
square, they could not suffir with this task, in particular in the event of strong defect of normality. For of
to discharge, without having recourse to other information (to obtain a criterion exact it would be necessary to be able
to consider conditionings spectral of spaces invariants and the angles which they make between them. It
who is sometimes difficult to obtain, even a posteriori!) a priori, the preceding property is used
supplemented by a criterion due to Z.Bai et al. [bib34]
(
)
B
PREC_SOREN
m
m
m
T
m
+
<
1,
max
,
~
E X
B
where
is the precision machine and
PREC_SOREN
is a key word initialized under
CALC_FREQ
. The user has
thus a quality control (partial) of the modes, that of which it did not lay out with the other methods
established in the code. Taking into account the various standards used, this error is different from
that resulting from the total postprocessing (cf [§2.9]) which is displayed in the columns results.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
56/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Note:
·
the technique of polynomial acceleration used is more effective than that of Tchebycheff,
since the latter is explicit and requires
m
products matrix-vector,
·
to avoid the deterioration of the vectors of Ritz (and thus of the approximate clean vectors) by
eigenvalues of very large modules (partners with the core of
B
in “shift and invert”)
a filtering of the type Ericsson & Ruhe [bib35] was established. More robust techniques
exist but they require information a priori concerning in particular the blocks of
Jordan associated with the core with the operator (cf Meerbergen & Spence, 1996),
·
this algorithm can be seen like the truncated shape of the algorithm
QR
implicit of
J.C.F. Francis (cf [Appendix 1]).
The following paragraph will clarify the choices which led to the alternative installation in the code.
6.5
Establishment in Code_Aster
6.5.1 ARPACK
The package of origin [bib29] (ARPACK for ARnoldi Package) coding the method is available in
freeware on Internet. These designers, D. Sorensen, R. Lehoucq and C.Yang of Rice University of
Houston, wanted it at the same time:
·
simple of access (FORTRAN77, “reverse communication”),
·
flexible (it is based on libraries LINPACK, LAPACK [27] and BLAS [bib36] (BLAS
is the acronym for BASIC Linear Algebra Subprograms)),
·
and rich in functionalities (decomposition of Schur, shifts flexible, many
spectral transformations).
Its effectiveness is multiplied by ten by the use of very optimized BLAS from level 2 and 3 and by the setting in
place “reverse communication”. The user is thus a Master of his structures of data and
of its procedures of processing concerning the operator and the scalar product of work. It is him which
provides to routines ARPACK this type of information. That thus makes it possible to manage as well as possible, with
tools and the procedures ASTER, the products matrix-vector, factorizations, resolutions of
system…
6.5.2 Adaptations of the algorithm of Sorensen
To deal with generalized modal problems, this package proposes a whole series of
spectral transformations, in reality or complex, square or not. Into square, the algorithm of
Sorensen is based on the Lanczos- couple
QL
(IRLM for Implicit Restarted Lanczos Method).
two approaches (square or not) are not designed besides to treat pseudo-products
scalars related to indefinite matrices.
Precisely, for at the same time circumscribing numerical problems involved in their properties enough
heterogeneous in the code and, in addition, to ensure itself of a better total robustness, we have
chosen to work in nonsymmetrical (IRAM with Arnoldi and
QR)
, on the couple operator of
following scalar work-product:
(
)
()
With
B
B
With
U
U
X y
y X
-
= -
=
-
µ
1
1
! “
# #
$
# #
!“$
,
T
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
57/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
One could have dealt with the problems of buckling in “buckling mode” via the same “shift and invert” and
the scalar pseudo-product introduced by
With
. But because of quasi-systematic introduction of
Lagranges, this matrix becomes indefinite even singular, which disturbs the process largely.
The same causes produce the same effects when, for a calculation of dynamics, one A resorts to
B
- scalar product. Rather than to modify all the package while introducing a scalar pseudo-product,
we thus chose a simple scalar product “Euclidean” more robust and much less
expensive.
It was thus necessary to modify the procedures of “reverse-communication” of the package, because it did not envisage
this option (with matrices standards one classically prefers to enrich the components with one
matric scalar product, even in nonsymmetrical). Contrary to the alternative of Newmann &
Pipano (cf [§5.5.1]) introduced for Lanczos, we deliberately placed ourselves in one
nonsymmetrical configuration. But in order to avoid as much as doing it can the problems of orthogonality
recurring, even into symmetrical, we would have chosen the version of IRAM using Arnoldi.
The disadvantage of this step is that it is necessary, in preliminary postprocessing of IRAM,
B
- orthonormaliser approximate clean vectors to find property 2 numerically
exploited by the modal recombinations. This stage does not disturb the base of clean modes
exhumed and it is very effectively carried out via the IGSM of Kahan-Parlett.
The taking into account of the limiting conditions and, in particular of the double dualisations, was led
as for Lanczos according to the procedure describes in [§2.2]. In particular, once the vector
initial in acceptable space one is applies the operator of work to him. This conventional process
allows to purge the vectors of Lanczos (and thus the vectors of Ritz) of the components of the core.
In addition, in certain configurations for which the number of frequencies requested
p
is equal to the number of ddls active (cf [§2.2]), one owed bluffer the algorithm which stopped in error
fatal! Indeed, it detected generally well space invariant awaited (from size
p
), but because of
particular structure of the clean vectors associated Lagranges (cf proof of property 4, [§2.5])
it had much evil to generate an initial vector which is proportional for them.
One would have needed particular processing taking account of the classification of these Lagranges, which
would have been all the more expensive as they are not fundamentally necessary to solve it
problem requested! All the spectral information being already present at deepest of the algorithm,
it is not thus the sorrow to complete the two remaining iterations (when the user asks
p
=
N
ddl-credits
, one imposes automatically
m= p+2
). It is enough to short-circuit the natural yarn of
the algorithm, to withdraw the interesting modes of Ritz and post-to treat them to return in space
initial.
Note:
This type of case of figure in which one seeks a number of clean modes very near to
ddls numbers leaves the “ideal” perimeter of use of this type of algorithm (cf [§2.8]). A good
QR
would be without any doubt more effective, but it is a good means of testing the algorithm.
6.5.3 Parameter setting
To be able to activate this method, the key word should be initialized
METHOD
with `
SORENSEN
'. Size of
subspace of projection is determined, either by the user, or empirically from
formulate:
(
)
(
)
m
p p
N
ddl active
=
+
-
min max
,
,
2
2
where
p
is the number of eigenvalues to calculate,
N
ddl-credits
is the number of active degrees of freedom of the system (cf [§2.2])
The user can always impose to him even dimension by indicating it with the key word
DIM_SOUS_ESPACE
key word factor
CALC_FREQ
.
The parameter of the IGSM of Kahan-Parlett (cf [Appendix 2])
PARA_ORTHO_SOREN
, the maximum number
iterations of the total process
,
NMAX_ITER_SOREN
, and the criterion of quality control of the modes,
PREC_SOREN
, are accessible by the user under the key word factor
CALC_FREQ
. When this last
key word is null, the algorithm initializes it with the precision machine (
= 2.22.E-16 on SGI).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
58/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
6.5.4 Display in the file message
The example below resulting from the list of case tests of the code (ssll103b) recapitulates the whole of the traces
managed by the algorithm. One finds in particular, for each critical load (or frequency),
the estimate of its quality via the standard of error.
Here, IRAM reiterated only only once and used 30 IGSM (in its first phase). The resolution
total consumed 91 products (makes some, less than that because of “implicit” introduction of the product
Euclidean scalar) matrix-vector and 31 inversions of system (just increase because the operator of
work is already factorized).
------------------------------------------------------------
THE NUMBER OF DDL
TOTAL EAST: 68
LAGRANGE EAST: 14
THE NUMBER OF ACTIVE DDL EAST: 47
------------------------------------------------------------
The SELECTED OPTION EAST:
PLUS_PETITE
THE VALUE OF SHIFT CRITICAL LOAD EAST: 0.00000E+00
------------------------------------------------------------
INFORMATION ON CALCULATION REQUIRES:
A NUMBER OF REQUESTS MODES: 10
THE DIMENSION OF REDUCED SPACE EAST: 30
=============================================
= METHOD OF SORENSEN (CODE ARPACK) =
= VERSION: 2.4 =
= DATE: 07/31/96 =
=============================================
A NUMBER OF RESTARTINGS = 1
A NUMBER OF PRODUCTS COP * X = 31
A NUMBER OF PRODUCTS B * X = 91
A NUMBER OF REORTHOGONALISATIONS (STAGE 1) = 30
A NUMBER OF REORTHOGONALISATIONS (STAGE 2) = 0
A NUMBER OF RESTARTINGS OF A NULL V0 = 0
------------------------------------------------------------
CRITICAL LOADS CALCULEES INF. AND SUP. ARE:
CHARGE_CRITIQUE_INF: - 9.96796E+06
CHARGE_CRITIQUE_SUP: - 6.80007E+05
------------------------------------------------------------
MODAL CALCULATION: METHOD Of SIMULTANEOUS ITERATION
METHOD OF SORENSEN
NUMBER CRITICAL LOAD NORMALIZES ERROR
1 - 6.80007E+05 5.88791E-12
2 - 7.04572E+05 1.53647E-12
3 - 7.09004E+05 1.16735E-12
4 - 7.10527E+05 1.72306E-12
5 - 7.11205E+05 2.41783E-12
6 - 7.11542E+05 7.88981E-13
7 - 7.11703E+05 5.71621E-13
8 - 1.50492E+06 1.17776E-11
9 - 6.02258E+06 2.42221E-11
10 - 9.96796E+06 3.55014E-12
------------------------------------------------------------
CHECKING A POSTERIORI OF THE MODES
IN the INTERVAL (- 1.00178E+07, - 6.76607E+05)
IT THERE A WELL 10 LOAD (S) CRITICAL (S)
------------------------------------------------------------
Example 6:
MODE_ITER_SIMULT
with `
SORENSEN'
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
59/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Note:
The introduction of this method made it possible to balance many cards of faults related to
multiplicities, of the clusters or search of eigenvalues of orders of magnitude very
different on which Lanczos and Bathe & Wilson stumbled.
Now let us recapitulate the parameter setting available of the operator
MODE_ITER_SIMULT
with this
option
METHOD = “SORENSEN”
.
6.5.5 Summary of the parameter setting
Key word factor
Key word
Default value
References
“DYNAMIC” TYPE_RESU
“DYNAMIC”
[§2.1]
“MODE_FLAMB”
[§2.1]
METHOD “SORENSEN”
“SORENSEN”
[§6.4]
CALC_FREQ
FREQ
[§4.4]
CHAR_CRIT
[§4. 4]
OPTION “PLUS_PETITE”
`PLUS_PETITE ''
[§4.4]
“TAPE”
[§4.4]
“CENTER”
[§4.4]
NMAX_FREQ
10
[§4.4]
DIM_SOUS_ESPACE
Calculated
[§6.5.3]
PREC_SOREN
0.
[§6.4], [§6.5.3]
NMAX_ITER_SOREN
20
[§6.4], [§6.5.3]
PARA_ORTHO_SOREN
0.717
[§6.5.3], [Appendix 2]
NPREC_SOLVEUR
8
[§2.6]
NMAX_ITER_SHIFT
5
[§2.6]
PREC_SHIFT
0.05
[§2.6]
SEUIL_FREQ
1.E-02
[§2.9]
VERI_MODE
STOP_ERREUR “YES”
“YES”
[§2.9]
“NOT”
[§2.9]
PREC_SHIFT
5.E-03
[§2.9]
THRESHOLD
1.E-06
[§2.9]
STURM “YES”
“YES”
[§2.9]
“NOT”
[§2.9]
Summary table 6.5.5-a: of the parameter setting of
MODE_ITER_SIMULT
with `
SORENSEN'
Note:
·
in V5, the user can specify the class of membership of his calculation by initializing it
key word
TYPE_RESU
. According to this value, one informs the vector
FREQ
or
CHAR_CRIT
,
·
one finds all the “tripaille” of parameters related to the preprocessings of the test of Sturm
(
NPREC_SOLVEUR
,
NMAX_ITER_SHIFT
,
PREC_SHIFT
) and with postprocessings of checking
(
SEUIL_FREQ
,
VERI_MODE
),
·
at the time of the first passages, it is strongly advised to modify only them
main parameters noted in fat. The others relate to more the mysteries of
the algorithm and they were initialized empirically with values standards,
·
in particular, to improve quality of a mode, the fundamental parameter is
dimension of the subspace
DIM_SOUS_ESPACE
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
60/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
7
Method of Bathe and Wilson (
METHOD = `JACOBI
)
7.1 Principle
The method of Bathe and Wilson is a method of simultaneous iterations which consists in extending
the algorithm of the iterations opposite. One works starting from the shifted problem
(
)
With X
B X
=
-
.
One distinguishes four parts in the algorithm [bib1], [bib4]:
·
to choose
p
independent initial vectors
X X
X
1
2
,
,…,
p
and to build the matrix
X
that they
generate,
·
to calculate the clean elements in the subspace of Ritz while solving
(
)
With
B U
WITH Q A Q
B Q BQ
-
=
=
=
I
I
T
T
0 where
and
. One returns then to space
initial (for the clean vectors) via the transformation
{}
[]
X QU
U
U
=
=
where
I
,
·
to test the convergence of the clean modes
I
.
7.2
Tests of convergence
The method of Bathe and Wilson converges towards
p
smaller eigenvalues provided that them
p
initial vectors are not
B
- orthogonal with the one of the clean vectors. In addition, matrices
With
and
B
tend towards diagonal matrices. For this reason and as the matrices are full,
one uses the method of Jacobi (cf [Appendix 3]) to find the elements clean of the subspace of
Ritz.
To test the convergence of the eigenvalues, one classifies them after each iteration by command
growing in absolute value and one looks at if, for each eigenvalue, the following test is checked
K
K
K
PREC BATHE
+
+
-
1
1
_
where the exponent
K
indicate the iteration count. If afterwards
NMAX_ITER_BATHE
iterations, one does not have
converged for all the eigenvalues, a message of alarm is transmitted in the file message.
7.3
Establishment in Code_Aster
7.3.1 Dimension of the subspace
If one wishes to calculate
p
eigenvalues, it is recommended to use under space of dimension
Q
higher. One will check convergence only for
R
smaller eigenvalues where
p R Q
. It
seem that
r= p
is not sufficient: one can find the good eigenvalues but the vectors
clean are not correct (convergence is slower for the clean vectors than for
eigenvalues).
(
)
R
p Q
=
+
/2
seem a good choice. For
Q
one usually takes [bib1]
(
)
Q
p
p
=
+
min
,
8 2
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
61/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
7.3.2 Choice of the initial vectors
To choose them
Q
initial vectors, one operates in the following way:
·
first vector such as
X
I
II
II
1
=
With
B
,
·
for the other vectors of 2 with
()
Q
-
1
X
X
2
1
1
2
0
1
0
0
1
0
0
=
- -
=
- -
-
line
line
I
I
Q
,…
where
I
are the indices the corresponding to smallest successive values of
With
B
II
II
,
·
the last
vector
X
Q
, random vector.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
62/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
7.3.3 Parameters in Code_Aster
To be able to use the method of Bathe and Wilson, the control should be chosen
MODE_ITER_SIMULT
and to select
METHOD = “JACOBI”
. Two parameters relating to convergence directly
method are accessible under the key word factor
CALC_FREQ
using the key words
PREC_BATHE
and
NMAX_ITER_BATHE
. One finds there also those managing the internal method of resolution
modal,
PREC_JACOBI
and
NMAX_ITER_JACOBI
.
Key word factor
Key word
Default value
References
“DYNAMIC” TYPE_RESU
“DYNAMIC”
[§2.1]
“MODE_FLAMB”
[§2.1]
METHOD “JACOBI”
“SORENSEN”
[Appendix 3]
CALC_FREQ
FREQ
[§4.4]
CHAR_CRIT
[§4. 4]
OPTION “PLUS_PETITE”
`PLUS_PETITE ''
[§4.4]
“TAPE”
[§4.4]
“CENTER”
[§4.4]
NMAX_FREQ
10
[§4.4]
DIM_SOUS_ESPACE
Calculated
[§7.3.1]
PREC_BATHE
1.E-10
[§7.2]
NMAX_ITER_BATHE
40
[§7.2]
PREC_JACOBI
1.E-12
[Appendix 3]
NMAX_ITER_JACOBI
12
[Appendix 3]
NPREC_SOLVEUR
8
[§2.6]
NMAX_ITER_SHIFT
5
[§2.6]
PREC_SHIFT
0.05
[§2.6]
SEUIL_FREQ
1.E-02
[§2.9]
VERI_MODE
STOP_ERREUR “YES”
“YES”
[§2.9]
“NOT”
[§2.9]
PREC_SHIFT
5.E-03
[§2.9]
THRESHOLD
1.E-06
[§2.9]
STURM “YES”
“YES”
[§2.9]
“NOT”
[§2.9]
Summary table 7.3.3-a: of the parameter setting of
MODE_ITER_SIMULT
with `
JACOBI'
Note:
·
in V5, the user can specify the class of membership of his calculation by initializing it
key word
TYPE_RESU
. According to this value, one informs the vector
FREQ
or
CHAR_CRIT
,
·
one finds all the “tripaille” of parameters related to the preprocessings of the test of Sturm
(
NPREC_SOLVEUR
,
NMAX_ITER_SHIFT
,
PREC_SHIFT
) and with postprocessings of checking
(
SEUIL_FREQ
,
VERI_MODE
),
·
at the time of the first passages, it is strongly advised to modify only them
main parameters noted in fat. The others relate to more the mysteries of
the algorithm and they were initialized empirically with values standards.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
63/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
8
Conclusion - Synthesis
The optimal perimeters of use of the modal operators of Code_Aster can be dissociated.
When it is a question of determining some eigenvalues (typically a half-dozen) or
to refine some estimates, the operator
MODE_ITER_INV
is completely indicated. It gathers
heuristic algorithms and those of type powers (cf [§3]), which were historically developed
first to solve generic modal problems. It couples a phase of localization of
eigenvalues (via a technique of bisection and a method of the secant) with an improvement
of these estimates and a calculation of the clean vectors associated by a method of iterations opposite
(crossbred, or not, of quotient of Rayleigh).
On the other hand, to capture a significant part of the spectrum, one A resorts to
MODE_ITER_SIMULT
.
This last federates the methods known as of “subspace” (Lanczos [§4], [§5], IRAM [§6], Bathe & Wilson
[§7]) which projects the operator of work in order to obtain an approximated spectrum of more reduced size (treated
then by a total method of type
QR
or Jacobi, cf [Appendix 1] and [Appendix 3]). In addition to theirs
numerical qualities (reduced calculation complexities and memory, facilitated coupling of techniques of
deflation, of restarting and acceleration…) and mathematics (good convergence, controls
quality of the spectrum of the projected operator…), it should be noted that they necessarily do not require
knowledge of the operator of work but that of its action on a vector (this characteristic is
very useful to treat large systems but it is not used in Code_Aster where one assembles
all matrices before treating them).
Until now, these algorithms stumbled regularly on the same shelves: detection
correct of multiple modes, modes of rigid body and a general way, processing of
packed spectrum. All this led to the appearance of “phantom” modes sometimes badly easy to detect
(modes corresponding to multiplicities missed and being able to generate correct residues with the direction
of Aster, and this, more especially as the criteria of the post-modal checks were sometimes permissive
(residue in 10
- 2
instead of the 10
- 6
current), decontaminated even insufficient (test of Sturm limited to the values
clean positive)) who cause distortions of results downstream from calculation, during projections on
base modal. To be been free from the recurring problems to this type of approach, one thus proposed
to enrich
MODE_ITER_SIMULT
(starting from V5) from the algorithm (“Implicit Restarted Arnoldi WILL GO”
[§6]).
This alternative of Arnoldi, initiated by D.C. Sorensen in 1992, makes real great strides for the resolution of
great modal systems on parallel supercomputers. Its framework of application is completely
General, it deals with as well the real problems as complex, square or not. He summarizes himself in one
succession of factorizations of Arnoldi whose results control restartings automatically
statics and implicit, via polynomial filters modelized by iterations
QR
implicit.
The algorithm tent will bring an elegant remedy for the recurring problems raised by the others
approaches. There is not any more a question to be posed concerning the strategy of reorthogonalisation,
frequency of the restarts, their establishments, the criteria of detection of possible phantom modes…
In short, numerically, IRAM gets a better total robustness while improving them
calculation complexities and memory (especially compared to simple Lanczos such as that of Newman &
Pipano) for a fixed precision.
From a functional point of view, the establishment of this method allowed a gain in calculation complexity
(observed) with minima of command 2. From now on, the user has a real control on the quality of the modes via one
suitable parameter setting. One, moreover, did not reinforce the severity of the post-modal checks and extended them
applicability.
The use of IRAM (per defect in
MODE_ITER_SIMULT
) is thus with advising in all them
case of figures, including for the search for some eigenvalues (one can even fire part
excellent properties of the method of the powers opposite concerning the calculation of vectors
clean and, to refine the result, while starting again
MODE_ITER_INV
with for estimate values
clean exhumed by
MODE_ITER_SIMULT
). Beyond their numerical specificities and
functional calculuses which are included in this document, one can synthesize the modal methods of
Code_Aster in the shape of table below (the default values are materialized in
fat).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
64/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Operator/
Perimeter
of application
Algorithm
Key word
Advantages
Disadvantages
MODE_ITER_INV
1
era
phase
(heuristics)
Calculation of some
modes
Bisection
“SEPARATE”
Calculation of some
modes
Bisection +
Secant
“ADJUSTS”
Better precision
Cost calculation
Improvement of
some estimates
Initialization by
the user
“NEAR”
Resumption of values
clean estimated
by another
process.
Cost calculation of this
phase quasi-no one
No the capture of
multiplicity
2
ième
phase
(method of
powers properly
said)
Basic method
Powers
opposite
“DIRECT”
Very good
construction of
clean vectors
Not very robust
Option of acceleration
Quotient of
Rayleigh
“RAYLEIGH”
Improve
convergence
Cost calculation
MODE_ITER_SIMULT
Calculation of part of
spectrum
Bathe & Wilson
“JACOBI”
Not very robust
Lanczos
(Newman- Pipano)
“TRI_DIAG”
Not very robust
IRAM (Sorensen)
“SORENSEN”
Increased robustness.
Better
calculation complexities and
memory.
Control
quality of the modes.
Summary table 8-a: of the modal methods of Code_Aster
IRAM made it possible to balance all the software faults related to the generalized modal problems,
but it seems that alone its variation per blocks or a version incorporating of the “purging and
lock " [bib23] can guarantee to us “an quasi-infallible” detection of the spectrum of a standard operator
(i.e too badly not conditioned).
A class of algorithm known as of “Jacobi-Davidson” [bib37] seems even more promising to treat
pathogenic cases. It uses an algorithm of the Lanczos type that it packaged via a method
of Rayleigh.
Let us conclude by noting that the algorithm WILL GO was not extended yet to the quadratic case in
modal operators of the code [R5.01.02].
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
65/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
9 Bibliography
[1]
K.J. BATHE: Broad solution methods for generalized eigenvalue problems in
structural engineering. California university Berkeley, 1971.
[2]
F. CHATELIN: Eigenvalues of matrices. Masson, 1988.
[3]
J.K. CULLUM & R.A. WILLOUGHBY: Broad Lanczos algorithms for symetric
eigenvalue computations. Flight 1 Theory, Birkhäuser 1985.
[4]
DHATT & TOUZOT: A presentation of the finite element method. Maloine
S.A. Edition, 1984.
[5]
A. DUBRULLE, R.S. MARTIN & J.H. WILKINSON: The Implicit
QL
Algorithm
pp241-248, in Handbook for automatic computation. Flight 2, Linear will algebra,
Springer-Verlag, 1971.
[6]
G.H. GOLUB & c.f. Van LOAN: Matrix computations, The Johns Hopkins
university near, 1989.
[7]
Y. HAUGAZEAU: Application of the theorem of Sylvester to the localization of
eigenvalues of
Ax
=
Bx
in the symmetrical case. RAIRO Analyzes numerical,
vol.14, n°1 pp25-41, 1980.
[8]
D. HO: Broad Tchebychev technical acceleration for scale nonsymmetric matrices.
Numerische mathematik. Flight 56, pp731-734.
[9]
J.F. IMBERT: Analyze structures by finite elements. Sup Aéro, CEPADUES
Editions.
[10]
C. LANCZOS: Year iteration method for the solution off the eigenvalue problem off
linear differential and integral operators. Newspaper. off research. off the national office
off standard, flight 45, n°4 pp255-282, Oct. 1950.
[11]
P. LASCAUX & R. THEODOR: Matric numerical analysis applied to the Art of
the engineer. Masson, 1986.
[12]
R.S. MARTIN, G. PETERS & J.H. WILKINSON: The
QR
Algorithm for Real
Hessenberg Matrices pp358-371, in Handbook for automatic computation. Flight.
2, Linear will algebra, Springer-Verlag, 1971.
[13]
R.S. MARTIN & J.H. WILKINSON: Similarity Reduction off has General Matrix to
Hessenberg From pp339-358, in Handbook for automatic computation. Vol. 2
Linear will algebra, Springer-Verlag, 1971.
[14]
Mr. NEWMANN & A. PIPANO: Modal Fast extraction in Nastran via FEER
computer programs, Nastran, to use manual, NASA Langley Research Center
pp485-506, 1977.
[15]
B. NOUR-OMID, B.N. PARLETT & R.L. TAYLOR: Lanczos versus subspace
iteration for solution off eigenvalue problems, Int. J. for numerical methods in
engineering. Vol. 19, pp859-871, 1983.
[16]
INTER-F OJALVO & Mr. NEWMANN: Vibrations modes off broad structures by year
automatic matrix-reduction method. AIAA, vol.8 n°7, 1970, pp1234-1239.
[17]
E.E. OSBORNE: One pre-conditioning off matrices. J. Assoc. Comput. Mach., flight 7,
pp338-345, 1960.
[18]
B.N. PARLETT: The symetric eigenvalue problem. Prentice hall, Englewoods
Cliffs, 1980.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
66/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
[19]
B.N. PARLETT & C. REINSCH: Balancing has matrix for calculation off eigenvalues
and eigenvectors. pp315-326, in Handbook for automatic computation. Vol. 2,
Linear will vlgebra, Springer-Verlag, 1971.
[20]
Y. SAAD: One the spleens off convergences off the Lanczos and the Block-Lanczos
methods. SIAM J. Numer. Anal., flight 17 N° 5, pp687-706, october 1980.
[21]
H.D. SIMON: The Lanczos algorithm with partial reorthogonalisation, Mathematics
off computation. Flight 42, n° 165, pp115-142, 1984.
[22]
Mr. SADKANE: With block Arnoldi-Tchebychev method for computing the leading
ergenpairs off broad sparse unsymmetric matrices. Numerische mathematik, flight 64,
pp181-193, 1993.
[23]
J.L. VAUDESCAL: Introduction to the methods of resolution of problems with
eigenvalues of big size. Note HI-72/00/01A, 2000.
[24]
J.H. WILKINSON & C. REINSCH: Handbook for automatic computation, flight 2,
Linear Algebra, Springer-Verlag, New York, 1971.
[25]
O. BOITEAU & B. QUINNEZ: Method of spectral analysis in Code_Aster.
Booklet of the Course of dynamics, March 2000.
[26]
D.C. PAIGE: Accuracy and effectivness off the Lanczos algorithm for the symmetric
eigenproblem. Linear will algebra and its applications, flight 34, 1980.
[27]
J.J. DONGARRA & Al: LAPACK' S users guide. SIAM, 1992.
[28]
J.J. DONGARRA & Al: EISPACK' S users guide. Springer verlag, flight 6, 1976.
[29]
R.B. LEHOUCQ, D.C. SORENSEN & C. YANG: ARPACK' S users guide. 1996.
[30]
D.C. SORENSEN: Implicit applications off polynomial filters in A k-step Arnoldi
method. SIAM J. Anal Matrix. Appl., flight 13, pp357-385, 1992.
[31]
W.E. ARNOLDI: The principle off minimized iterations in the solution off the matrix
eigenvalue problem. Quarter. Appl. Maths, flight 9, pp17-19, 1951.
[32]
Y. SAAD: Variation one Arnoldi' S method for computing eigenelements off broad
unsymetric matrices. Linear will algebra and its applications, flight 34, pp269-2395, 1980.
[33]
Y. SAAD: Broad Numerical meyhods for eigenvalue problems. Manchester
university, close series in algorithms and architectures for advanced scientific
computing, 1991.
[34]
Z. BAI, J. DEMMEL & A.M. KENNEY: One computing condition numbers for the
nonsymmetric eigenproblem. ACM transactions one mathematical software, flight 19,
pp202-223, 1993.
[35]
T. ERICSSON & A. RUHE: Spectral The transformation Lanczos method for the
numerical solution off broad sparse generalized symmetric eigenvalue problems.
Mathematics off computations, flight 35, pp1251-1268, 1980.
[36]
J.J. DONGARRA & Al: Year extended set off FORTRAN BASIC linear will algebra
subprograms. ACM transactions one mathematical software, flight 14, pp1-17, 1998.
[37]
R.B. MORGAN & D.S. SCOTT: Preconditionning the Lanczos algorithm for sparse
symmetric eigenvalue problems. SIAM J. Sci. Comp, flight 14, 1993.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
67/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Appendix 1 General information on the algorithm
QR
A1.1 Principle
Algorithms of the type
QR
(cf [§2.8]) were had a presentiment of per H. Rutishauser (1958) and were formalized
jointly by J.C. Francis and V.N. Kublanovskaya (1961). This fundamental method is
often implied in the other approaches adapted better to deal with the problems of large
sizes (in particular methods of projection).
For K 1,… to make
=
(factorization
),
=
;
Fine loops.
K
K
K
K +1
K
K
=
H
Q R
QR
H
R Q
Algorithm 1.1:
QR
theoretical
The process leads repeatedly towards a matrix
H
K
triangular higher (or triangular
by blocks) whose diagonal terms are the eigenvalues of the initial operator
H
=
H
1
. The notation
H
is not in fact not innocent, because there is any interest with transforming beforehand orthogonally (
manner to modify only the shape of the shifté operator and not his spectrum) the operator of work
With
in the form of higher Hessenberg, that is to say into arithmetic real
H
Q AQ
1
0
0
=
T
This can be carried out via various orthogonal transformations (Householder, Givens,
Gram-Schmidt…) and their cost (about
O (10n
3
/3)
) negligible is compared with the gain that they
allow to realize with each iteration of the total process:
O (N
2
)
(with Householder or
Fast-Givens one has more precisely
O (2n
2
)
against
O (4n
2
)
with simple Givens)
against
O (N
3
)
. It
gain of a command magnitude can be even improved when the operator is tridiagonal symmetrical
(it is the case of Lanczos with a true scalar product):
O (20n)
.
Convergence towards a simple triangular matrix is carried out only if all the eigenvalues
are distinct modules and that if the initial matrix is not “pathologically” too poor in
clean directions. The convergence of
I
ème
mode (arranged classically by order descending of
modulate) is carried out then in:
max
J I
J
I
,
what can prove very slow if no complementary process is implemented.
Note:
·
determination of the spectrum of the matrix of Rayleigh with the alternative of Newmann & Pipano
(cf [§5.5.1]) is carried out via one
QR
(or one
QL
into symmetrical) simple of this type (with, with
precondition, balancing). The only parameter accessible by the user is the maximum number
acceptable iterations
NMAX_ITER_QR
,
·
for IRAM (cf [§6.4]), this calculation is carried out via a method
QR
with double explicit shift
whereas the polynomial filters managing the restarts use one
QR
with double implicit shift.
By prudence, no parameter is accessible for the user in Code_Aster!
·
the method, the algorithm should not be confused
QR,
and one of its conceptual tools,
factorization
QR
,
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
68/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
·
this class of algorithm is very much used to determine the complete spectrum of an operator,
because it is very robust (it is the reference in this field). However it is very
greedy memory from there what places makes its use crippling on great systems,
·
the perimeter of application of the algorithm
QR
is much more general than that of Jacobi
(it is the second standard algorithm providing all the spectrum of an operator for can
that one is ready to entirely store it) who is limited to the square matrices.
To accelerate the convergence of the simple algorithm which can be very slow (in the presence of clusters
for example) a multitude of alternatives, based on the choice of shifts answering certain criteria,
were born.
A1.2 strategy of the shift
This strategy consists in causing artificially a phenomenon of deflation (cf [§3.1], [§5.4.1])
within the matrix of work. That offers triple favors:
·
of
to be able to isolate a real eigenvalue even two combined complex eigenvalues,
·
all
in
reducing the size of the problem to be treated,
·
and
in
accelerating convergence.
In its version with simple explicit shift, the method is rewritten then in the following form:
For K 1,… to make
To choose the shift,
=
(factorization
),
=
+
;
Fine loops.
K
K
K
K
K +1
K
K
=
=
-
µ
µ
µ
S
H
I
S
Q R
QR
H
R Q
I
K
,
Algorithm 1.2:
QR
with simple explicit shift
Note:
This process spreads intuitively with several shifts. One builds then, for each
total iteration
K
, as many auxiliary matrices
S
K
I
that of shift
µ
I
.
The convergence of the process is largely improved in the direction where under-diagonal terms
cancel themselves asymptotically in:
max
J I
J
I
-
-
µ
µ
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
69/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
In theory, if this shift
µ
is eigenvalue of the problem, then deflation is exact. In practice, them
effects of flare disturb this phenomenon, it is what is called the property of direct instability of
the algorithm. The main difficulty lies in the choice of (or of) the shifts. In addition, one does not keep
not the same shift for all the iterations. One must change some when it is associated a value
clean converged. Indeed, it will have numerically caused its ousting of the spectrum of work in
causing a deflation with the preceding iteration.
Since the Sixties a whole zoology of shifts developed. While simplest
use the last diagonal term
H
N, N
, that of J.H. Wilkinson [bib24] consists in determining
analytically the eigenvalue
µ
diagonal block
H
H
H
H
N
N
N
N
N N
N N
-
-
-
-
1
1
1
1
,
,
,
,
nearest to this term. This technique makes it possible to even obtain a quadratic convergence
cubic (in the symmetrical case) and it proves particularly effective to capture modes
double or of the distinct eigenvalues of the same module. However this strategy can appear
ineffective in the presence of combined complex clean modes.
The same author then proposed an alternative including a double shift corresponding to both
complex eigenvalues
µ
1
and
µ
2
(beforehand given) of the accused block. But in addition to
numerical difficulties with réfreiner appearance of invading complex components (which in theory
do not take place to be), one has in addition much evil to preserve the character of “Hessenberg
higher " of the matrix of work. As well as possible, that slowed down the convergence of the algorithm
QR
, with
worse, that distorts its operation completely.
In order to bearing with these underhand numerical disadvantages, J.C. Francis developed a version with
double implicit shift. It is very well clarified in [bib6] (pp377-81) and one will be satisfied just here
to summarize philosophy of it.
To minimize the effects of flares, would have to be constituted directly the auxiliary matrix
S
K
resulting
simultaneous application of the two shifts
µ
1
and
µ
2
(
)
(
)
S
H
H
I
K
K
K
=
-
+
+
2
1
2
1 2
µ
µ
µ µ
before factorizing it in form
QR
and to build the new one reiterated
H
k+1
S
Q R
H
Q H Q
K
K
K
K
K
T
K
K
=
=
+
,
.
1
But only the cost of the initial assembly (about
O (N
3
)
) the tactics make inoperative. This
alternative, based on the theorem
Q
- implicit, consists in applying to the matrix of work of
particular transformations of Householder allowing to find a matrix of Hessenberg
“primarily” equalizes with
H
k+1
, i.e. of the type:
(
)
~
H
E H
E
E
K
K
+
-
+
=
±
±
1
1
1
with = diag 1,… 1
.
The matrix of work thus preserves, at lower cost, its particular structure and its spectrum.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
70/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
All these alternatives are in fact very sensitive to the techniques of balancing established for
préconditionner the initial operator. The following paragraph will summarize this technique of “setting with
the scale " (“balancing” for the Anglo-Saxons) of the terms of the matrix of work which is very
employee in modal calculation.
A1.3 Balancing
It is a question of mitigating the effects of round by attaching the perimeter of expansion of the terms of
the operator of work, i.e. to prevent that they do not become too small or, on the contrary, too large.
On this subject, E.E. Osborne [bib17] (1960) noticed that generally the error on the calculation of
clean elements of
With
is about
With
2
where
is the precision machine. It proposed then
to transform the initial matrix into a matrix
~
WITH D AD
=
-
1
(with
D
a diagonal matrix),
such as:
~
With
With
2
2
<<
.
In fact one calculates a succession of matrices repeatedly
With
D
With
D
K
K
K
K
=
-
-
-
-
1
1
1
1
such as:
With
With
K
K
F
checking
has
has
I
I
2
2
=
with
has
I
and
has
I
, respectively, them
I
ème
column and line of
With
F
.
Note:
·
this technique was generalized to any matric standard induced by the discrete standard of
L
Q
with
L
Q
the whole of the complex continuations
(
U
N
)
N
such as
U
N
N
Q
<
,
·
its employment is very widespread in scientific computation and in particular among the direct solveurs of
system of linear equations.
The base of calculation of the computer, noted
, intervenes in the determination of the terms of the matrices
D
K
.
In order to minimize the round-off errors, one chooses the elements of
D
K
so that they are powers
of this base.
Are
R
K
and
C
K
p
standards (in practice one often takes
p= 2
), respectively of the line and
column
I
matrix
With
k-1
(
I
is selected in such way that
[]
I
K
N
- -
1
1
). By supposing that
R C
K
K
0
, it is shown whereas there is a single signed entirety
such as:
2
1
2
1
-
+
<
R
C
K
K
.
That is to say
F
=
, the matrix is defined
D
K
such as:
(
) (
)
(
)
D
Id
E E
C
R
C
R
Id
K
I iT
K
p
K
p
kp
kp
F
F
F
=
+
-
+
<
+
(
)
/
1
if
if not
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
71/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
where
0
1
<
is a constant and
E
I
I
ème
canonical vector. One builds then repeatedly:
D
D D
With
D WITH
D
K
K
K
K
K
K
K
=
=
-
-
-
1
1
1
,
.
,
while initializing
D
0
=
Id
. The process stops as soon as
D
Id
K
=
.
Note:
·
before carrying out balancing, a search of the isolated eigenvalues is carried out in
detecting the presence of lines and quasi-null columns (all the terms are null, except
that placed on the diagonal). When there are some one can, by carrying out permutations
suitable, to put the matrix of work in the following form blocks:
*
*
*
*
*
*
*
*
0
0
0
0
0
Y
Z
X
T
It is then necessary to carry out balancing only on the central block
X
because terms
diagonal of the two higher triangular matrices are eigenvalues of the matrix
of work.
·
in Code_Aster, one chose
p
=
1
and
=
0 95
.
(cf [bib19]),
·
before applying the method
QR
, one starts by balancing the matrix of work and then,
one transforms it in the form of Hessenberg higher. Once calculation
QR
carried out, one
go up vectors of Ritz to the clean vectors via the reverse of the transformations
orthogonal due to the setting in Hessenberg form and balancing. It is this process which
is set up as well in Lanczos in IRAM. However in addition to the fact that it
require the storage of these orthogonal matrices, it is also and especially extremely
sensitive to the effects of flares. Thus, it would be largely preferable to estimate the vectors
clean via a method of the powers opposite initialized by the eigenvalues
exhumed.
A1.4 method QL
Factorization
QL
of a matrix
With
consist with orthonormaliser its vectors columns of
starting with the last (contrary to the algorithm
QR
who begin with the first) giving thus
a matrix
L
triangular lower regular. It is shown besides that the algorithm
QL
applied to
the operator
With
invertible is equivalent to the algorithm
QR
applied to
With
-
*
(transposed matrix
combined reverse of
With
).
The method installation in Code_Aster is identical to the method of simple shift of
J.H. Wilkinson seen previously by adapting the search for this shift to the miner 2x2 highest.
Note:
Contrary to
QR
who captures the eigenvalues by order ascending of module, one obtains here
preferentially dominant modes, then others, by order descending of module.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
72/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Appendix 2 Orthogonalization of Gram-Schmidt
A2.1 Introduction
One saw on several occasions in this document that the quality of orthogonalization of a family of
vectors is crucial for the good unfolding of the algorithms and the quality of the modes obtained.
This task is besides to permanently realize from where the importance of a fast algorithm.
The simple algorithms of orthonormalisation are deduced from the conventional process of Gram-Schmidt
but they are often “conditionally stable” (the quality of their work depends on
matric conditioning of the family with orthonormaliser). This defect of robustness can prove
problems to treat situations particularly badly conditioned. One prefers to them then
more expensive but robust algorithms, containing projections or of rotations: transformations
spectral of Householder and Givens.
In practice, for the establishment of the method WILL GO in Code_Aster, we retained one
iterative version of the process of Gram-Schmidt (the IGSM of Kahan-Parlett [bib18]) It carries out one
good compromise between the robustness and calculation complexity since it is unconditionally
stable and orthogonaliser allows except for the precision machine, in, to the maximum, twice more
time that conventional Gram-Schmidt (GS).
In the following paragraphs we will detail the operation of the basic algorithm, thus
that that of these two main alternatives. First was installation in
MODE_ITER_INV
and
second is used in
MODE_ITER_SIMULT
. Comparative very percussion of these methods is
declined in [bib23] (pp33-36) starting from a very simple example.
A2.2 Algorithm of Gram-Schmidt (GS)
Being given
K
vectors independent of
N
,
()
X
I I
K
=
1,
one wishes to obtain
K
orthonormal vectors
(compared to an unspecified scalar product)
()
y
I I
K
=
1,
space which they generate. In others
terms, one wishes to obtain another orthonormal family generating same space. The process
of conventional orthonormalisation of Gram-Schmidt is as follows:
()
For I 1, K to make
For J 1, I - 1 to make
Calculation of R
Fine loops.
Calculation of
Calculation of
Calculation of
Fine loops.
ji
=
=
=
=
-
=
=
= -
y X
y
X
y
y
y
y
J
I
I
I
ji
J
I
J
II
I
I
I
II
R
R
R
,
;
%
,
%,
%;
,
1
1
Algorithm 2.1: Algorithm of Gram-Schmidt (GS)
This process is simple but very unstable because of the round-off errors, which causes to produce
nonorthogonal vectors. In particular when the initial vectors are almost dependant
that creates important variations magnitude in the second stage of the process
%
,
y
X
y
I
I
ij
J
I
J
R
=
-
= -
1
1
From a numerical point of view, the management of these variations is very difficult.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
73/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Note:
·
the problem is completely similar to that met by the inversions of systems put in
place in the code (cf [§2.6]),
·
while noting
Q
the matrix generated by
()
y
I I
K
=
1,
, one explicitly was thus built
factorization
QR
initial matrix
X
bound to
()
X
I I
K
=
1,
. It is in fact the goal of very proceeded
of orthogonalization.
A2.3 Algorithm of Modified Gram-Schmidt (GSM)
In order to évincer these instabilities numerical one reorganizes the preceding algorithm. Mathematically
equivalent with the preceding process, the aforementioned is then much more robust because it avoids the variations of
magnitude important between the vectors handled in the algorithm.
In the initial process, the orthogonal vectors
y
I
are obtained without taking account of
i-1
preceding orthogonalizations. With Modified Gram-Schmidt, one orthogonalise more
gradually by taking account of preceding deteriorations according to the process below.
()
For I 1, K to make
For J 1, I - 1 to make
Calculation of R
Calculation of
Fine loops.
Calculation of
Calculation of
Fine loops.
ji
=
= =
=
=
-
=
=
%
,
, %,
%
%
;
%,
%;
y
X
y y
y
y
y
y
y
y
I
I
J
I
I
I
ji
J
II
I
I
I
II
R
R
R
Algorithm 2.2: Algorithm of Modified Gram-Schmidt (GSM)
The orthonormality of the basic vectors is much better with this process and it can be even obtained with
the precision machine and with a constant (depend on conditioning on
X
) near. However, for
to treat situations particularly badly conditioned, this “conditional” stability can prove
quickly problematic, from where the recourse to the following iterative algorithm.
Note:
GSM is twice more effective than a method of Householder to obtain a factorization
QR
initial matrix
X
. It only requires
O (2nk
2
)
operations (with
N
the number of lines
matrix).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
74/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
A2.4 Algorithm of Iterative Gram-Schmidt (IGSM)
To ensure itself nevertheless of orthogonality except for the precision machine, one recommends to realize
one second orthogonalization. And if, at the end of the latter, orthogonality is not assured it
is not any more the sorrow to start again, the handled quantities are then certainly very close and
their variations oscillate around zero. This approach is based on a theoretical analysis due to
W. Kahan and taken again by B.N. Parlett [bib18] (cf pp105-110).
()
()
For I 1, K to make
For J 1, I - 1 to make
Calculation of R
Calculation of
If
then
Exit loops in J
If not
Calculation of R
Calculation of
If
then
Exit loops in J
If not
,
Passage to I following;
End if.
Fine loops.
Calculation of
Calculation of
ji
ji
I
=
= =
=
=
-
=
=
-
=
%
,
, %,
~
%
,
~
%
%
~,
;
~
, ~,
~
~
~
~
,
~
~
~
%
~
~,
;
%
%,
y
X
y y
y
y
y
y
y
y
y
y y
y
y
y
y
y
y
y
y
0
y
y
I
I
J
I
I
I
ji
J
I
I
I
I
J
I
I
I
ji
J
I
I
I
I
II
I
R
R
R
I
I
II
R
=
%;
y
Fine loops.
Algorithm 2.3: Algorithm of Iterative Gram-Schmidt of type Kahan-Parlett (IGSM)
During the use of IRAM in Code_Aster, the criterion
is parameterized by the key word
PARA_ORTHO_SOREN
(cf [§6.5]). It is shown that its interval of validity is
[1.2e, 0.83-
]
with
precision machine and one generally allot value 0.717 to him (by defect).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
75/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
Note:
·
more the value of the parameter
is large, less the reorthogonalisation starts, but
that affects the quality of the process,
·
contrary to the version “house” developed for Lanczos (cf [§5.5.1]) here criteria
stops concentrate on the standards of the vectors orthogonalized rather than on the products
scalars inter-vectors which tend more to reflect the effects of flare. That, joined to
suppression of the iterations higher than two, therefore useless, can explain the addition
of effectiveness of the version of Kahan-Parlett,
·
according to D.C.Sorensen the paternity of this method is rather to allot to J. Daniel and Al
(paper of 1976 subjected to Mathematics off Computation, vol.30, pp772-795).
Appendix 3 Method of Jacobi
A3.1 Principle
Method of Jacobi [bib4], [bib6], [bib11] allows to calculate all the eigenvalues of one
generalized problem whose matrices are definite positive and symmetrical (the matrices obtained with
each iteration by the method of Bathe &Wilson check these properties; cf [§7]). It consists with
to transform the matrices
With
and
B
problem
With U
B U
=
in diagonal matrices, while using
successively orthogonal similar transformations (matrices of rotation of Givens).
process can be schematized in the following way:
With
With
B
B
With
Q WITH Q
B
Q B Q
With
Q
With
Q
B
Q
B
Q
1
1
2
1
1 1
2
1
1 1
1
1
1
1
1
1
=
=
=
=
=
=
-
-
-
-
-
-
T
T
K
K
T
K
K
K
K
T
K
K
…
With
With
B
B
K
K
D
K
K
D
stamp
diagonal
stamp
diagonal
Algorithm 3.1: Process of Jacobi
The eigenvalues are given by
()
=
-
WITH B
D
D
1
that is to say
=
With
B
II
D
II
D
and stamps it clean vectors
check:
X Q Q
Q
With
With
With
=
1 2
11
22
1
1
1
…
…
/
/
/
K
D
D
N
D
Each matrix
Q
K
is selected so that a term
()
I J
,
nondiagonal and not no one of
With
K
or of
B
K
that is to say no one after the transformation (for the choice of this term, cf [§0]).
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
76/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
A3.2 Some choices
In this algorithm, one realizes that the important points are as follows:
·
How to choose the terms to be cancelled?
·
How to measure the diagonal character of the matrices when
K
tends towards the infinite one?
·
How to measure convergence?
A3.2.1 nondiagonal Terms to cancel
For the choice of the terms to be cancelled, there are several methods:
·
the first consists with each stage
K
, to choose the largest element modulates some not
diagonal of the matrix
With
K
or
B
K
and to cancel it by carrying out a rotation (cf [§0]). This choice
ensure the convergence of the method of Jacobi but is relatively expensive (search for
the maximum element),
·
the second solution consists in successively cancelling all the nondiagonal elements of
these matrices while following the natural command
With
With
With
13
1
23
K
N
K
K
,…,
,
,…
. When one arrives at
With
N
N
K
-
1,
, one
start again the cycle. This method converges slowly.
An alternative of this method, consists very of following the natural command of the terms, to cancel
only those which are higher than a precision
K
data. At the end of a cycle, one decreases
value of this criterion and one starts again. It is this strategy which is used in Code_Aster.
A3.2.2 Test of convergence
To test the convergence and the diagonal character of the matrices, one operates thus. It is checked that all
coupling coefficients defined by:
F
F
I
J
ij
ij
ij
II
jj
ij
II
jj
With
B
With
WITH A
B
B B
=
=
are lower than a given precision (diagonal character of the matrices). One also controls
convergence of the eigenvalues via the indicator
F
Max
I
I
K
I
K
I
K
=
-
-
-
1
1
so that there remains lower than a given precision
jaco
.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
77/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
A3.2.3 Algorithm in Code_Aster
The algorithm implemented in Code_Aster is summarized with:
Initialization of the matrix of the vectors specific to the matrix identity.
Initialization of the eigenvalues.
To define the precision of necessary dynamic convergence.
To define the precision total
glob
.
For each cycle k=1,
N
jaco
_ max
To define the dynamic tolerance:
()
K
dyn
K
=
,
L
=
0
,
For each line i= 1, N
For each column j= 1, N
Calculation of the coupling coefficients
F
F
I
J
ij
K L
ij
K L
ij
K L
II
K L
jj
K L
ij
K L
II
K L
jj
K L
With
B
With
WITH A
B
B B
,
,
,
,
,
,
,
,
=
=
,
If
F
ij
K L
With
,
or
F
ij
K L
K
B
,
then
Calculation of the coefficients of the rotation of Givens,
Transformation of the matrices
With
K L
,
and
B
K L
,
,
Transformation of the clean vectors,
L
L
= +
1
End if.
Fine loops.
Fine loops.
Calculation of the eigenvalues
I
K
II
K L
II
K L
=
With
B
,
,
.
Calculation of
F
Max
I
ik
ik
ik
=
-
-
-
1
1
.
Calculation of the total coupling coefficients
F
Max
F
Max
I J
I J
ij
K L
II
K L
jj
K L
I J
I J
ij
K L
II
K L
jj
K L
With
B
With
WITH A
B
B B
=
=
,
,
,
,
,
,
,
,
If
F
F
F
glob
glob
glob
With
B
and
and
then
Correction of the clean modes (let us divide by
B
II
K
),
Exit;
End if.
Fine loops.
Algorithm 3.2: Method of Jacobi established in Code_Aster
One notes
N
jaco
_ max
the maximum number of allowed iterations.
Note:
Matrices
With
and
B
being symmetrical, only the half of the terms is stored in form
of a vector.
Code_Aster
®
Version
5.0
Titrate:
Algorithm of resolution for the generalized problem
Date:
02/03/01
Author (S):
O. BOITEAU
Key:
R5.01.01-C
Page:
78/78
Manual of Reference
R5.01 booklet: Modal analysis
HI-75/01/001/A
A3.2.4 Stamps rotation
One seeks with each stage to cancel the terms being in position
I
and
J
()
I
J
matrices
With
K L
,
and
B
K L
,
by multiplying them by a matrix of rotation which with the following form:
G
G
G
L
ij
ji
L
N
has
B
=
=
=
=
1
1,
;
;
other terms being null. One wishes to have
With
B
ij
K L
ij
K L
,
,
+
+
=
=
1
1
0
what leads to:
(
)
(
)
has
B has
B
has
B has
B
II
K L
ij
K L
jj
K L
II
K L
ijk L
jj
K L
With
With
With
B
B
B
,
,
,
,
,
,
+ +
+
=
+ +
+
=
1
0
1
0
because
With
G WITH G
B
G B G
K L
T
K, L
K L
T K L
,
,
,
+
+
=
=
1
1
and
. If the two equations are proportional one a:
has
and B
ijk L
jj
K L
=
= -
0
With
With
,
,
if not:
has
C
D
B
C
D
=
= -
2
1
with:
C
C
II
K L
ij
K L
II
K L
ij
K L
jj
K L
ij
K L
jj
K L
ijk L
1
2
=
-
=
-
With
B
B
With
With
B
B
With
,
,
,
,
,
,
,
,
()
C
D
C
sign C
C
C C
II
K L
jj
K L
II
K L
jj
K L
3
3
3
3 2
1 2
2
2
=
-
=
+
+
With
B
B
With
,
,
,
,
Note:
·
if
D
=
0
then one is in the case proportional,
·
if the matrix
B
is definite positive then
C
C C
3 2
1 2
2
+
is positive [bib11] what gives one well
feel with the parameter
D
.
Summary A3.2.5 of the parameter setting
To reach the parameters which intervene in the method of Jacobi, the user must select
following key words:
·
dyn
key word
PREC_JACOBI
under the key word factor
CALC_FREQ
,
·
N
jaco
_ max
key word
NMAX_ITER_JACOBI
under the key word factor
CALC_FREQ
.
Total precision
glob
is equal to the precision of convergence requested in the method of
Bathe and Wilson. One thus has
glob
bath
=
(cf [§7]).