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

Question

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}]

Answer 1

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}}} *)