Skip to main content

Nonlinear eigenvalue and frequency response problems in industrial practice

Abstract

Background

We discuss the numerical solution of large scale nonlinear eigenvalue problems and frequency response problems that arise in the analysis, simulation and optimization of acoustic fields. We report about the cooperation with the company SFE in Berlin. We present the challenges in the current industrial problems and the state-of-the-art of current methods.

Results

The difficulties that arise with current off-the-shelf methods are discussed and several industrial examples are presented.

Conclusions

It is documented that industrial cooperation is by no means a one-way street of transfer from academia to industry but the challenges arising in industrial practice also lead to new mathematical questions which actually change the mathematical theory and methods.

1 Background

1.1 Introduction

Traffic noise emissions by transport vehicles, such as cars, trains or airplanes are one of the key factors restricting the quality of life in urban areas. Acoustic waves in dynamically moving vehicles arise from many different sources, such as, for example, the noise of engines or the vibrations of the structure due to external excitations like road contact or head wind. The reduction of such noise emissions is therefore an important factor in the design and the economic success of new products.

To minimize potential noise emissions already in an early design phase of new products (avoiding the costly production of prototypes), it is necessary to use simulation, optimization and control techniques based on mathematical models of the vehicle. Performing these tasks requires mathematical models that describe the acoustic field associated with the structure of the complete vehicle including its interaction with the environment such as, for example, the air in and around the vehicle or the road or rail contact. Furthermore, these models must allow to identify and influence potential noise sources, and as a vision for the future they must allow to minimize the emissions.

Although much research is carried out in universities and industrial research and development departments, the model based minimization of the noise emissions of a complete car or train (not to mention airplane) including the majority of external and internal excitations is still a vision for the future. To achieve this, a joint effort is needed that includes the identification of sources, the construction of adequate mathematical models that incorporate all possible sources for acoustic waves, the analysis of these models, concerning robustness of their descriptions and their potential for optimization techniques, the development of numerical methods for the simulation and optimization of these models, as well as the implementation of these techniques as production software on modern high performance computers for industrial use.

1.1.1 Modeling

The modeling of acoustic fields inside or outside of dynamically moving vehicles typically uses coupled systems of partial differential equations (PDEs), for example, for the generation of noise by vibrating parts, surface contact, engine noise or head wind. These methods are well established in the engineering community using commercial finite element packages [1, 2]. However, the techniques for the resulting systems of PDEs still have many deficiencies, in particular the development of appropriate solvers for the solution of the linear and nonlinear systems and eigenvalue problems that have to be solved after discretization. These turn out to be extremely large and ill-conditioned - that is, extremely sensible with respect to perturbations in the problem-defining data - when a reasonably fine three-dimensional model is used; they may consist of hundreds of millions of equations.

An even greater challenge is to use these models and methods within an optimization loop. There is no chance to use classical off-the-shelf optimization methods for these problems, the problem size is just too large. Instead, one currently applies model reduction techniques [3, 4] to approximate the given fine model that is used for the simulation by a rather crude model that is used in the optimization.

To make such an approach viable within a design environment, where not only the geometry and topology of the vehicle and its material parameters are subject to changes, but also the interaction with the environment is rather complicated, it is necessary to develop integrated mathematical techniques for the modeling, simulation, and optimization that make as much use as possible of the properties of the underlying physical model and to transfer this into numerical methods and production codes.

1.1.2 Content of the paper

In this paper we will discuss some of the work and some of the challenges in a long-term cooperation with the company SFE GmbH in Berlin, Germany, which produces software for the simulation and optimization of acoustic fields.

The cooperation involves the development (and implementation on current high-performance computers) of numerical methods for the solution of linear systems and nonlinear eigenvalue problems arising from discretized partial differential equations modeling noise emissions of cars and trains.

The paper is organized as follows.

Section 1.2 briefly describes the modeling of acoustic fields inside a car as coupling of fluid and structure vibrations. In Section 2.1 we study the direct frequency response problem, that is, the numerical solution of a series of large, sparse, complex symmetric linear systems with a varying parameter. We will address in particular in-core/out-of-core storage methods, preconditioning and parallel execution aspects.

Section 2.2 treats modal reduction to reduce the dimension of the generated models. This requires the computation of all eigenvalues of a large sparse nonlinear eigenvalue problems in a given region of the complex plane.

The described problems and results demonstrate many challenges in the transfer of state-of-the-art numerical methods into the industrial practice, and show that there is mutual benefit from such a cooperation between industry and academia.

1.2 Mathematical modeling of interior car acoustics

To model the propagation of acoustic waves inside a car, as illustrated in Figure 1, the three-dimensional lossless wave equation is used, which can be obtained from the continuity equation (conservation of mass)

ρ ˜ t +( ρ ˜ v)=0,
(1)

together with the Euler equations of fluid dynamics

ρ ˜ ( v t +(v)v)= p ˜ .
(2)

Here v=v(x;y;z;t) is the particle velocity, ρ ˜ = ρ ˜ (x;y;z;t) is the particle density, and p ˜ = p ˜ (x;y;z;t) denotes the pressure, depending on Cartesian coordinates x y z and time t. Under the simplifying assumptions that there is no temperature change, that the fluid is inviscid, that is, no shear forces occur, that the influence of external forces is restricted to those coming from displacements of the structure at the boundaries, that the fluid is adiabatic, that is, there is no heat exchange during compression, that we have an ideal gas, that is, ρ= p c 2 , where c is the speed of sound, and finally that (v)v and ρ v t are small, we can make the Taylor expansions

p ˜ = p 0 + p ( x ; y ; z ; t ) , with  p 0 p ( p 0 = 10 6 p ) , ρ ˜ = ρ 0 + ρ ( x ; y ; z ; t ) , with  ρ 0 ρ .

