Application Center - Maplesoft

App Preview:

Classroom Tips and Techniques: An Undamped Coupled Oscillator

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


Classroom Tips and Techniques: An Undamped coupled Oscillator


Robert J. Lopez

Emeritus Professor of Mathematics and Maple Fellow





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.








interface(typesetting = extended)



Data Definitions


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


KE := (1/2)*(sum(m[j]*(diff(x[j](t), t))^2, j = 1 .. 3)); PE := 1/2*(sum(k[j]*x[j](t)^2, j = 1 .. 3)+k[4]*(x[2](t)-x[1](t))^2+k[5]*(x[3](t)-x[2](t))^2+k[6]*(x[3](t)-x[1](t))^2)


The Lagrangian L = KE-PE is then


Lt := KE-PE


so that the Euler-Lagrange equations for the problem are


EL := convert(remove(has, EulerLagrange(Lt, t, [x[1](t), x[2](t), x[3](t)]), K), list)


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


for j to 3 do E || j := -select(has, EL, m[j])[] end do

k[1]*x[1](t)-k[4]*(x[2](t)-x[1](t))-k[6]*(x[3](t)-x[1](t))+m[1]*(diff(diff(x[1](t), t), t))

k[2]*x[2](t)+k[4]*(x[2](t)-x[1](t))-k[5]*(x[3](t)-x[2](t))+m[2]*(diff(diff(x[2](t), t), t))

k[3]*x[3](t)+k[5]*(x[3](t)-x[2](t))+k[6]*(x[3](t)-x[1](t))+m[3]*(diff(diff(x[3](t), t), t))


If the following matrices are defined, the system can be written in the form M*`#mover(mi("u",fontstyle = "normal",fontweight = "bold"),mrow(mo(".."),mo("⁢")))`+K*u = 0, where u is a vector whose components are the unknowns x[j](t).


SymbolicK := GenerateMatrix([E1, E2, E3], [seq(x[j](t), j = 1 .. 3)])[1]; SymbolicM := DiagonalMatrix([m[1], m[2], m[3]])

SymbolicK := Matrix(3, 3, {(1, 1) = k[1]+k[6]+k[4], (1, 2) = -k[4], (1, 3) = -k[6], (2, 1) = -k[4], (2, 2) = k[2]+k[5]+k[4], (2, 3) = -k[5], (3, 1) = -k[6], (3, 2) = -k[5], (3, 3) = k[6]+k[3]+k[5]})

Matrix(%id = 106753616)


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.


params := {k[1] = 1000, k[2] = 1000, k[3] = 1000, k[4] = 1000, k[5] = 1000, k[6] = 1000, m[1] = 1/2, m[2] = 1., m[3] = 3/2}


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.


M := eval(SymbolicM, params)
K := eval(SymbolicK, params)

