Modal analysis is a set of mathematical techniques for determining the modal characteristics of a dynamical system. Modal analysis can be used to determine the resonant frequencies and vibration responses of a structure and is closely related to failure analysis. The technique often assumes a system to behave linearly as an approximation so that matrix algebra can be used, for which computers are especially well suited. This assumption turns out to be pretty good, in fact, as unperson points out, since stable equilibrium configurations correspond to local minima in potential energy, which can be approximated locally by parabolas, the hallmarks of the simple harmonic oscillator.

Consider three masses m1 = m, m2 = M, and m3 = m constrained to move along a single direction, say, the x axis, with m1 and m3 on either side of m2. Now join m1 and m3 separately to m2 with two springs having spring constants equal to k. Let the three generalized coordinates of the masses be x1, x2 and x3, each of which are zero at the equilibrium positions of the masses.

The total kinetic energy of the system is

T = ½m(dx1/dt)2 + ½M(dx2/dt)2 + ½m(dx3/dt)2

From Hooke's law, the total potential energy of the system is

V = ½k(x2-x1)2 + ½k(x3-x2)2

V = ½kx12 + kx22 + ½kx32 - kx1x2 - kx2x3

So the Lagrangian for the system is L = T - V:

L = ½m(dx1/dt)2 + ½M(dx2/dt)2 + ½m(dx3/dt)2 - ½kx12 - kx22 - ½kx32 + kx1x2 + kx2x3

Now, assuming there aren't any damping forces present, we can obtain the equations of motion from

d/dt (∂L/∂x'i) - ∂L/∂xi = 0

Where x'i is taken to be dxi/dt. So,

∂L/∂x'1 = mx1
∂L/∂x1 = kx1 - kx2

∂L/∂x'2 = Mx2
∂L/∂x2 = 2kx2 - kx1 - kx3

∂L/∂x'3 = mx3
∂L/∂x3 = kx3 - kx2

This results in a set of three equations:

md2x1/dt2 - kx1 + kx2 = 0
Md2x2/dt2 - 2kx2 + kx1 + kx3 = 0
md2x3/dt2 - kx3 + kx2 = 0

Now, we let x = [x1, x2, x3]T and x" = [d2x1/dt2, d2x2/dt2, d2x3/dt2]T, and define the following matrices:

```    /m 0 0\
M = |0 M 0|
\0 0 m/

/ k  -k  0\
K = |-k  2k -k|
\ 0  -k  k/
```

Here M is the mass matrix and K is the stiffness matrix, although in general these matrices may be in units other than that of mass or stiffness. With these definitions, we may write our equations of motion as:

Mx" - Kx = 0

Now we may make a series of substitutions in order to simplify calculations. First, we define the matrix square root. The square root of a matrix A, if defined, is A½ such that A½ A½ = A. A convenient way of finding the square root of a matrix is by noting that if A is diagonalizable by an orthogonal matrix P, such that PTAP = D where D is a diagonal matrix, then A½ = PD½PT. This is because PPT = PTP = I and PD½PTPD½PT = PD½D½PT = PDPT = A.

For those unfamiliar with the diagonalization process, D is a diagonal matrix with diagonal entries being the eigenvalues of A. The columns of P are the eigenvectors of A corresponding to the eigenvalues of A. The columns must be in the same order as their corresponding eigenvectors in D.

Now D½ is easy to calculate because only the diagonal terms in D are present, so we take their square roots individually. Finally, let D be the inverse of D½.

With these in place, let x = Mq, or equivalently let q = M½x. In our example. We do not need to use diagonalization to calculate M because M is already a diagonal matrix. Therefore substitute into the equation of motion to obtain

MMq" - KMq = 0

Left multiply by M to obtain

MMMq" - MKMq = 0

q" - MKMq = 0

Now, let K = MKM, which is the modified stiffness matrix. At this point we obtain

q" - Kq = 0

This is an eigenvalue problem and can be solved directly, but to further simplify the problem, we can turn use some more substitution. First, we must diagonalize K, i.e. find a matrix P, such that PTKP = D where D is a diagonal matrix. Then let q = Pr so that

Pr" - KPr = 0

This time, left multiply by PT so that:

PTPr" - PTKPr = 0

r" - Dr = 0

With the problem simplified thusly, we can easily express the solutions to the equations of motion, since by changing coordinates from xi to qi and then to ri, we've managed to decouple the systems completely (D is diagonal), so that we may write:

r"i - λi ri = 0

Writing λi = ωi2, where λi is the ith eigenvalue of K and diagonal element of D, and where { ωi } is the set of natural or resonant frequencies of the system. We can solve each of these equations independently, so that ri = Ai cos(ωit + φi). Any general set of coordinates in which the system oscillates at a single frequency in each coordinate is called a set of normal coordinates. Thus {ri} represent the normal coordinates for the system. We can get back our original coordinates by applying the transformations in reverse:

x = Mq
q = Pr

So,

x = MPr

The matrix S = MP is called the matrix of modal shapes, since it encodes characteristic shape of each mode. Writing x = Sr, and denoting the column vectors of S as Si we get that

x = A1S1cos(ω1t + φ1) + A2S2cos(ω2t + φ2) + A3S3cos(ω3t + φ3) + ...

Thus, the net displacement of each mass, the components of x, are linear superpositions of the displacements brought about by each mode of oscillation. For instance, in our example of three masses connected to eachother with springs, we would find that one mode for which S1 = [ -1,0,1 ] T. This mode is antisymmetric with the outer masses oscillating 90° out of phase with one another, and the inner mass remaining stationary. This mode has its own frequency ω1 and is independent of and uninfluenced by the other modes.

In order to solve the problem in its entirety, then, we need determine the values of Ai and φi. This can be done as long as the initial conditions for the system are given, or any other independent conditions that fix these constants. Normally the initial conditions are given in terms of the system's initial displacements and velocities.

In the above development we used an approach that allowed us to find both the modal shapes and natural frequencies. If we were only after the frequencies ωi we would have used the equation Mx" - Kx = 0 directly, by assuming that xi was of the form Ai cos(ωit + φi). In this method we consider each mode separately, by assuming that the amplitudes of all the other modes are zero (by some appropriate initial conditions). We can then substitute into the above equation to obtain

ωi2Mx + Kx = 0

i2M + K)x = 0.

Since non-trivial solutions require that det(ωi2M + K) = 0, we can find { ωi } as the roots of the characteristic equation.