$P(X_1 < \min(X_i,\ldots, X_n))$ across different normal random variables I have a set of mutually independent normal distributions $X_1$ to $X_5$ (with means and standard deviations) which represent finishing times for swimmers over a certain distance. The actual data is as follows:
$$X_1(60, 3.0)$$
$$X_2(61, 1.0)$$
$$X_3(58, 2.3)$$
$$X_4(63, 2.4)$$
$$X_5(61, 1.7)$$
So swimmer 1 ($X_1$) has a mean finishing time of 60 seconds with a standard deviation of 3.0 seconds.
Question 1: What is the probability of an event where $X_i$ finishes first. e.g.
$$P(X_1 \lt X_i, i=2,\ldots,n)$$
Question 2: If I calculate this for all swimmers, can I simply order the results to determine the most probable finishing order?
This is not homework.
Based on the answers to this Cross Validated question, I have tried to solve this problem based on the first answer. i.e.
$$\Pr(X_1 \le X_i, i=2,\ldots,n) = \int_{-\infty}^{\infty} \phi_1(t) [1 - \Phi_2(t)]\cdots[1 - \Phi_n(t)]dt$$
Where  $\phi_i$ is the PDF of $X_i$ and $\Phi_i$ is its CDF.
Based on this formula, the results I obtained were:
$$\Pr(X_1 \le X_i, i=2\ldots5) = 0.259653$$
$$\Pr(X_2 \le X_i, i=1,3\ldots5) = 0.214375$$
$$\Pr(X_3 \le X_i, i=1\ldots2, 4\ldots5) = 0.611999$$
$$\Pr(X_4 \le X_i, i=1\ldots3, 5) = 0.0263479$$
$$\Pr(X_5 \le X_i, i=1\ldots4) = 0.0697597$$
However, the probabilities add to 1.182135 when they should add to 1.0. I’m not sure if the formula is incorrect or my implementation of the integral (I used Excel and the trapezoidal method).
I also attempted to solve the problem using Dillip’s method (from the above mentioned question) as follows:
\begin{align*}
P(X_1 < \max X_i) &=  P\{(X_1 < X_2) \cup \cdots \cup (X_1 < X_n)\\
&\leq \sum_{i=2}^n P(X_1 < X_i)\\
&= \sum_{i=2}^n Q\left(\frac{\mu_1 - \mu_i}{\sqrt{\sigma_1^2 + \sigma_i^2}}\right)
\end{align*}
However, the probability results were much greater than 1 in most cases so abandoned this approach. By the way, what exactly does $\max X_i$ mean?
Any assistance in calculating the probability would be appreciated.
 A: The answers above don't provide a closed-form solution; nor do the ones to the related question here. I will try to give an analytical answer. To fix notation, let:
$X_i$ be the independent random variables, with $X_i \sim N(\mu_i,\sigma_i^2)$
$Y_{i-1}:= X_1 - X_i, i =2,\ldots, n$. 
Finally, let $e\in R^{n-1}$ be such that $(e)_i=1$.
Then
$P(X_1 \le X_i, i=2,\ldots,n) = P(Y_i \ge 0, i=2,\ldots,n)$
It is $Y_i  \sim N(\mu_1-\mu_i, \sigma_1^2+\sigma_i^2)$ and moreover (straightforward calculation) $\text{cov} (Y_i, Y_j)=\sigma_1^2$.
so the covariance matrix $\Sigma$ of $Y$ is given by
$\Sigma = \text{diag}(\sigma_2^2,\ldots, \sigma_n^2) + \sigma_1^2 ee'$
Define $\nu_{i-1} := \mu_i-\mu_1, i =2,\ldots, n$. Now we use the standard affine transformation of a multivariate standard normal to obtain an arbitrarily distributed multivariate normal. Let $\xi \sim N(0, I)$ be a $(n-1)$-dimensional standard normal. We have $Y =_\text{dist} \nu + Q^{1/2} \xi$. 
The square root $Q^{1/2}$ is uniquely identified.
The probability we want to estimate becomes
$P(Y \ge 0) = P(\xi \ge -Q^{-1/2} \nu) = \prod_{i=1}^{n-1} \bar \Phi( (-Q^{-1/2} \nu)_i)$
where the last equality follows from the independence of the $\xi_i$ an d $\bar \Phi$ is the complement of the cumulative distribution of the standard normal.
Last observation: $Q^{1/2}$ is trivial in the case where $\kappa:=\sigma_2=\ldots=\sigma_n$. I am omitting the derivation, but point out that there is one eigenvector $e$ with eigenvalue $\kappa^2+\sigma_1^2$. The other eigenvector lie in $[e]^\perp$ and have all eigenvalues equal to $\kappa^2$. 
A: You had two questions

Question 1: What is the probability of an event where $X_i$ finishes first.

Your proposed solution look correct but, as you say, you clearly have an error in implelentation as the probablities do not add up to $1$. M. Berk has shown that it is likely to be an issue with Swimmer 2.

Question 2: If I calculate this for all swimmers, can I simply order the results to determine the most probable finishing order? 

Not quite - if you have two swimmers with the same mean time in the middle of the group then the chance of one beating the other is $\frac12$, but the one with the higher standard deviation is more likely to be the overall winner: in your particular example, the chance of $X_2$ beating $X_5$ is $0.5$ but $X_5$ is more likely than $X_2$ to be the overall winner, so $X_5$ is less likely than $X_2$ to be third because of its higher standard deviation.
Using simulation, the most likely finishing order is $X_3,X_1,X_2,X_5,X_4$ (with a probability about 9.0%) above $X_3,X_1,X_5,X_2,X_4$ (about 8.0%), $X_1,X_3,X_2,X_5,X_4$ (about 6.1%), $X_3,X_2,X_5,X_1,X_4$ (about 5.8%), $X_1,X_3,X_5,X_2,X_4$ (also about 5.8%), $X_3,X_5,X_1,X_2,X_4$ (about 4.5%) and other less likely outcomes.  Your idea would have predicted 31524, the second most likely outcome. 
A: The formula at https://stats.stackexchange.com/a/44142/919 is correct.
Below is R code to implement it (as a function p.max). The code to input the data of the question and apply p.max is short:
X <- data.frame(Mean = c(60, 61, 58, 63, 61),
                SD = c(3, 1, 2.3, 2.4, 1.7))
X$p <- sapply(seq_len(nrow(X)), p.max, mu=-X$Mean, sigma=X$SD)

(Note how supplying the negatives of the times enables p.max to compute the chances of being the smallest.)
This plot illustrates the five densities, colored according to the winning chances:

The sum of the computed chances is within $10^{-10}$ of equaling $1.$
#
# Given X[i] ~ Normal(mu[i], sigma[i]),
# compute the chance that the jth component is largest.
# `...` is passed to `integrate`.
#
p.max <- function(j, mu, sigma, Z=7, ...) {
  if (is.infinite(mu[j])) {
    if (mu[j] < 0) {
      if (all(is.infinite(mu) & mu < 0)) return(1 / length(mu))
      return(0)
    }
    return(1 / sum(is.infinite(mu) & mu > 0))
  }
  #
  # Choose units that make X[j] standard Normal.
  #
  mu <- (mu[-j] - mu[j]) / sigma[j]
  sigma <- sigma[-j] / sigma[j]
  #
  # Set integration limits.
  #
  xlim <- range(c(-Z, Z, mu-Z*sigma, mu+Z*sigma))
  #
  # Integrate.
  #
  h <- Vectorize(function(x)
    exp(dnorm(x, log=TRUE) + sum(pnorm(x, mu, sigma, log.p=TRUE))))
  integrate(h, xlim[1], xlim[2], ...)$value
}

