What's the difference between principal component analysis and multidimensional scaling? How are PCA and classical MDS different? How about MDS versus nonmetric MDS? Is there a time when you would prefer one over the other? How do the interpretations differ?
 A: Two types of metric MDS
The task of metric multidimensional scaling (MDS) can be abstractly formulated as follows: given a $n\times n$ matrix  $\mathbf D$ of pairwise distances between $n$ points, find a low-dimensional embedding of data points in $\mathbb R^k$ such that Euclidean distances between them approximate the given distances: $$\|\mathbf x_i - \mathbf x_j\|\approx D_{ij}.$$
If "approximate" here is understood in a the usual sense of reconstruction error, i.e. if the goal is to minimize the cost function called "stress": $$\text{Stress} \sim \Big\|\mathbf D - \|\mathbf x_i - \mathbf x_j\|\Big\|^2,$$ then the solution is not equivalent to PCA. The solution is not given by any closed formula, and must be computed by a dedicated iterative algorithm.
"Classical MDS", also known as "Torgerson MDS", replaces this cost function by a related but not equivalent one, called "strain": $$\text{Strain} \sim \Big\|\mathbf K_c - \langle\mathbf x_i, \mathbf x_j\rangle\Big\|^2,$$ that seeks to minimize reconstruction error of centered scalar products instead of distances. It turns out that $\mathbf K_c$ can be computed from $\mathbf D$ (if $\mathbf D$ are Euclidean distances) and that minimizing reconstruction error of $\mathbf K_c$ is exactly what PCA does, as shown in the next section.
Classical (Torgerson) MDS on Euclidean distances is equivalent to PCA
Let the data be collected in matrix $\mathbf X$ of $n \times k$ size with observations in rows and features in columns. Let $\mathbf X_c$ be the centered matrix with subtracted column means.
Then PCA amounts to doing singular value decomposition $\mathbf X_c = \mathbf {USV^\top}$, with columns of $\mathbf{US}$ being principal components. A common way to obtain them is via an eigendecomposition of the covariance matrix $\frac{1}{n}\mathbf X_c^\top \mathbf X^\vphantom{\top}_c$, but another possible way is to perform an eigendecomposition of the Gram matrix $\mathbf K_c = \mathbf X^\vphantom{\top}_c \mathbf X^\top_c=\mathbf U \mathbf S^2 \mathbf U^\top$: principal components are its eigenvectors scaled by the square roots of the respective eigenvalues.
It is easy to see that $\mathbf X_c = (\mathbf I - \frac{1}{n}\mathbf 1_n)\mathbf X$, where $\mathbf 1_n$ is a $n \times n$ matrix of ones. From this we immediately get that $$\mathbf K_c = \left(\mathbf I - \frac{\mathbf 1_n}{n}\right)\mathbf K\left(\mathbf I - \frac{\mathbf 1_n}{n}\right) = \mathbf K - \frac{\mathbf 1_n}{n} \mathbf K - \mathbf K \frac{\mathbf 1_n}{n} + \frac{\mathbf 1_n}{n} \mathbf K \frac{\mathbf 1_n}{n},$$ where $\mathbf K = \mathbf X \mathbf X^\top$ is a Gram matrix of uncentered data. This is useful: if we have the Gram matrix of uncentered data we can center it directly, without getting back to $\mathbf X$ itself. This operation is sometimes called double-centering: notice that it amounts to subtracting row means and column means from $\mathbf K$ (and adding back the global mean that gets subtracted twice), so that both row means and column  means of   $\mathbf K_c$ are equal to zero.
Now consider a $n \times n$ matrix $\mathbf D$ of pairwise Euclidean distances with $D_{ij} = \|\mathbf x_i - \mathbf x_j\|$. Can this matrix be converted into $\mathbf K_c$ in order to perform PCA? Turns out that the answer is yes.
Indeed, by the law of cosines we see that 
\begin{align}
D_{ij}^2 = \|\mathbf x_i - \mathbf x_j\|^2 &= \|\mathbf x_i - \bar{\mathbf x}\|^2 + \|\mathbf x_j - \bar{\mathbf x}\|^2 - 2\langle\mathbf x_i - \bar{\mathbf x}, \mathbf x_j - \bar{\mathbf x} \rangle \\ &= \|\mathbf x_i - \bar{\mathbf x}\|^2 + \|\mathbf x_j - \bar{\mathbf x}\|^2 - 2[K_c]_{ij}.
\end{align} So $-\mathbf D^2/2$ differs from $\mathbf K_c$ only by some row and column constants (here $\mathbf D^2$ means element-wise square!). Meaning that if we double-center it, we will get $\mathbf K_c$: $$\mathbf K_c = -\left(\mathbf I - \frac{\mathbf 1_n}{n}\right)\frac{\mathbf D^2}{2}\left(\mathbf I - \frac{\mathbf 1_n}{n}\right).$$
Which means that starting from the matrix of pairwise Euclidean distances $\mathbf D$ we can perform PCA and get principal components. This is exactly what classical (Torgerson) MDS does: $\mathbf D \mapsto \mathbf K_c \mapsto \mathbf{US}$, so its outcome is equivalent to PCA.
Of course, if any other distance measure is chosen instead of $\|\mathbf x_i - \mathbf x_j\|$, then classical MDS will result in something else.
Reference: The Elements of Statistical Learning, section 18.5.2.
A: Uhm... quite different. In PCA, you are given the multivariate continuous data (a multivariate vector for each subject), and you are trying to figure out if you don't need that many dimensions to conceptualize them. In (metric) MDS, you are given the matrix of distances between the objects, and you are trying to figure out what the locations of these objects in space are (and whether you need a 1D, 2D, 3D, etc. space). In non-metric MDS, you only know that objects 1 and 2 are more distant than objects 2 and 3, so you try to quantify that, on top of finding the dimensions and locations.
With a notable stretch of imagination, you can say that a common goal of PCA and MDS is to visualize objects in 2D or 3D. But given how different the inputs are, these methods won't be discussed as even distantly related in any multivariate textbook. I would guess that you can convert the data usable for PCA into data usable for MDS (say, by computing Mahalanobis distances between objects, using the sample covariance matrix), but that would immediately result in a loss of information: MDS is only defined up to location and rotation, and the latter two can be done more informatively with PCA.
If I were to briefly show someone the results of non-metric MDS and wanted to give them a rough idea of what it does without going into detail, I could say:

