Can I use DEigensystem to solve coupled ordinary differential equations such as a line of coupled harmonic oscillators?

I am working on systems of multiple coupled harmonic oscillators, where we might in the simplest case have a single line of masses each coupled with springs to its neighbor masses, with something like the following system of 4 ODEs:

m1 x1''[t] =  k(x2[t] - x1[t]) + k( 0    - x1[t])
m2 x2''[t] =  k(x3[t] - x2[t]) + k(x1[t] - x2[t])
m3 x3''[t] =  k(x4[t] - x3[t]) + k(x2[t] - x3[t])
m4 x4''[t] =  k( 0    - x4[t]) + k(x3[t] - x4[t])


With the two end springs attached to fixed positions, the whole line of oscillators acts like a discrete version of a stretched string. It exhibits normal modes, which are the finite-dimensional eigenvectors, oscillating at frequencies which are the corresponding eigenvalues. The eigenvalue appears as omega in the exponent of the solution vector, for example,

{x1[t_] := a1 Exp[I omega t], x2[t_]:= a2 Exp[I omega t]...} and so on.


Although DSolve works for these, I want to experiment with DEigensystem. Unfortunately almost all examples for DEigensystem in the documentation are for partial differential equations. I have looked for a good example of DEigensystem with coupled ODEs, but I can't seem to find one, and my attempts to call the function have failed with the function and its arguments returned unchanged. Can DEigensystem be used for that kind of system of ODEs? If so, does anyone have an example?

Here is a working example using DSolve:

(* Begin working code example *)

Needs["VariationalMethods"]

\[Omega]Substitutions = First[Solve[{Sqrt[k]/Sqrt[m] == \[Omega]1, Sqrt[k + 2*kprime]/Sqrt[m] == \[Omega]2}, {k, kprime}]]

ke = (1/2)*m*Derivative[1][x1][t]^2 +(1/2)*m*Derivative[1][x2][t]^2;
pe = (1/2)*k*x1[t]^2 + (1/2)*kprime*(x2[t] - x1[t])^2 + (1/2)*k*x2[t]^2;
lagrangian = ke - pe /. \[Omega]Substitutions

eulerLagrEquations = EulerEquations[lagrangian, {x1[t], x2[t]}, t]

equationsAndInitialConditions = Join[{eulerLagrEquations, x1[0] == x10, x2[0] == x20, Derivative[1][x1][0] == v10, Derivative[1][x2][0] == v20}];

solutions = FullSimplify[DSolve[equationsAndInitialConditions, {x1[t], x2[t]}, t]]

numericSolutions = N[solutions] /. {\[Omega]1 -> 0.1, \[Omega]2 -> 0.142857, m -> 1, v10 -> 0, v20 -> 0, x10 -> 0, x20 -> 0.05};
Plot[{x1[t] /. numericSolutions, x2[t] + 0.1 /. numericSolutions}, {t, 0, 200}, PlotRange -> {-0.15, 0.15}]

• Please provide a small complete system, perhaps three oscillators. Also, where is the eigenvalue in your equation? Commented Jun 6, 2018 at 3:50
• Note that use of EulerEquations requires the VariationalMethods package. You may wish to edit your code to include Needs["VariationalMethods"]. Commented Jun 6, 2018 at 13:16
• I don't understand what you are doing with the substitution in terms of $x_i$. If you do that with the same $\omega$ then you get a system of non differential equations. Whereas you seem to be doing something else in the example you give, which doesn't have an eigenvalue. Commented Jun 6, 2018 at 13:21

While the coupled oscillator system involves differential equations, and can be solved using eigenvector methods, it is not a "differential eigenvalue problem" in the usual sense. Specifically, the problems that DEigensystem can solve are finding the values of $\lambda$ such that $$\mathcal{L} u = \lambda u \tag{1}$$ has solutions subject to some boundary conditions. Here, $\mathcal{L}$ is a linear differential operator, and $u$ is some vector of functions on a given domain.

In your case, you have a system of the form $$M \frac{d^2}{dt^2} \begin{bmatrix} x_1\\x_2 \end{bmatrix} = K \begin{bmatrix} x_1\\x_2 \end{bmatrix},$$ where $M$ and $K$ are matrices. Since there is no arbitrary parameter $\lambda$ for which your equation takes the form (1) (unless $K$ is proportional to the identity matrix), it's not of the correct form to be handled by DEigensystem.

I suppose that it's possible that you could view this problem as a "generalized differential eigenvalue problem", in the same way we can view certain matrix problems as generalized eigenvalue problems. Mathematica can solve generalized eigenvalue problems involving matrices; but it does not appear to have this capacity for differential problems. Perhaps in a later version this will be addressed.

Honestly, if you want to solve this problem as an eigenvalue problem, it's probably easier to extract the mass matrix and potential matrix from ke and pe directly:

pecoeffmat = Normal[CoefficientArrays[pe, {x1[t], x2[t]}][[3]]];
potmat = 1/2 (pecoeffmat + Transpose[pecoeffmat])
kecoeffmat = Normal[CoefficientArrays[ke, {x1'[t], x2'[t]}][[3]]];
kinmat = 1/2 (kecoeffmat + Transpose[kecoeffmat])
Eigensystem[{potmat, kinmat}]

(* {{(k + kprime)/2, -(kprime/2)}, {-(kprime/2), (k + kprime)/2}} *)
(* {{m/2, 0}, {0, m/2}} *)
(* {{k/m, (k + 2 kprime)/m}, {{1, 1}, {-1, 1}}} *)

• Mathematica can and did solve the problem in the working example I provided. DSolve found the complete solution set, and I am not unhappy with the results. I just wanted to see if I could get DEigensystem to work with this, and consider whether it might be the better tool to use. Commented Jun 6, 2018 at 15:03
• The code you suggest is interesting and helpful. Thank you. Commented Jun 6, 2018 at 16:07
• But now are you sure DEigensystem can't be used for this somehow? The documentation seems to suggest something like it: "DEigensystem[eqns,u,t,{x,y,[Ellipsis]}[Element][CapitalOmega],n] gives the eigenvalues and eigenfunctions for solutions u of the time-dependent differential equations eqns." Commented Jun 6, 2018 at 18:16
• @RalphDratman: It looks like in that context, DEigensystem is looking for PDEs involving time on a $D$-dimensional domain Ω, for a total of $D + 1$ independent variables. Your system has only 1 independent variable ($t$), so viewed that way, your domain Ω consists of a zero-dimensional set — a single point. The $x_i(t)$ are then functions from $\mathbb{R} \times \Omega \simeq \mathbb{R} \to \mathbb{R}$. Commented Jun 6, 2018 at 18:38
• I suppose it might be possible to define Ω this way and invoke DEigensystem, but I wouldn't be surprised if it throws an error. (I'd try it out myself, but I'm running a long calculation on my kernel at the moment...) Commented Jun 6, 2018 at 18:39