Then from the Euler equation we obtain ρ 0 ( v t )=p or by differentiation ρ 0 ( v t )=Δp, and from the continuity equation we get ρ t + ρ 0 v=0. Altogether we obtain the highly simplified system of partial differential equations

1 c 2 2 p t 2 + ρ 0 v t =Δp+ ρ 0 v t =0

to which we have to add appropriate initial and boundary conditions. In a closed structure like the interior of the car, we can use the boundary conditions that are obtained from the displacement of the structure to obtain the fluid structure interaction, in an open structure like the acoustic waves emitted by the car to the outside, we have to incorporate appropriate far field conditions or transparent boundary conditions [5].

Figure 1
figure 1

FEM model of a car structure and illustration of the propagation of acoustic waves inside a car.

Let u be the vector of displacements of the structure on the surface. Then v= u t and thus with the outer normal ν we get

ν ρ 0 2 u t 2 =νp.

To obtain a variational formulation, we multiply the equations with a test function w, and use integration by parts by integrating over control volumes V with surface elements S. This gives

V 1 c 2 w 2 p t 2 dV+ V (w)pdV= ρ 0 S νw 2 u t 2 dS,

or equivalently

V 1 ρ 0 c 2 w 2 p t 2 dV+ V 1 ρ 0 (w)pdV= S νw 2 u t 2 dS.

One of the most difficult tasks in the modeling is the incorporation of appropriate damping and absorption, because this depends in a rather complicated way on the materials used inside the car, the surface geometry of the interior and this is also one of the terms, where the influence of the optimization is strong.

In the model treated here, damping and absorption was realized by an additional term

S w r ρ 0 2 c 2 p t dS,

where r=r(α) is depending on the damping properties of the material (described by a parameter vector α). We then obtain

V 1 ρ 0 c 2 w 2 t 2 pdV+ S w r ρ 0 2 c 2 p t dS+ V 1 ρ 0 (w)pdV= S νw 2 u t 2 dS.

To this variational formulation we can directly apply finite element discretization in space. This leads (in time) to a second order system of implicit differential equations in descriptor form

M f p ¨ d + D f p ˙ d + K f p d + D s f u ¨ d =0,
(3)

where M f = M f T is a positive definite mass matrix, K f = K f T is a positive definite stiffness matrix, D f (α) is a symmetric positive semi-definite damping/absorption matrix, and D s f describes the fluid structure coupling.

It should be noted that for the fluid model, the finite element discretization applied to the weak form of the partial differential equation is used. In contrast to this in the model describing the vibration of the structure, in industrial and engineering practice, typically a direct discrete finite element model is employed. This leads to one of the challenges that we will discuss below.

For the displacement vector u d in the different finite elements of the discretized structure, assuming linear material laws, one directly obtains a linear second order system of differential-algebraic equations given by

M s u ¨ d + D s u ˙ d + K s u d = f e + f p ,
(4)

where f e is a (discrete) external load and f p is the pressure load. Here M s is a symmetric positive definite mass matrix, and D s is a symmetric positive semi-definite damping matrix; both are real. The matrix K s has the form K s = K 1 (ω)+ı K 2 with real symmetric K 1 K 2 , where K 1 is the positive semi-definite stiffness matrix. It is often frequency dependent to model nonlinear material behavior. The matrix K 2 models hysteretic damping, that is, damping that is proportional to the displacement (instead of the velocity), but is in phase with velocity [6]. Typically, the matrix K 2 is of very small rank.

Although M s is positive definite in theory, the matrix encountered in practice is highly singular due to the fact that rotational masses are omitted. On the positive side, it is block diagonal with small blocks.

It remains to further discuss the coupling of the fluid part and the structure part of the system. The term f p originates from the pressure load F p = S pνdS. Finite element modeling/discretization yields f p = D s f T p d and hence,

M s u ¨ d + D s u ˙ d + K s u d D s f T p d = f e .
(5)

If we combine all the equations then we obtain a second order system of differential-algebraic equations

[ M s 0 D s f T M f ] [ u ¨ d p ¨ d ] + [ D s 0 0 D f ] [ u ˙ d p ˙ d ] + [ K s ( ω ) D s f T 0 K f ] [ u d p d ] = [ f s 0 ] .
(6)

Typically in this model the format of M s is a factor 1,000-10,000 larger than the format of M f , the mass and stiffness matrices are highly dependent on the geometry and topology of the car, and the type of finite elements that are used. Since the structure is essentially modeled with a fine (pretty much uniform) mesh, the matrices have dimensions of several millions. The stiffness matrix depends also on the material parameters, while the damping and coupling matrices depend on geometry, topology and material parameters.

2 Results and discussion

2.1 Direct frequency response

Usually one of the first tasks in the analysis of acoustic fields is to solve the frequency response problem, that is, to simulate the behavior of the acoustic field under excitations (typically of the structure). The classical approach for the frequency response analysis of linear vibrational systems is to perform a Fourier ansatz

[ u d p d ] = [ u ˆ p ˆ ] e ı ω t , f s = f ˆ (ω) e ı ω t ,

which then leads to a frequency dependent linear system of equations

( ω 2 [ M s 0 D s f T M f ] + ı ω [ D s 0 0 D f ] + [ K s ( ω ) D s f T 0 K f ] ) [ u ˆ ( ω ) p ˆ ( ω ) ] = [ f ˆ ( ω ) 0 ] .
(7)

This linear system, which (for a reasonable structure) has several millions of equations, typically has to be solved for a large frequency range ω=0,,1,000 Hz in small frequency steps and, depending on the type of excitation, for many right hand sides (load vectors).

Based on the frequency response analysis it is then possible to detect places where the excitation leads to large noise emissions (hot spots), and this approach can be used to improve or even optimize the frequency response behavior within an optimization loop.

