1
$\begingroup$

How can I write a program in mathematica that produce the solution of a linear chain of harmonic oscillators?

$$m\frac{d^2x_i}{dt^2}=k(x_{i-1}-2x_i+x_{i+1})$$

where $i=0...N+1$, $x_0$ and $x_{N+1}$ are constant (the walls).

It is simple to solve and find the eigenmodes using DSolve for cases like $i=2$ and get an exact solution.

I would like a program that receives the number of masses $n$ and can produce the trajectories of all masses, and all the eigenmodes.

$\endgroup$
1
  • 2
    $\begingroup$ What have you tried? $\endgroup$ May 21, 2018 at 16:49

2 Answers 2

4
$\begingroup$

This is one of a million ways to do that:

n = 20;
k0 = 1.;
m0 = RandomReal[{1, 10}, n];
k = SparseArray[{
 Band[{1, 1}] -> -2. k0, 
 Band[{1, 2}] -> k0, 
 Band[{2, 1}] -> k0
}, {n, n}, 0.];
m = DiagonalMatrix[SparseArray[m0]];
xx = Table[x[i][t], {i, 1, n}];

eq = Join[
   Thread[(m.D[xx, t, t])[[2 ;; -2]] == (k.xx)[[2 ;; -2]]],
   {x[1][t] == 1., x[n][t] == n},
   Thread[xx == Range[n] /. t -> 0],
   Thread[
    D[xx, t] == 0.25 Sin[Subdivide[0., 2. Pi, n - 1]] /. t -> 0]
   ];

sol = NDSolve[eq, xx, {t, 0, 100}];
Plot[Evaluate[xx /. sol[[1]]], {t, 0, 100}]

enter image description here

$\endgroup$
1
  • $\begingroup$ how can u retrieve the eigenmodes using your example $\endgroup$
    – jarhead
    May 21, 2018 at 18:20
4
$\begingroup$

It is also possible to use a vectorized NDSolve approach. First, I replicate @Henrik's answer, but add a random seed so that my results can be replicated:

SeedRandom[1]
n = 20;
k0 = 1.;
m0 = RandomReal[{1, 10}, n];
k = SparseArray[
    {
    Band[{1, 1}] -> -2. k0, 
    Band[{1, 2}] -> k0, 
    Band[{2, 1}] -> k0
    },
    {n, n},
    0.
];
m = DiagonalMatrix[SparseArray[m0]];
xx = Table[x[i][t], {i, 1, n}];

eq = Join[
    Thread[(m.D[xx, t, t])[[2 ;; -2]] == (k.xx)[[2 ;; -2]]],
    {x[1][t] == 1., x[n][t] == n},
    Thread[xx == Range[n] /. t -> 0],
    Thread[D[xx, t] == 0.25 Sin[Subdivide[0., 2. Pi, n - 1]] /. t -> 0]
];

sol = NDSolve[eq, xx, {t, 0, 100}];
Plot[Evaluate[xx /. sol[[1]]], {t, 0, 100}]

enter image description here

The vectorized approach:

sol2 = NDSolveValue[
    {
    X''[t] == ArrayPad[IdentityMatrix[n-2], 1] . Inverse[m] . k . X[t],
    X[0] == Range[n],
    X'[0] == 0.25 Sin[Subdivide[0.,2. Pi,n-1]]
    },
    X,
    {t, 0, 100}
];

I use a truncated identity matrix to set the variation of the first and last elements of the vector to zero, so that they are constant. And a visualization:

Plot[
    Evaluate @ Table[Indexed[sol2[t], i], {i,20}],
    {t, 0, 100}
]

enter image description here

$\endgroup$
1
  • $\begingroup$ this is great thanks!, do u know how can I plot the modes, and which of them are "used" as a function of time? $\endgroup$
    – jarhead
    May 22, 2018 at 13:36

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