We propose a numerical algorithm for the computation of multi-marginal optimal transport (MMOT) problems involving general probability measures that are not necessarily discrete. By developing a relaxation scheme in which marginal constraints are replaced by finitely many linear constraints and by proving a specifically tailored duality result for this setting, we approximate the MMOT problem by a linear semi-infinite optimization problem. Moreover, we are able to recover a feasible and approximately optimal solution of the MMOT problem, and its sub-optimality can be controlled to be arbitrarily close to 0 under mild conditions. The developed approximation scheme leads to a numerical algorithm which can compute a feasible and approximately optimal solution of the MMOT problem with arbitrarily small sub-optimality. Besides the approximately optimal solution, the algorithm also computes upper and lower bounds for the optimal value of the MMOT problem. The difference between the computed bounds provides an explicit sub-optimality bound for the computed approximately optimal solution. We demonstrate the proposed algorithm in three numerical experiments involving an MMOT problem that stems from fluid dynamics, the Wasserstein barycenter problem, and a large-scale MMOT problem with 100 marginals. We observe that our algorithm is capable of computing high-quality solutions of these MMOT problems and the computed sub-optimality bounds are much less conservative than their theoretical upper bounds in all the experiments. Moreover, we compare our algorithm with existing regularization-based algorithms to showcase its advantages.