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