Given the measures of similarity or dissimilarity that we have, we are trying to map our objects/subjects in such a way that the 'cities' they make up have distances between them that are as close to these similarity measures as we can make them. We could only map them perfectly in $n$-dimensional space, though, so I am representing the two most informative dimensions here -- kinda like what you would do in PCA if you showed a picture with the two leading principal components.

A: Comparison: "Metric MDS gives the SAME result as PCA"- procedurally- when we look at the way SVD is used to obtain the optimum. But, the preserved high-dimensional criteria is different. PCA uses a centered covariance matrix while MDS uses a gram matrix obtained by double-centering distance matrices.
Will put the difference mathematically: PCA can be viewed as maximizing $Tr(X^T(I-\frac{1}{n}ee^T)X)$ over $X$ under constraints that $X$ is orthogonal, thereby giving axes/principal components. In multidimensional scaling a gram matrix(a p.s.d matrix that can be represented as $Z^TZ$) is computed from euclidean distance between rows in $X$ and the following is minimized over $Y$. minimize: $||G-Y^TY||_{F}^{2}$.
A: PCA yields the EXACT same results as classical MDS if Euclidean distance is used.
I'm quoting Cox & Cox (2001), p 43-44:

There is a duality between a principals components analysis and PCO [principal coordinates analysis, aka classical MDS] where dissimilarities are given by Euclidean distance.

The section in Cox & Cox explains it pretty clearly:


*

*Imagine you have $X$ = attributes of $n$ products by $p$ dimensions, mean centered

*PCA is attained by finding eigenvectors of the covariance matrix ~ $X'X$ (divided by n-1) -- call the eigenvectors $\xi$, and eigenvalues $\mu$.

