Classroom Tips and Techniques: An Undamped coupled Oscillator
Robert J. Lopez
Emeritus Professor of Mathematics and Maple Fellow
Maplesoft

Introduction


By coincidence, right after completing last month's article Simultaneous Diagonalization and the Generalized Eigenvalue Problem, a longtime Maple user sent me a Maple worksheet whose purpose was to solve a threedegreeoffreedom springmass system using exactly the mathematics in the asyet 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, secondorder, constantcoefficient differential equations. In addition, we aim to show how the use of mathmode input simplifies the notation, and illuminates the connection between the underlying mathematics, and Maple's expression of it.


Data Definitions


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


The Lagrangian is then
>


so that the EulerLagrange equations for the problem are
>


where we have removed the first integrals that the EulerLagrange command generates. Except for the zeros on the righthand 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 hardcoded 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 floatingpoint 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 hardcoded as assignments.
>




dsolve Solution


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.


Structure


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.


Generalized Eigenvalue Problem and Its Eigenpairs


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 .
>




Simultaneous Diagonalization


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.
•

Form the matrix .


>




>



•

Sort the three lists by ordering the eigenvalues from smallest to largest.


>



•

Construct the matrix whose columns are the eigenvectors, now both unit and sorted.


>



•

Form the matrix .


>



•

Form the vector whose elements are the sorted eigenvalues. These are the squares of the modal angular frequencies.


>



•

Form the vector whose elements are the modal angular frequencies.


>



Table 1 Obtaining the matrix



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, floatingpoint numbers implemented in the hardware of the computer. These are different from "software floats," the Maple floatingpoint 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.
•

Show that .


>



•

Show that .


>



•

Show that ; i.e., that is an orthogonal matrix.


>



Table 2 Verifications: , , and



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:


Comparing Solutions


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.


References


1.

Numerical Analysis: A Practical Approach, 3rd ed., Melvin J. Marin and Robert J. Lopez, Wadsworth Publishing Company, 1991.

3.

Applied Linear Algebra, 3rd ed., Ben Noble and James W. Daniel, PrenticeHall Publishing Company, 1988.


Legal Notice: © Maplesoft, a division of Waterloo Maple Inc. 2012. Maplesoft and Maple are trademarks of Waterloo Maple Inc. This application may contain errors and Maplesoft is not liable for any damages resulting from the use of this material. This application is intended for noncommercial, nonprofit use only. Contact Maplesoft for permission if you wish to use this application in forprofit activities.