Before we discuss solution methods let us introduce a useful modification of (7). For nonzero frequencies we may multiply the second block row by ω 1 and introduce the new variable v ˆ (ω)= ω 1 p ˆ (ω). Then we obtain a new linear system

( ω 2 [ M s 0 0 M f ] + ı ω [ D s ı D s f T ı D s f T D f ] + [ K s ( ω ) 0 0 K f ] ) [ u ˆ ( ω ) v ˆ ( ω ) ] = [ f ˆ ( ω ) 0 ] ,
(8)

which now has complex symmetric coefficients and, hence allows the use of memory- and flop-saving structured methods. To simplify notation, we write this as

P(ω)x(ω):=( ω 2 M+ıωD+K(ω))x(ω)=b(ω).
(9)

Since this system of equations has to be solved for many right hand sides and over a large frequency range, and all this within an optimization loop, a very efficient solver is necessary. This solver has to be able to recycle information from nearby problems (when stepping through the frequency range and modifying parameters in an optimization) and has to work efficiently in a modern multi-processor, multi-core hardware environment. Altogether, this is really a lot to ask for from a linear system solver and makes the use of black box solvers extremely difficult.

2.1.1 State-of-the-art in linear system solvers

Modern techniques for the solution of linear systems in industrial applications are typically a combination of methods, that use the best of each of the various available classes of methods.

The first class of methods are highly efficient direct solution methods that use Gaussian elimination with partial pivoting combined with other techniques from graph theory to achieve optimal performance, make efficient use of the sparsity structure to avoid fill-in, and save storage. Many of them are even implemented for the use on multiprocessor machines. Well known packages include UMFPACK [7], PARDISO [8], MUMPS [9], WSMP [10], or the HSL collection (for example, HSL_MA78 [11]) to name a few. In view of possibly many right hand sides they are a clear option, despite the fact that for the size and type of problems considered in industrial practice, the storage demands are so high that they typically can only work in an out-of-core setting, that is, the matrix factors are stored on hard disk instead of main memory. Another difficulty is that a new factorization has to be computed for every frequency or every modification in the optimization loop. Finally also the bad scaling and the fact that the problems get increasingly ill-conditioned for large frequencies presents a real challenge, because the desired solution accuracy may not be realizable. Some of the packages provide specialized routines for complex symmetric systems, which have the potential to half the storage and computational requirements, but may also increase the complexity of pivoting strategies.

The second class of methods are iterative methods of the Krylov subspace type like the generalized minimal residual method (GMRES) [12], the Bi-conjugate gradient method in its stabilized form (BICGSTAB) [13], or the quasi-minimal residual method (QMR) [14]. A variant of the conjugate gradient method for complex symmetric systems like (9) is CSYM [15]. In principle these methods would be very good candidates, since the linear systems are sparse and matrix vector multiplications are relatively cheap. Also, if f ˆ and K are independent of ω, then one can exploit the shift invariance property of Krylov subspaces [16, 17]. However, often f ˆ or K do depend on ω and without a good preconditioner the convergence of iterative methods (in particular for large frequencies) is dramatically slow or not realizable.

The third class of potential methods (exploiting the fact that there is a partial differential equation in the background) are multigrid or multilevel methods [1820]. Unfortunately, despite the fact that they are efficiently used in the solution of wave propagation problems, such as Helmholtz or Maxwell equations [21], they cannot be used in a simple way in the described acoustic field problems, since the discrete FEM modeling of the structure does not provide a nice hierarchy of basis functions. If such a hierarchy was available then these methods would lead to very good preconditioners for the Krylov subspace methods. However, with the current state of discrete FEM modeling (being engineering practice for decades) it is extremely difficult to incorporate them into the current code solvers. This is even true for the algebraic multigrid methods like the Ruge-Stüben method [22] or hypre [23], due to the difficulty of having to solve for a large frequency range, many right hand sides and all this inside an optimization loop. It is needless to say that also the incorporation of these solvers in a multiprocessor, multi-core industrial software tool is another challenge.

The fourth class of methods which have made tremendous progress in the last decades are the adaptive finite element methods, see, for example, [24]. These refine the computational grid according to a priori and a posteriori error estimates and if implemented properly avoid globally fine meshes and therefore the extremely large linear systems and eigenvalue problems in the first place. However, the construction, analysis, and implementation of such methods for their use inside an industrial package, including the treatment of a full frequency range and many right hand sides, is still in its infancy and requires a major research effort which fortunately is currently addressed in several research projects world-wide.

When our cooperation with SFE started, we essentially made the above analysis of the available methods and realized the following major obstacles.

  • The problems are badly scaled and get increasingly ill-conditioned when ω grows;

  • for some parameter constellations the system becomes exactly singular with possibly inconsistent right hand sides;

  • typically there are many right hand sides;

  • classical off-the-shelf iterative methods do not work well, or fail completely;

  • direct solution methods have to work by storing the factors in out-of-core memory and cannot easily recycle information from previous frequency or optimization steps;

  • no multilevel or adaptive grid refinement is applicable, the methods must be matrix based;

  • the matrices M, K 1 are often slightly indefinite because of rounding errors during their creation.

2.1.2 What did we do?

With all the obstacles of the problem and all the deficiencies of the current methods, the development and implementation of an appropriate linear solver for acoustic field frequency response problems is a major challenge. Furthermore, as in all industrial projects, there was a deadline. In such a circumstance, the only possibility is to compromise between the optimality of a given method for a given problem, the provability of success, and the practical needs in the industrial environment.

Our compromise, since we had to deal with given matrices (as opposed to partial differential equations), was to develop and implement (in the SFE software environment) a preconditioned subspace recycling Krylov subspace method which has the following features:

The initial search space for any frequency ω is chosen to contain the solutions for the last few frequencies, an approach also used in [25]. For the preconditioner, whose construction and application represent the dominant cost of the algorithm, we had to fulfill the requirement that it had to run in a distributive setting, store its results out-of-core, support complex arithmetic, in particular support complex symmetric systems. This is a tall order to ask for but fortunately we were able to base the preconditioner on existing software, using the complex LD L T factorization code of the MUMPS package [9]. The preconditioner was constructed from the symmetric real part of the linear system (9), that is, from

P ˜ (ω):= ω 2 [ M s 0 0 M f ] + [ K 1 0 0 K f ] .

For small values of ω up to ca 602π corresponding to a frequency of 60 Hz, typically only 3-6 iteration steps per frequency were necessary and the preconditioner could be kept unchanged for the whole frequency range. But the number of iteration steps per frequency grows substantially the larger ω gets, so that more and more new preconditioners have to be computed. For large frequencies close to 1,000 Hz, the preconditioned Krylov method tends to be extremely slow or not convergent at all. As a consequence we constructed a hybrid method, where a certain frequency onwards MUMPS is used as a direct solver for the system with the full matrix P(ω). Furthermore, it was decided to make use of the multi-core environment and to treat several linear systems simultaneously in a distributed fashion. One group of processors solves the systems corresponding to the lower frequencies using the described iterative recycling method. Meanwhile, the other processor groups treat the high frequency systems with a full direct solve for each value of ω.

In Figure 2 we present computation times of our solver for a model with 551,388 degrees of freedom, solved on a 4 socket machine equipped with Intel Xeon X5670 processors (resulting in 48 virtual CPUs) using 7.9 GB of RAM under Linux 2.6.34. All factorizations were performed out-of-core. It can be observed that computation times scale almost perfectly with the number of systems treated concurrently and that an increase in the number of right hand sides causes only a mild increase in computation time. Moreover, times can be greatly reduced by distributing the solution (LD L T decomposition) of every individual system on multiple processors, although in this case scaling is not perfect.

Figure 2
figure 2

Computation times of the described solver for 1, 2, 4, and 8 concurrent freqencies with different number of processors per frequency (×1 or ×6) and different number of right hand sides (RHS).

2.1.3 Evaluation of our approach

Since there are commercial packages available, a question that may arise is why we need new numerical methods and solvers in the first place. From the point of view of the company SFE this is clear, they wanted their own solver in their product, and the solution method including the implementation was satisfactory for their needs. But in many respects the current methods are unsatisfactory from a scientific point of view. The discussed industrial problems present some of the current grant challenges in the field and the commercial packages are by no means able to solve these problems in a completely satisfactory way. To always get a more-or-less reasonable solution, and to be efficient, they often compromise for accuracy of the results which is certainly inadequate from a scientific point of view. Thus the cooperation with the industrial partner triggered new research questions for the academic world.

It became clear during the project that the concept of recycling in iterative methods, that is, the reuse of information that was already computed for other frequencies or in the course of an optimization procedure is not well-enough understood. As a consequence, motivated by the project with SFE, a new research project [26] in the DFG Research Center Matheon was started to further investigate the basic mathematical principles and to understand how they can be incorporated into new efficient methods [27].

As a second major obstacle, we identified the use of discrete FEM modeling in structural engineering. It would be much easier if adaptive FEM would be usable for the discussed class of problems, and it would also be good to have a grid hierarchy that allows the use of efficient multilevel preconditioners [1820]. Despite the high research activity in this field this is a major challenge for the described problem classes and almost no analysis or methods are available. We will discuss this in more detail in the section devoted to eigenvalue computation, but again here we see a need for a stronger cooperation between academia and industry in this area of structural engineering, to transfer new ideas that are developed now into the industrial practice.

2.2 Modal reduction

We have seen in the previous section that due to the fine mesh used for the discrete FEM modeling of the structure, the linear systems have very large dimensions. Typically, however, one is interested in damping the low frequencies associated with the eigenvalues in the neighborhood of 0 (and the imaginary axis) of the complex symmetric matrix function

Q(λ):= λ 2 M+λD+K=Q ( λ ) T ,
(10)

with M D K as in (9). The methods discussed below assume that P is quadratic in λ for the eigenvalue computation, that is, the nonlinear dependency of K on the frequency is ignored. One easily verifies that if λ 0 is an eigenvalue of Q(λ), that is, P( λ 0 )x=0, then x T is the corresponding left eigenvector but essentially no other general properties of the eigenvalue problem are available, in particular, there is unfortunately no immediate variational property for the eigenvalues/eigenfunctions available, as there is for the undamped case [28].

Since one is interested only in part of the spectrum, a natural idea is to identify a space associated with eigenvectors and generalized eigenvectors associated with the important eigenvalues, and to project the problem into this subspace. This is a model reduction approach called modal reduction which, if efficiently implemented, can save a huge amount of computing time and storage, despite the fact that it is partially heuristic.

2.2.1 Complex symmetric quadratic eigenvalue problems

In view of the desire to do modal reduction, another task in the cooperation with SFE GmbH was the computation of a small number, say , eigenvalues in a trapezoidal region near zero, see Figure 3(left), and the corresponding subspace S spanned by the eigenvectors and generalized eigenvectors associated with these eigenvalues, that is, S =span{ s 1 ,, s }, where the s i form an orthonormal basis of this subspace, that is, S =[ s 1 ,, s ] is an isometry.

Figure 3
figure 3

A typical trapezoidal region - within which all eigenvalues are sought - at beginning and after three shifts have been processed.

The projected system would have the form

Q (λ):= λ 2 M +λ D + K := λ 2 S T M S +λ S T D S + S T K S
(11)