*MDS is attained by first converting $X$ into distance matrix, here, Euclidean distance, i.e., $XX'$, then finding the eigenvectors -- call the eigenvectors $v$, and eigenvalues $\lambda$.

*p 43: "It is a well known result that the eigenvalues of $XX'$ are the same as those for $X'X$, together with an extra n-p zero eigenvalues."  So, for $i < p$, $\mu_i$ = $\lambda_i$

*Going back to definition of eigenvectors, consider the $i^{th}$ eigenvalues.  $X'Xv_i = \lambda_i v_i$

*Premultiply $v_i$ with $X'$, we get $(X'X)X'v_i = \lambda_i X'v_i$

*We also have $X'X \xi_i = \mu_i \xi_i$.  Since $\lambda_i = \mu_i$, we get that $\xi_i = X'v_i$ for $i<p$.

A: Classic Torgerson's metric MDS is actually done by transforming distances into similarities and performing PCA (eigen-decomposition or singular-value-decomposition) on those. [The other name of this procedure (distances between objects -> similarities between them -> PCA, whereby loadings are the sought-for coordinates) is Principal Coordinate Analysis or PCoA.] So, PCA might be called the algorithm of the simplest MDS.
Non-metric MDS is based on iterative ALSCAL or PROXSCAL algorithm (or algorithm similar to them) which is a more versatile mapping technique than PCA and can be applied to metric MDS as well. While PCA retains m important dimensions for you, ALSCAL/PROXSCAL fits configuration to m dimensions (you pre-define m) and it reproduces dissimilarities on the map more directly and accurately than PCA usually can (see Illustration section below).
Thus, MDS and PCA are probably not at the same level to be in line or opposite to each other. PCA is just a method while MDS is a class of analysis. As mapping, PCA is a particular case of MDS. On the other hand, PCA is a particular case of Factor analysis which, being a data reduction, is more than only a mapping, while MDS is only a mapping.
As for your question about metric MDS vs non-metric MDS there's little to comment because the answer is straightforward. If I believe my input dissimilarities are so close to be euclidean distances that a linear transform will suffice to map them in m-dimensional space, I will prefer metric MDS. If I don't believe, then monotonic transform is necessary, implying use of non-metric MDS.

A note on terminology for a reader. Term Classic(al) MDS (CMDS) can have two different meanings in a vast literature on MDS, so it is ambiguous and should be avoided. One definition is that CMDS is a synonym of Torgerson's metric MDS. Another definition is that CMDS is any MDS (by any algorithm; metric or nonmetric analysis) with single matrix input (for there exist models analyzing many matrices at once - Individual "INDSCAL" model and Replicated model).

Illustration to the answer. Some cloud of points (ellipse) is being mapped on a one-dimensional mds-map. A pair of points is shown in red dots.

Iterative or "true" MDS aims straight to reconstruct pairwise distances between objects. For it is the task of any MDS. Various stress or misfit criteria could be minimized between original distances and distances on the map: $\|D_o-D_m\|_2^2$, $\|D_o^2-D_m^2\|_1$, $\|D_o-D_m\|_1$. An algorithm may (non-metric MDS) or may not (metric MDS) include monotonic transformation on this way.
PCA-based MDS (Torgerson's, or PCoA) is not straight. It minimizes the squared distances between objects in the original space and their images on the map. This is not quite genuine MDS task; it is successful, as MDS, only to the extent to which the discarded junior principal axes are weak. If $P_1$ explains much more variance than $P_2$ the former can alone substantially reflect pairwise distances in the cloud, especially for points lying far apart along the ellipse. Iterative MDS will always win, and especially when the map is wanted very low-dimensional. Iterative MDS, too, will succeed more when a cloud ellipse is thin, but will fulfill the MDS-task better than PCoA. By the property of the double-centration matrix (described here) it appears that PCoA minimizes $\|D_o\|_2^2-\|D_m\|_2^2$, which is different from any of the above minimizations.
Once again, PCA projects cloud's points on the most advantageous all-corporal saving subspace. It does not project pairwise distances, relative locations of points on a subspace most saving in that respect, as iterative MDS does it. Nevertheless, historically PCoA/PCA is considered among the methods of metric MDS.
