1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You should try DSolve rather than Solve $\endgroup$ – Tugrul Temel Oct 10 '18 at 20:21
1
$\begingroup$

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.}}
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – felimz Oct 11 '18 at 18:43
  • $\begingroup$ @felimz: I edited the answer showing the two symbolic eigenvalues. $\endgroup$ – Tugrul Temel Oct 11 '18 at 19:38
  • $\begingroup$ 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)) $\endgroup$ – felimz Oct 11 '18 at 20:02
  • $\begingroup$ @felimz: I edited the answer, and please see my answer for the derivation of the eigenvalues and the vectors. $\endgroup$ – Tugrul Temel Oct 11 '18 at 20:59

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.