and would still be complex symmetric. In the context of model reduction, the requirements for the reduced model are that the projected system is a good approximation to the large scale system for a large frequency range, and also for a large set of parameter variations. Furthermore, it has to be of small enough dimension, so that classical methods for nonlinear non-smooth optimization can be applied. Again this is a lot to ask, since currently for large scale problems there are really no methods available that guarantee to obtain all the eigenvalues of (10) in a specified region of the complex plane. For small dense problems one could employ the sign function method, [29, 30] but this would require storing full dense inverses of matrices of the given class, which is certainly not possible in the described acoustic field problems.

The currently used industrial techniques typically solve a simplified problem, for example, the eigenvalues and associated eigenvectors in the desired region for the undamped/uncoupled problem are used for the projection or an algebraic multilevel substructuring method (AMLS) is used that exploits the structure of the matrices as they arise form the discrete FEM [31]. All these techniques are partially heuristic, since there is no guarantee that all the desired eigenvalues are captured or that the eigenvalues from the projected system are close to those of the original physical system or of the FEM discretized system, that is, in general, currently no error bounds are available. To generate such error bounds is a major challenge for the academic community and research in this direction has started in the project [32] with first results in [33, 34].

Furthermore, the numerical solution of the large scale nonlinear eigenvalue problem (10) itself also presents a mathematical challenge, for several reasons.

First of all the mass matrix is singular, so the problem has eigenvalues at ∞ which often cause convergence problems for the iterative methods, second of all, there are also several (typically six) eigenvalues at 0, corresponding to the six free degrees of motion, three each for translation and rotation, of the whole structure. In industrial practice the eigensolver is furthermore also often used for model verification, by checking whether there are exactly 6 eigenvalues at 0. Then, if the model is flawed, the singularity may be even higher.

The major challenge, however, is that we want the eigenvalues near 0 and it is well known that classical iterative methods for large sparse eigenvalue problems, like the implicitly restarted Arnoldi method [35], typically converge fast only to the eigenvalues at the periphery of the spectrum [36]. Thus, either a shift-and-invert technique that transforms the problem and maps the desired eigenvalue to the periphery is necessary or other techniques like Newton’s method [37, 38] or the Jacobi-Davidson method [39, 40] have to be used. All these methods require the solution of linear systems of the form ( λ ˆ 2 M+ λ ˆ D+K)x=b with a given shift-point λ ˆ . But this is the problem that we wanted to avoid in the first place and hence a vicious circle is closed.

2.2.2 What did we do?

Having assessed the various options and their potential advantages and drawbacks, we decided to develop a new method based on the following concepts and to implement it into the SFE environment.

First of all, since the methods that work directly on the quadratic eigenvalue problems, for example, [28, 41, 42], are not yet as mature as those for linear eigenvalue problems, the problem is linearized by introducing a new variable λx and turning (10) from a quadratic into a linear eigenvalue problem.

Since there is no method available that is guaranteed to find all eigenvalues in a region at feasible costs, we had to revert to a partially heuristic approach. We are looking for eigenvalues in the vicinity of target points which are scattered inside the region of interest. To find eigenvalues close to a particular target σ, also called a shift, we use the shift-linearize-invert-ansatz [42] that requires computation of the largest eigenvalues of the matrix

A:= [ Q ( σ ) 1 0 0 I ] [ 2 σ M + D M I 0 ] .

This is done by the block Arnoldi method [43] which requires the application of the matrix A to given blocks of vectors during the generation of the Krylov subspace. This involves a solution of a linear system with the complex symmetric matrix Q(σ). We again use the direct solver MUMPS to compute a complex LD L T factorization.

The block Arnoldi method searches for eigenvectors of A within the block Krylov subspace

K m (A,B):=span{B,AB, A 2 B,, A m 1 B},
(12)

which is known to contain, for increasing m, increasingly good approximations to eigenvectors corresponding to eigenvalues of A of maximum absolute value - which correspond to the eigenvalues of Q(λ) closest to σ. In the classical Arnoldi method, the starting block B is a vector. Otherwise, when B consists of several, say n b , columns, the resulting method is the block Arnoldi method. The matrix B can be chosen in an (almost) arbitrary way, but if eigenvector or subspace approximations are known, the use of these will speed up convergence. Among the advantages of using the block method are that clusters of eigenvalues are handled better [44]. Furthermore, all vector-vector operations become matrix-matrix operations which can be implemented much more efficiently by use of BLAS level 3 routines [45] and the matrix-vector multiplication with the large matrix A can be carried out with a block of vectors. Furthermore, the heuristic part of the algorithm is more reliable in the block case, see below.

In the actual implementation of the block Arnoldi algorithm the terms A i B are never explicitly formed. Instead an orthonormal basis Q m =[ V 1 , V 2 ,, V m ] of K m (A,B) is generated by setting V 1 to be an orthonormal basis of span(B) followed by iteratively forming V i + 1 by orthonormalization of A V i against V 1 ,, V i . Classically, a variant of the Gram-Schmidt method is used for this task [35], but any other orthonormalization procedure may be employed. We use block Householder reflections [46, 47] that are very stable and attain high efficiency by using BLAS-3 operations. Collecting all the V i and the orthonormalization coefficients H i , j results in the well known Arnoldi relation A Q m = Q m H m + V i + 1 H m + 1 , m E m T with an (m×m)-block Hessenberg matrix H m = ( H i , j ) i , j = 1 m . The eigenvalues of H m , called Ritz values, are used as approximations of eigenvalues of A. Likewise, the eigenvectors of H m , multiplied with Q m , are used as approximate eigenvectors of A.

