Rayleigh–Ritz method

Description of method

The Rayleigh–Ritz method allows for the computation of Ritz pairs ( λ ~ i , x ~ i ) which approximate the solutions to the eigenvalue problem

Where A ∈ C N × N .

The procedure is as follows:

1. Compute an orthonormal basis V ∈ C N × m approximating the eigenspace corresponding to m eigenvectors
2. Compute R ← V ∗ A V
3. Compute the eigenvalues of R solving R v i = λ ~ i v i
4. Form the ritz pairs ( λ ~ i , x ~ i ) = ( λ ~ i , V v i )

One can always compute the accuracy of such an approximation via ∥ A x ~ i − λ ~ i x ~ i ∥

If a Krylov subspace is used and A is a general matrix, then this is the Arnoldi Algorithm.

The method in calculus of variation

In this technique, we approximate the variational problem and end up with a finite dimensional problem. So let us start with the problem of seeking a function y = y0(x) that extremizes an integral I(y). Assume that we are able to approximate y(x) by a linear combination of certain linearly independent functions of the type:

y(x) ≈ φ0(x) + c1φ1(x) + c2φ2(x) + ··· + cN φN (x) (1)

where we will need to determine the constant coefficients c1, ··· cN .

The selection of which approximating function φ(x) to use is arbitrary except for the following considerations:

a) If the problem has boundary conditions such as fixed end points, then φ0(x) is chosen to satisfy the problem’s boundary conditions, and all other φivanish at the boundary. (This should remind the reader of the method of eigenfunction expansion for inhomogeneous partial differential equation. The functions φi are not necessarily eigenfunctions though.)

b) In those problems where one knows something about the form of the solution then the functions φi(x) can be chosen so that the expression (1) will have that form.

By using (1) we essentially replace the variational problem of finding an arc y(x) that extremizes I to finding a set of constants c1, ··· , cN that extremizes I(c1, c2, ··· , cN ). We solve this problem as a finite dimensional one i. e. by solving

∂I ∂ci = 0, i = 1, ··· , N (2)

The procedure is to first determine an initial estimate of c1 by the approximation y ≈ φ0+c1φ1. Next, the approximation y ≈ φ0+c1φ1+c2φ2 is used (with c1 being redetermined). The process continues with y ≈ φ0 + c1φ1 + c2φ2 + c3φ3 as the third approximation and so on. At each stage the following two items are true:

a) At the i th stage, the terms c1 ··· , ci−1 that have been previously determined are redetermined

b) The approximation at the i th stage

y ≈ φ0 + c1φ1 + ··· + ciφi (3)

will be better or at least no worse than the approximation at the i − 1st stage

y ≈ φ0 + c1φ1 + ··· + ci−1φi−1 (4)

Convergence of the procedure means that

lim N→∞ � φ0 + � N i=1 ciφi � = y0(x) (5)

where y0(x) is the extremizing function.

In many cases one uses a complete set of functions e. g. polynomials or sines and cosines. A set of functions φi(x) (i = 1, 2, ···) is called complete over [a, b] if for each Riemann integrable∗∗ function f(x), there is a number N (depending on , c1, ··· , cN ) such that

max [a,b] [f − � N i=1 ciφi] 2 < (6)

The above outlined procedure can be extended in a number of ways. For example, more than one independent variable may be involved. So for the problem of

min I = � � R F(x, y, w, wx, wy)dydx (7)

subject to w = h(s) on the boundary Γ of R (8)

where h(s) is some prescribed function and s is the arc length along Γ. Analogous to (1) we write

w(x, y) = φ0(x, y) + c1φ1(x, y) + ··· + cN φN (x, y) (9)

and φ0 satisfies (8) and φi(x) i = 1, 2, 3 ··· are zero on Γ. We could also extend the procedure to functions involving higher derivatives, more independent variables, etc.

Example: Mechanical Engineering

Typically in mechanical engineering it is used for finding the approximate real resonant frequencies of multi degree of freedom systems, such as spring mass systems or flywheels on a shaft with varying cross section. It is an extension of Rayleigh's method. It can also be used for finding buckling loads and post-buckling behaviour for columns.

The following discussion uses the simplest case, where the system has two lumped springs and two lumped masses, and only two mode shapes are assumed. Hence M = [m₁, m₂] and K = [k₁, k₂].

A mode shape is assumed for the system, with two terms, one of which is weighted by a factor B, e.g. Y = [1, 1] + B[1, −1]. Simple harmonic motion theory says that the velocity at the time when deflection is zero, is the angular frequency ω times the deflection (y) at time of maximum deflection. In this example the kinetic energy (KE) for each mass is 1 2 ω 2 Y 1 2 m 1 etc., and the potential energy (PE) for each spring is 1 2 k 1 Y 1 2 etc. For continuous systems the expressions are more complex.

We also know, since no damping is assumed, that KE when y=0 equals the PE when v=0 for the whole system. As there is no damping all locations reach v=0 simultaneously.

so, since KE = PE,

Note that the overall amplitude of the mode shape cancels out from each side, always. That is, the actual size of the assumed deflection does not matter, just the mode shape.

Mathematical manipulations then obtain an expression for ω , in terms of B, which can be differentiated with respect to B, to find the minimum, i.e. when d ω / d B = 0 . This gives the value of B for which ω is lowest. This is an upper bound solution for ω if ω is hoped to be the predicted fundamental frequency of the system because the mode shape is assumed, but we have found the lowest value of that upper bound, given our assumptions, because B is used to find the optimal 'mix' of the two assumed mode shape functions.

There are many tricks with this method, the most important is to try and choose realistic assumed mode shapes. For example, in the case of beam deflection problems it is wise to use a deformed shape that is analytically similar to the expected solution. A quartic may fit most of the easy problems of simply linked beams even if the order of the deformed solution may be lower. The springs and masses do not have to be discrete, they can be continuous (or a mixture), and this method can be easily used in a spreadsheet to find the natural frequencies of quite complex distributed systems, if you can describe the distributed KE and PE terms easily, or else break the continuous elements up into discrete parts.

This method could be used iteratively, adding additional mode shapes to the previous best solution, or you can build up a long expression with many Bs and many mode shapes, and then differentiate them partially.