M := Matrix(3, 3, {(1, 1) = 1/2, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1., (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 3/2})

Matrix(%id = 106757648)


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


inits := {x[1](0) = .5, x[2](0) = -0.5e-1, x[3](0) = 0, (D(x[1]))(0) = 0, (D(x[2]))(0) = 0, (D(x[3]))(0) = 0}



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


DEs := eval({E1, E2, E3}, params)

{3000*x[1](t)-1000*x[2](t)-1000*x[3](t)+(1/2)*(diff(diff(x[1](t), t), t)), 3000*x[2](t)-1000*x[1](t)-1000*x[3](t)+1.*(diff(diff(x[2](t), t), t)), 3000*x[3](t)-1000*x[2](t)-1000*x[1](t)+(3/2)*(diff(diff(x[3](t), t), t))}


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.


Nsol := dsolve(`union`(DEs, inits), numeric)


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


odeplot(Nsol, [t, x[1](t)], t = 0 .. 1)

odeplot(Nsol, [t, x[2](t)], t = 0 .. 1)

odeplot(Nsol, [t, x[3](t)], t = 0 .. 1)


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 M to be positive definite, which in this case, it certainly is.


The differential equations are  M*`#mover(mi("u",fontstyle = "normal",fontweight = "bold"),mrow(mo(".."),mo("⁢")))`+K*u = 0. A matrix S exists, and with it, both M and K are simultaneously diagonalized. Simultaneous diagonalization of the matrices K and M means that S^T*M*S = I and S^T*K*S = Lambda, so that M = 1/(S^T*S) and K = Lambda/(S^T*S). Hence, the differential equation becomes


`#mover(mi("u",fontstyle = "normal",fontweight = "bold"),mrow(mo(".."),mo("⁢")))`/(S^T*S)+Lambda*u/(S^T*S) = 0 


`#mover(mi("u",fontstyle = "normal",fontweight = "bold"),mrow(mo(".."),mo("⁢")))`/S+Lambda*u/S = 0 


upon multiplying from the left by S^T. Defining X = u/S, or equivalently, u = S*X, puts the differential equation into the uncoupled form diff(X(t), t, t)+Lambda*X(t) = 0, from which a solution is immediately available. The square roots of the diagonal elements in Lambda are the angular frequencies for the normal modes of oscillation of this undamped vibrating system.


Alternatively, if solutions of the form u = e^(omega*t)*v are sought, where v is a constant vector, the differential equation becomes omega^2*M*v+K*v = 0, a generalized eigenvalue problem.


The matrix S is obtained as follows. Obtain M = LL^T, the Cholesky decomposition of M, where L is a lower triangular factor. Form the matrix K(1/L)^T/L 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 H, 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 S is given by S = (1/L)^T*H.


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.


temp := map(simplify, [Eigenvectors(K, M, output = list)])
Eigenpairs := sort(temp[], proc (a, b) options operator, arrow; is(a[1] < b[1]) end proc)


[950.2505048, 1, {Vector(3, {(1) = -.4873561197, (2) = -.5632820932, (3) = -.6672310665})}]

[3351.809967, 1, {Vector(3, {(1) = .2295974482, (2) = .8232250697, (3) = -.5192162328})}]

[6697.939528, 1, {Vector[column](%id = 112165364)}]


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


Lambda := `<,>`(seq(Eigenpairs[j][1], j = 1 .. 3))
lambda := map(sqrt, Lambda)

Lambda := Vector(3, {(1) = 950.2505048, (2) = 3351.809967, (3) = 6697.939528})

Vector[column](%id = 112166836)



Simultaneous Diagonalization


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


L := LUDecomposition(M, method = Cholesky)

Typesetting:-mrow(Typesetting:-mi("L", italic = "true", mathvariant = "italic"), Typesetting:-mo("&Assign;", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.2777778em", rspace = "0.2777778em"), Typesetting:-mverbatim("-I'RTABLEG6"6%"*/Lk$H-I'MATRIXG6"6#7%7%$"0[l="y1rq!#:$""!""!$""!""!7%$""!""!$"""""!$""!""!7%$""!""!$""!""!$"0f"Rr[uC7!#9I'MatrixG6$%*protectedGI(_syslibG6""))


Form the matrix K(1/L)^T/L and obtain its eigenpairs, sorted. The eigenvectors are automatically normalized by Maple's Eigenvectors command. Form H, a matrix of the sorted eigenvectors. Then define S = (1/L)^T*H. These steps are summarized in Table 1.



Form the matrix K(1/L)^T/L.

Temp1 := 1/L.K.(1/L)^%T


Obtain the eigenpairs of this matrix, each in the form [a, b, c], where a is the eigenvalue, b is its algebraic multiplicity, and c is a set of any corresponding eigenvectors.

Temp2 := simplify([Eigenvectors(Temp1, output = list)])


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

EP := sort(Temp2[], proc (a, b) options operator, arrow; is(a[1] < b[1]) end proc)


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

H := Matrix([seq(EP[k][3][], k = 1 .. 3)])


Form the matrix S = (1/L)^T*H.

S := (1/L)^%T.H


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


Lambda := `<,>`(seq(EP[k][1], k = 1 .. 3))

Vector[column](%id = 128405000)


Form the vector lambda whose elements are the modal angular frequencies.


lambda := map(sqrt, Lambda)

Vector[column](%id = 84079216)

Table 1   Obtaining the matrix S


The matrix S is



Matrix([[-.463866737142020, -.218078018204323, 1.31805524469171], [-.536133262900000, -.781921981900000, -.318055244700000], [-.635072147896599, .493165964854146, -.141906183042453]])


The number of digits in the elements of S is a function of the calculation that produced L. The numerical routine that returned the Cholesky decomposition of M 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.


map(evalf[10], S)

Matrix(%id = 127654960)


Table 2 contains some verifications of the preceding calculations.



Show that S^T*M*S = I.

map(`@`(simplify, fnormal), S^%T.M.S, 9)

Matrix(%id = 108118508)


Show that S^T*K*S = Lambda.

map(`@`(simplify, fnormal), S^%T.K.S, 8)

Matrix(%id = 110401892)


Show that H^T*H = I; i.e., that H is an orthogonal matrix.

map(`@`(simplify, fnormal), H^%T.H, 9)

Matrix(%id = 110402276)

Table 2   Verifications: S^T*M*S = I, S^T*K*S = Lambda, and H^T*H = I


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 `#mover(mi("X",fontstyle = "normal",fontweight = "bold"),mrow(mo("&period;&period;"),mo("&InvisibleTimes;")))`+Lambda*X = 0, where here, Lambda is a diagonal matrix whose diagonal elements are the components of the vector Lambda.


Xg := `~`[`*`](`<,>`(seq(a[j], j = 1 .. 3)), `~`[cos](lambda*t))+`~`[`*`](`<,>`(seq(b[j], j = 1 .. 3)), `~`[sin](lambda*t))

Vector[column](%id = 110403236)


Note the use of the elementwise operator, the tilde. The construct `~`[cos](lambda*t) or `~`[sin](lambda*t) multiplies each element of the vector lambda by t, and then computes the cosine or sine of that product. The vector of cosines is multiplied elementwise by a vector of constants a[j]; the vector of sines, by a vector of constants b[j].


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


x0 := `<,>`(.5, -0.5e-1, 0)
xdot0 := `<,>`(0, 0, 0)NULL

x0 := Vector(3, {(1) = .5, (2) = -0.5e-1, (3) = 0})

Vector[column](%id = 99240132)


Recall that X = u/S, where X is the solution of the uncoupled system, and u is the solution of the original coupled system. Since X(t) is known, but u(0) and (D(u))(0) are given, the constants a[j] and b[j] in X(t) are determined by the equations X(0) = u(0)/S and (D(X))(0) = (D(u))(0)/S. The Equate and solve commands are instumental in setting up an solving these algebraic equations.


A := solve(Equate(eval(Xg, t = 0), 1/S.x0))B := solve(Equate(eval(map(diff, Xg, t), t = 0), 1/S.xdot0))NULL

{a[1] = -0.8916002112e-1, a[2] = -0.1542340544e-1, a[3] = .3454165734}

{b[1] = 0., b[2] = 0., b[3] = 0.}


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


X := eval(Xg, `union`(A, B))

Vector[column](%id = 110408868)


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 u = S*X, obtained below with all the digits of the hardware floats in S, but displayed with just five digits.


u := S.X; map(evalf[5], u)

Vector[column](%id = 131019932)


Test that this solution satisfies the initial conditions:


map(`@`(simplify, fnormal), eval(u, t = 0), 9)
map(`@`(simplify, fnormal), eval(map(diff, u, t), t = 0))

Vector(3, {(1) = .5000000000, (2) = -0.5000000000e-1, (3) = 0.})

Vector[column](%id = 105421420)



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 .1 <= t and t <= 1.


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


Nsol := dsolve(`union`(DEs, inits), numeric, output = listprocedure)


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



[t(1) = 1., (x[1](t))(1) = HFloat(0.48462071608167273), (diff(x[1](t), t))(1) = HFloat(-5.399770918272315), (x[2](t))(1) = HFloat(-0.0660550136699385), (diff(x[2](t), t))(1) = HFloat(1.5672045432948187), (x[3](t))(1) = HFloat(-0.003032171633109441), (diff(x[3](t), t))(1) = HFloat(2.03737405500597)]


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


Y1 := rhs(Nsol[2])
Y2 := rhs(Nsol[4])

Y3 := rhs(Nsol[6])


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


seq(LinearAlgebra[Norm](map(`@`(simplify, fnormal), evalf(eval(`<,>`(Y1(t), Y2(t), Y3(t))-u, t = (1/10)*j*Pi)), 7)), j = 1 .. 10)

0., 0., 0., 0., 0., 0., 0., 0., 0., 0.


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







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



Advanced Engineering Mathematics, Robert J. Lopez, Addison-Wesley Publishing Company, 2001. Now available from Maplesoft as the ebook Advanced Engineering Mathematics with Maple.



Applied Linear Algebra, 3rd ed., Ben Noble and James W. Daniel, Prentice-Hall 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 non-commercial, non-profit use only. Contact Maplesoft for permission if you wish to use this application in for-profit activities.