A drawback of the block Arnoldi method is the necessity to store the basis Q m which grows for increasing m. A typical way out of this problem is to restart the algorithm at some point. Instead of restarting with a single new starting block it is possible to restart with a whole Arnoldi relation. For the vector Arnoldi method, two such implicit restarting schemes are common, one using a filter polynomial [35] and one using a reordering of the Schur form of H m [48]. We are using the latter approach as the first one does not elegantly generalize to the block case [46].

The heuristic part of our method is that we assume that the block Arnoldi method really finds all eigenvalues in a circle around the shift σ. While quite often the eigenvalues do indeed appear and converge in the order of the distance to the shift, it is not rare that one or a group of eigenvalues converge slower than other farer away eigenvalues. However, in these cases usually the missed eigenvalues are present as (yet) unconverged Ritz values. Therefore, we use the unconverged Ritz value that is closest to the shift as the radius of a circle that is trusted to contain no missed eigenvalues.

In the implementation, we repeatedly run the block Arnoldi method for different shifts, possibly several at once in a distributed setting. Figure 3(right) shows the situation after a few iterations. Large parts of the trapezoidal region are covered, leaving only some small remaining regions to be searched. New shifts are placed inside the largest such white regions, until the whole trapezoidal region of interest is covered by trusted circles.

Of course, with such a covering approach an eigenpair could be computed by more than one Arnoldi run for different shifts. For that reason the freshly discovered eigenpairs have to be checked for being copies of already previously found pairs. To achieve this we consider a new eigenpair ( λ , x ) being a copy, if [ x T , λ x T ] T is almost linearly dependent to the span of the vectors [ x old T , λ old x old T ] T corresponding to every known eigenvalue λ old sufficiently close to λ .

We applied the solver to a car model, discretized by a regular mesh of 35 mm leading to 219,432 degrees of freedom. We were looking for eigenvalues within the triangle bordered by the lines Im(λ)>20Re(λ), Im(λ)>20Re(λ), and Im(λ)<f2π, for f{50,100,150,200,250}. Table 1 lists the number of eigenvalues within these triangles, the number of used shifts, that is, used matrix factorizations, and the number of overall Arnoldi iterations, that is, matrix-block products. The test was performed on a PC with an Intel Core2 Duo E6850 CPU clocked at 3.0 GHz, with 4 GB RAM. One shift was addressed at a time using one processor. The block size was 5.

Table 1 Number of required shifts and iterations to find all eigenvalues between 0 and {50,,250} Hz.

2.2.3 Evaluation of our approach

Our approach certainly does not use many novel ideas, but instead builds on mature and proven concepts. The method can run in a distributed setting processing several shifts at once, each one running in several processes. The matrix factorizations can be kept out-of-core. Moreover, since the basis vectors are explicitly orthogonalized by a stable Householder scheme, the generated block Arnoldi basis is orthonormal to working precision. This is in sharp contrast to the basis generated by the non-symmetric Lanczos method which typically looses linear independence after enough iterations and leads to so-called spurious eigenvalues [49, 50].

In our experiments the heuristic choice of the radii of trusted circles worked very well. In thousands of test cases with problems from SFE as well as randomly generated problems, it happened only once that this approach missed an eigenvalue. This happened with a block size of one, that is, with the vector Arnoldi method. The method worked fine for this case if a block size of at lest two was used. In our implementation the default block size is 8.

Unfortunately, the improved robustness of the block-Arnoldi method compared with the standard Arnoldi method comes at the price of increased memory requirements to store the basis. This downside is somewhat mitigated, however, by the use of restarts.

As another downside of our approach, so far only quadratic eigenvalue problems can be solved, for truly nonlinear problems with ω dependence in K the methods is not directly applicable. On the other hand, the method can be applied to fluid or structure subsystems separately, or to the complete coupled system, and it tolerates a singular mass matrix.

3 Conclusions. The two-way street of industrial cooperation

We discussed the modeling, simulation and (at least the vision) of the optimization of acoustic fields. As an example of our industrial cooperation we studied frequency domain computations and modal reduction methods for the acoustic field inside a car, as well as the challenges of the resulting linear systems and eigenvalue problems. One lesson from the project was that the methods to be used in industrial practice cannot be constructed and implemented using textbook approaches. For instance, very often tricks have to be employed that increase the efficiency of the computation, but that lack a full mathematical understanding. This can lead to misunderstandings between engineers, programmers and mathematicians. A stronger communication and cooperation between these groups is necessary to address the described challenges. If this is achieved then all sides benefit from a cooperation. The work with SFE GmbH on the frequency response problem for interior acoustic field computation started out as a clear transfer project (a one way street), with the idea to transport the knowledge and know-how about current linear system and eigenvalue solver technologies available in the academic environment into an industrial software environment.

But as we have discussed, already early on in the cooperation a lot of new research topics appeared that could not be treated within the current project (two-way-street). Examples include recycling methods for a sequence of slowly changing linear systems, updating of preconditioners, guaranteed location of all eigenvalues inside a region of the complex plane, to name a few.

