Having trouble working with two mass, three spring dynamic system

Question

I am following a dynamic analysis example https://www.math24.net/mass-spring-system/ and trying to implement it in Mathematica.

However, I am having trouble even getting simple properties like the determinant and the eigenvalues of the set of differential equations prior to attempting to solve it.

I can obtain the set of differential equations with this code:

ClearAll["Global`*"]
L = 
  1/2 (-k3 x2[t]^2 + m2 x2'[t]^2 - k2 (-x1[t] + x2[t])^2 - k1 x1[t]^2 + m1 x1'[t]^2);

Lpdx1 = D[L, x1'[t]];
Lpx1 = D[L, x1[t]];
Lpdx2 = D[L, x2'[t]];
Lpx2 = D[L, x2[t]];

DiffEqRules = 
  Apart[
    FullSimplify[
      {Solve[D[Lpdx1, t] == D[Lpx1], x1''[t]], 
       Solve[D[Lpdx2, t] == D[Lpx2], x2''[t]]}]];

DiffEqs = DiffEqRules[[All, 1]][[All, 1]] /. Rule -> Equal

{(x1'')[t] == -(((k1 + k2) x1[t])/m1) + (k2 x2[t])/m1, 
 (x2'')[t] == (k2 x1[t])/m2 - ((k2 + k3) x2[t])/m2}

Anything that helps me get any further in this example will be appreciated.

Answer 1

The following code may be of help. Given the parameter values:

k1 = 3;
k2 = 1;
k3 = 2;
m1 = 4;
m2 = 5;

DSolve[{D[Lpdx1, t] == D[Lpx1], D[Lpdx2, t] == D[Lpx2], x1[0] == 0, 
  x1'[0] == 1, x2[0] == 0, x2'[0] == 2}, {x1, x2}, t]

(* output *)
{{x1 -> Function[{t}, Sqrt[2] Sin[t/Sqrt[2]]], 
  x2 -> Function[{t}, 2 Sqrt[2] Sin[t/Sqrt[2]]]}}

If you run the following code:

ClearAll[a1, a2, b1, b2, k1, k2, k3, m1, m2];
L = 1/2 (-k3 x2[t]^2 + m2 x2'[t]^2 - k2 (-x1[t] + x2[t])^2 - 
     k1 x1[t]^2 + m1 Derivative[1][x1][t]^2);

Lpdx1 = D[L, x1'[t]];
Lpx1 = D[L, x1[t]];
Lpdx2 = D[L, x2'[t]];
Lpx2 = D[L, x2[t]];

DSolve[{D[Lpdx1, t] == D[Lpx1], D[Lpdx2, t] == D[Lpx2], x1[0] == a1, 
  x1'[0] == a2, x2[0] == b1, x2'[0] == b2}, {x1, x2}, t]

You get a symbolic solution (which is very long one, therefore I do not give it here).

To get the eigenvalues, try this:

Eigenvalues[{{Lpx1, Lpdx1}, {Lpx2, Lpdx2}}]

which yields

{1/2 (-k1 x1[t] - k2 x1[t] + k2 x2[t] + 
m2 Derivative[1][x2][
  t] - \[Sqrt]((k1 x1[t] + k2 x1[t] - k2 x2[t] - 
     m2 Derivative[1][x2][t])^2 - 
   4 (-k2 m1 x1[t] Derivative[1][x1][t] + 
      k2 m1 x2[t] Derivative[1][x1][t] + 
      k3 m1 x2[t] Derivative[1][x1][t] - 
      k1 m2 x1[t] Derivative[1][x2][t] - 
      k2 m2 x1[t] Derivative[1][x2][t] + 
      k2 m2 x2[t] Derivative[1][x2][t]))), 
1/2 (-k1 x1[t] - k2 x1[t] + k2 x2[t] + 
m2 Derivative[1][x2][
  t] + \[Sqrt]((k1 x1[t] + k2 x1[t] - k2 x2[t] - 
     m2 Derivative[1][x2][t])^2 - 
   4 (-k2 m1 x1[t] Derivative[1][x1][t] + 
      k2 m1 x2[t] Derivative[1][x1][t] + 
      k3 m1 x2[t] Derivative[1][x1][t] - 
      k1 m2 x1[t] Derivative[1][x2][t] - 
      k2 m2 x1[t] Derivative[1][x2][t] + 
      k2 m2 x2[t] Derivative[1][x2][t])))}

Maybe this:

eqs = {
        x1''[t] == -(((k1 + k2) x1[t])/m1) + (k2 x2[t])/m1,
        x2''[t] == (k2 x1[t])/m2 - ((k2 + k3) x2[t])/m2
      };

(* Jacobian for the right-hand side of the "eqs" *)
jac = D[eqs[[All, 2]], {{x1[t], x2[t]}}]

(* Find the equilibrium points *)
eqPoints = Solve[eqs /. {x1''[t] -> 0, x2''[t] -> 0}, {x1[t], x2[t]}]

(* Find the eigenvalues and eigenvectors of the Jacobian *)
{vals, vecs} = Eigensystem[N[jac /. eqPoints[[1]]]]

Then the two eigenvalues are:

{(0.5 (-1. k2 m1 - 1. k3 m1 - 1. k1 m2 - 1. k2 m2 + 
    1. Sqrt[-4. (1. k1 k2 + 1. k1 k3 + 1. k2 k3) m1 m2 + (1. k2 m1 + 
        1. k3 m1 + 1. k1 m2 + 1. k2 m2)^2]))/(
 m1 m2), -((
  0.5 (1. k2 m1 + 1. k3 m1 + 1. k1 m2 + 1. k2 m2 + 
     1. Sqrt[-4. (1. k1 k2 + 1. k1 k3 + 1. k2 k3) m1 m2 + (1. k2 m1 + 
         1. k3 m1 + 1. k1 m2 + 1. k2 m2)^2]))/(m1 m2))

Then the two eigenvectors are:

{{(0.5 (1. k2 m1 + 1. k3 m1 - 1. k1 m2 - 1. k2 m2 + 
     1. Sqrt[-4. (1. k1 k2 + 1. k1 k3 + 1. k2 k3) m1 m2 + (1. k2 m1 + 
         1. k3 m1 + 1. k1 m2 + 1. k2 m2)^2]))/(k2 m1), 1.}, {(
  0.5 (1. k2 m1 + 1. k3 m1 - 1. k1 m2 - 1. k2 m2 - 
     1. Sqrt[-4. (1. k1 k2 + 1. k1 k3 + 1. k2 k3) m1 m2 + (1. k2 m1 + 
         1. k3 m1 + 1. k1 m2 + 1. k2 m2)^2]))/(k2 m1), 1.}}