# Having trouble working with two mass, three spring dynamic system

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.

• You should try DSolve rather than Solve Oct 10 '18 at 20:21

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.}}
• How would you get the eigenvalues of the system, as the example does next? I am having a hard time manipulating the equations of motion into the auxiliary equation format and taking its determinant to obtain the eigenfrequencies. Oct 11 '18 at 18:43
• @felimz: I edited the answer showing the two symbolic eigenvalues. Oct 11 '18 at 19:38
• Thanks for the expedient solution, but note that the Eigenvalues that are output from this solution are in terms of x1[t] and x2[t] and their derivatives. The eigenfrequencies should be a function of mass and stiffness only, and are time and position invariant: ω^2 = 1/2*(((k1+k2)/m1+(k2+k3)/m2)±(((k1+k2)/m1−(k2+k3)/m2)^2+4k2^2/(m1*m2))^(1/2)) Oct 11 '18 at 20:02
• @felimz: I edited the answer, and please see my answer for the derivation of the eigenvalues and the vectors. Oct 11 '18 at 20:59