References

  1. ANSYS, Inc.: Programmer’s manual for ANSYS. (2007)[http://www1.ansys.com/customer/content/documentation/110/ansys/aprog110.pdf]

  2. MSC. Software Corporation: MSC.Nastran 2005 Installation and Operations Guide. (2004)[http://www.mscsoftware.com/Products/CAE-Tools/MSC-Nastran.aspx]

  3. Antoulas AC: Approximation of large-scale dynamical systems. In Advances in Design and Control. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; 2005.

    Google Scholar 

  4. Benner P., Mehrmann V., Sorensen D. (Eds): Dimension Reduction of Large-Scale Systems. Springer, Heidelberg; 2005.

    MATH  Google Scholar 

  5. Schmidt F, Deuflhard P: Discrete transparent boundary conditions for the numerical solution of Fresnel’s equation. Comput. Math. Appl. 1995, 29: 53–76. 10.1016/0898-1221(95)00037-Y

    Article  MATH  MathSciNet  Google Scholar 

  6. Bishop R: The treatment of damping forces in vibration theory. Jl R. Aeronaut. Soc. 1955, 59: 738.

    Google Scholar 

  7. Davis TA: Algorithm 832: UMFPACK V4.3 - an unsymmetric-pattern multifrontal method. ACM Trans. Math. Softw. 2004,30(2):196–199. 10.1145/992200.992206

    Article  MATH  Google Scholar 

  8. Schenk O, Gärtner K: Solving unsymmetric sparse systems of linear equations with PARDISO. Future Gener. Comput. Syst. 2004,20(3):475–487. [http://www.sciencedirect.com/science/article/B6V06–49NXY7J-F/2/e8260e5d8f19639019cddea4776c024c] 10.1016/j.future.2003.07.011

    Article  Google Scholar 

  9. MUMPS: MUltifrontal Massively Parallel Solver: Users’ Guide, (2009) http://graal.ens-lyon.fr/MUMPS

  10. Gupta, A.: WSMP: Watson Sparse Matrix Package, Part II. Direct solution of general sparse systems. IBM research report, IBM T. J. Watson Research Center. http://www.alphaworks.ibm.com/tech/wsmp (2000)

    Google Scholar 

  11. Reid, J.K., Scott, J.A.: An efficient out-of-core multifrontal solver for large-scale unsymmetric element problems. Technical Report RAL-TR-2007-014, SFTC Rutherford Appleton Laboratory (2007)

    Google Scholar 

  12. Saad Y, Schultz MH: GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 1986,7(3):856–869. 10.1137/0907058

    Article  MATH  MathSciNet  Google Scholar 

  13. van der Vorst HA: Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 1992,13(2):631–644. 10.1137/0913035

    Article  MATH  MathSciNet  Google Scholar 

  14. Freund RW, Nachtigal NM: QMR: a quasi-minimal residual method for non-Hermitian linear systems. Numer. Math. 1991,60(3):315–339.

    Article  MATH  MathSciNet  Google Scholar 

  15. Bunse-Gerstner A, Stöver R: On a conjugate gradient-type method for solving complex symmetric linear systems. Linear Algebra Appl. 1999,287(1–3):105–123. 10.1016/S0024-3795(98)10091-5

    Article  MATH  MathSciNet  Google Scholar 

  16. Gu GD, Simoncini V: Numerical solution of parameter-dependent linear systems. Numer. Linear Algebra Appl. 2005,12(9):923–940. 10.1002/nla.442

    Article  MATH  MathSciNet  Google Scholar 

  17. Simoncini V, Perotti F: On the numerical solution of $(\lambda^{2}A+\lambda B+C)x=b$(λ2A+λB+C)x=b and application to structural dynamics. SIAM J. Sci. Comput. 2002,23(6):1875–1897. (electronic) (electronic) 10.1137/S1064827501383373

    Article  MATH  MathSciNet  Google Scholar 

  18. Hackbusch W: Multigrid Methods and Applications. Springer-Verlag, Berlin; 1985.

    Book  Google Scholar 

  19. McCormick SF: Multilevel Adaptive Methods for Partial Differential Equations. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; 1989.

    Book  MATH  Google Scholar 

  20. McCormick SF: Multilevel Projection Methods for Partial Differential Equations. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; 1992.

    Book  MATH  Google Scholar 

  21. Ernst, O.G.: Fast numerical solution of exterior Helmholtz problems with radiation boundary condition by imbedding. Ph.D. thesis, Stanford University (1994)

    Google Scholar 

  22. Ruge JW, Stüben K: Algebraic multigrid. In Multigrid Methods. SIAM, Philadelphia, PA; 1987:73–130.

    Chapter  Google Scholar 

  23. Falgout RD, Jones JE, Yang UM: The design and implementation of hypre, a library of parallel high performance preconditioners. In Numerical Solution of Partial Differential Equations on Parallel Computers. Springer, Berlin; 2006:267–294.

    Chapter  Google Scholar 

  24. Bangerth W, Rannacher R: Adaptive Finite Element Methods for Differential Equations. Birkhäuser, Basel, Ch.; 2003.

    Book  MATH  Google Scholar 

  25. Fischer PF: Projection techniques for iterative solution of $A\underline{x}=\underline{b}$Ax̲=b̲ with successive right-hand sides. Comput. Methods Appl. Mech. Eng. 1998,163(1–4):193–204. 10.1016/S0045-7825(98)00012-7

    Article  MATH  Google Scholar 

  26. Numerical methods for large-scale parameter-dependent systems. Project C29, DFG Research Center Matheon, Berlin http://www.matheon.de/research/show_project.asp?id=172

  27. Gaul, A., Gutknecht, M.H., Liesen, J., Nabben, R.: Deflated and augmented Krylov subspace methods: Basic facts and a breakdown-free deflated MINRES. Preprint 749, DFG Research Center Matheon, Mathematics for key technologies, Berlin (2011).

    Google Scholar 

  28. Mehrmann V, Voss H: Nonlinear eigenvalue problems: a challenge for modern eigenvalue methods. Mitt. der Ges. f. Angewandte Mathematik und Mechanik 2005, 27: 121–151.

    MathSciNet  Google Scholar 

  29. Bai Z, Demmel J: Using the matrix sign function to compute invariant subspaces. SIAM J. Matrix Anal. Appl. 1998, 19: 205–225. (electronic) (electronic) 10.1137/S0895479896297719

    Article  MATH  MathSciNet  Google Scholar 

  30. Higham NJ: Functions of matrices. Theory and Computation. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; 2008.

    Book  MATH  Google Scholar 

  31. Bennighof JK: Adaptive multi-level substructuring method for acoustic radiation and scattering from complex structures. In Computational Methods for Fluid/Structure Interaction. Edited by: Kalinowski A.J.. AMSE, New York; 1993:25–38.

    Google Scholar 

  32. Adaptive solution of parametric eigenvalue problems for partial differential equations. Project C22, DFG research center Matheon, Berlin.[http://www.matheon.de/research/show_project.asp?id=126]

  33. Carstensen, C., Gedicke, J., Mehrmann, V., Miedlar, A.: An adaptive homotopy approach for non-selfadjoint eigenvalue problems. Preprint 718, DFG Research Center Matheon, Mathematics for key technologies, Berlin (2010)

    Google Scholar 

  34. Miedlar, A.: Inexact adaptive finite element methods for elliptic PDE eigenvalue problems. Ph.D. thesis, Technische Universität Berlin, Inst. f. Mathematik (2011)

    Google Scholar 

  35. Lehoucq RB, Sorensen DC, Yang C: ARPACK Users’ Guide. Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; 1998.

    Book  Google Scholar 

  36. Saad Y: Numerical Methods for Large Eigenvalue Problems. Manchester University Press, Manchester; 1992.

    MATH  Google Scholar 

  37. Schwetlick, H., Schreiber, K.: A primal-dual Jacobi-Davidson-like method for nonlinear eigenvalue problems. Internal Report ZIH-IR-0613, Techn. Univ. Dresden, Zentrum für Informationsdienste und Hochleistungsrechnen (2006)

    Google Scholar 

  38. Schwetlick, H., Schreiber, K.: Nonlinear Rayleigh functionals. Linear Algebra and Its Applications. In Press, Corrected Proof (2010). doi:10.1016/j.laa.2010.06.048

    Google Scholar 

  39. Sleijpen GLG, Booten AGL, Fokkema DR, Van der Vorst HA: Jacobi-Davidson type methods for generalized eigenproblems and polynomial eigenproblems. BIT Numer. Math. 1996,36(3):595–633. 10.1007/BF01731936

    Article  MATH  MathSciNet  Google Scholar 

  40. Sleijpen GLG, Van der Vorst HA: A Jacobi-Davidson iteration method for linear eigenvalue problems. SIAM J. Matrix Anal. Appl. 1996,17(2):401–425. 10.1137/S0895479894270427

    Article  MATH  MathSciNet  Google Scholar 

  41. Beyn, W.: An integral method for solving nonlinear eigenvalue problems. ArXiv e-prints (2010) Beyn, W.: An integral method for solving nonlinear eigenvalue problems. ArXiv e-prints (2010)

  42. Tisseur F, Meerbergen K: The quadratic eigenvalue problem. SIAM Rev. 2001,43(2):235–286. 10.1137/S0036144500381988

    Article  MATH  MathSciNet  Google Scholar 

  43. Arnoldi WE: The principle of minimized iteration in the solution of the matrix eigenvalue problem. Q. Appl. Math. 1951, 9: 17–29.

    MATH  MathSciNet  Google Scholar 

  44. Lehoucq, R., Maschhoff, K.: Implementation of an implicitly restarted block Arnoldi method. Preprint MCS-P649–0297, Argonne National Laboratory, Argonne, IL (1997)

    Google Scholar 

  45. Dongarra JJ, Croz JD, Hammarling S, Duff I: A set of level 3 basic linear algebra subprograms. ACM Trans. Math. Softw. 1990, 16: 1–17. 10.1145/77626.79170

    Article  MATH  Google Scholar 

  46. Baglama J: Augmented block Householder Arnoldi method. Linear Algebra Appl. 2008,429(10):2315–2334. 10.1016/j.laa.2007.12.021

    Article  MATH  MathSciNet  Google Scholar 

  47. Schreiber R, Van Loan C: A storage-efficient WY representation for products of Householder transformations. SIAM J. Sci. Stat. Comput. 1989, 10: 53–57. 10.1137/0910005

    Article  MATH  MathSciNet  Google Scholar 

  48. Stewart GW: A Krylov-Schur algorithm for large eigenproblems. SIAM J. Matrix Anal. Appl. 2001/02,23(3):601–614.

    Article  MathSciNet  Google Scholar 

  49. Cullum J, Willoughby RA: Lanczos and the computation in specified intervals of the spectrum of large, sparse real symmetric matrices. In Sparse Matrix Proceedings 1978 (Sympos. Sparse Matrix Comput., Knoxville, Tenn., 1978). SIAM, Philadelphia, Pa; 1979:220–255.

    Google Scholar 

  50. Freund RW, Gutknecht MH, Nachtigal NM: An implementation of the look-ahead Lanczos algorithm for non-Hermitian matrices. SIAM J. Sci. Comput. 1993, 14: 137–158. 10.1137/0914009

    Article  MATH  MathSciNet  Google Scholar 

Download references

Acknowledgements

We thank H. Zimmer and A. Hilliges from our industrial partner SFE for the permission to use the figures obtained in the cooperation, and we thank Tobias Brüll, Elena Teidelt, Matthias Miltenberger for their help in the implementation. Moreover, we thank two anonymous referees for their valuable comments on a preliminary version of this paper.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Volker Mehrmann.

Additional information

Competing interests

Our research was supported by Deutsche Forschungsgemeinschaft, via the DFG Research Center Matheon, Mathematics for Key Technologies, in Berlin and by SFE GmbH, Voltastr. 5, 13355 Berlin.

Authors’ contributions

All authors participated in designing and implementing the solvers as well as in writing the manuscript. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Mehrmann, V., Schröder, C. Nonlinear eigenvalue and frequency response problems in industrial practice. J.Math.Industry 1, 7 (2011). https://doi.org/10.1186/2190-5983-1-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/2190-5983-1-7

Keywords