By coincidence, right after completing last month's article Simultaneous Diagonalization and the Generalized Eigenvalue Problem, a long-time Maple user sent me a Maple worksheet whose purpose was to solve a three-degree-of-freedom spring-mass system using exactly the mathematics in the as-yet unreleased article. While editing the worksheet, it became clear that although the theory was easy enough to articulate, implementing an efficient solution of a real problem in Maple was something worth detailing.

Consequently, this month's article is a solution of the initial value problem I was sent. The purpose of the article is an articulation of the details of the calculations, showing how the theory developed in last month's article leads to a solution of a system of linear, second-order, constant-coefficient differential equations. In addition, we aim to show how the use of math-mode input simplifies the notation, and illuminates the connection between the underlying mathematics, and Maple's expression of it.

We take as given the following definitions of  and , the kinetic and potential energies of the system.

The Lagrangian  is then

so that the Euler-Lagrange equations for the problem are

where we have removed the first integrals that the EulerLagrange command generates. Except for the zeros on the right-hand sides, the individual equations are

If the following matrices are defined, the system can be written in the form , where u is a vector whose components are the unknowns .

The parameters appearing in the mass and stiffness equations are entered in set of equations, not hard-coded as assigned values. This keeps all the symbols as variables, and does not conflict the calculations with floats where they are not needed.

To obtain the numeric mass and stiffness matrices, evaluate the symbolic forms at the parameter values. By including one mass value as a float, all calculations involving these matrices will be executed in floating-point form. The exact representations of the eigenpairs for these matrices suffer badly from expression swell. This is to be expected, since the underlying math is the exact solution of a cubic characteristic equation.

The initial data for the problem is likewise entered as a set of equations, not hard-coded as assignments.

Although the system is linear, with constant coefficients, an exact solution suffers from expression swell because, in essence, an exact solution needs to access the exact roots of a cubic equation. A numeric solution is therefore sought. The specific ODEs for this problem are then

Invoke the dsolve command with the numeric option. The adroit use of sets allows the equations and the initial data to be coalesced by set union.

Obtain graphs of each dependent variable. Note the use of the odeplot command from the plots package.

The motion is highly oscillatory, so we expect the angular frequencies (which are the square roots of the generalized eigenvalues) to be relatively large.

This section sketches the algorithm for the simultaneous diagonalization of the mass and stiffness matrices. This algorithm requires both matrices to be symmetric, and for  to be positive definite, which in this case, it certainly is.

The differential equations are  . A matrix  exists, and with it, both M and K are simultaneously diagonalized. Simultaneous diagonalization of the matrices  and  means that  and , so that  and . Hence, the differential equation becomes

or

upon multiplying from the left by . Defining , or equivalently, , puts the differential equation into the uncoupled form , from which a solution is immediately available. The square roots of the diagonal elements in  are the angular frequencies for the normal modes of oscillation of this undamped vibrating system.

Alternatively, if solutions of the form  are sought, where v is a constant vector, the differential equation becomes , a generalized eigenvalue problem.

The matrix  is obtained as follows. Obtain , the Cholesky decomposition of , where  is a lower triangular factor. Form the matrix  and obtain its eigenpairs. Apparently, it is useful to sort the eigenpairs for physical reasons, but the mathematics underlying these calculations does not require this.

Form , a matrix whose columns are the normalized eigenvectors just obtained. (Maple returns normalized eigenvectors by default, so this step does not need to be implemented in the calculations below.) The required matrix  is given by .

The References section at the end of this document point to three different sources where this algorithm is derived.

Maple's Eigenvalues and Eigenvectors commands will solve the generalized eigenvalue problem. There's more work involved in sorting the eigenpairs than there is in finding them.

The generalized eigenvalues can be put into a vector , and their square roots, into a vector .

The following calculations implement the algorithm for finding  that simultaneously diagonalizes the mass and stiffness matrices. Begin by obtaining the Cholesky decomposition of .

Form the matrix  and obtain its eigenpairs, sorted. The eigenvectors are automatically normalized by Maple's Eigenvectors command. Form , a matrix of the sorted eigenvectors. Then define . These steps are summarized in Table 1.

The matrix  is

The number of digits in the elements of  is a function of the calculation that produced . The numerical routine that returned the Cholesky decomposition of  works in "hardware floats," that is, floating-point numbers implemented in the hardware of the computer. These are different from "software floats," the Maple floating-point numbers that can have an arbitrary number of digits. For a software float, the evalf command can be used to change the number of digits. Doing this with a hardware float requires a bit more care.

Table 2 contains some verifications of the preceding calculations.

To morph these matrices and vectors into an actual solution of the given IVP, write the general solution of the transformed (i.e., uncoupled) system , where here,  is a diagonal matrix whose diagonal elements are the components of the vector .

Note the use of the elementwise operator, the tilde. The construct  or  multiplies each element of the vector  by , and then computes the cosine or sine of that product. The vector of cosines is multiplied elementwise by a vector of constants ; the vector of sines, by a vector of constants .

Define vectors of initial values for the original initial value problem.

Recall that , where X is the solution of the uncoupled system, and u is the solution of the original coupled system. Since  is known, but  and  are given, the constants  and  in  are determined by the equations  and . The Equate and solve commands are instumental in setting up an solving these algebraic equations.

Substitute these values for the coefficients into Xg, the general solution of the uncoupled system.

This is the specific solution of the uncoupled system. It has to be transformed to the solution of the original coupled system. This final solution is , obtained below with all the digits of the hardware floats in , but displayed with just five digits.

Test that this solution satisfies the initial conditions:

Ordinarily, a graph would provide acceptable evidence that two solutions were equivalent. Given the highly oscillatory behavior of the solution of this initial value problem, we instead compare the values of the two solutions (one from dsolve, and one from uncoupling) at ten equispaced points in the interval .

Once again obtain the numeric solution provided by the dsolve command, but this time, set the output to be a list of procedures.

To see what the list of procedures is, call the solution at some given time.

Extract the procedures that give the values of the three dependent variables describing position.

At ten equispaced points in the interval , calculate the norm of the difference between the numeric solution and the vector solution obtained by the simultaneous diagonalization process detailed in this document.

To seven places, the solutions agree at these ten points.