Probability that no couple sits together in a circle Suppose $n$ couples are seated in a circle, and let $P_n$ be the probability that no couple is sitting together.
Using Inclusion-Exclusion, I believe it can be shown that 
$\hspace{.2 in}\displaystyle P_n=1-\sum_{i=1}^n (-1)^{i+1}\binom{n}{i}2^i\frac{(2n-1-i)!}{(2n-1)!},$
and I would like to find out how to prove that  $\displaystyle\lim_{n\to\infty}P_n=\frac{1}{e}\;\;$ (or show that this is not the case).

Here are some numerical values:
$P_3\approx.267,\;P_4\approx.295,\;P_5\approx.310\;, P_6\approx.320,\;P_7\approx.327,\;P_8\approx.332,\;P_9\approx.336,\;\;P_{10}\approx.340$
For a related question, see Showing probability no husband next to wife converges to $e^{-1}$
 A: Verification of the Sum
How many ways for $k$ of the $n$ couples to sit together:
$$
\overbrace{2^kk!\binom{n}{k}}^1\left[\vphantom{\binom{2n-k}{k}}\right.\overbrace{\binom{2n-k}{k}}^2+\overbrace{\binom{2n-k-1}{k-1}}^3\left.\vphantom{\binom{2n-k}{k}}\right]\overbrace{(2n-2k)!\vphantom{\binom{2n-k}{k}}}^4\tag{1}
$$
Explanation:
$1$: the number of ways to select the $k$ couples: $\binom{n}{k}$
$\phantom{1\text{: }}$times the number of ways to arrange them among themselves: $2^kk!$
$2$: the number of ways to arrange $k$ couples and $2n-2k$ singles
$3$: the number of ways to arrange $k-1$ couples and $2n-2k$ singles
$\phantom{1\text{: }}$this is the number of ways to have a couple across the end of the line
$4$: the number of ways to arrange the $2n-2k$ singles among themselves
Noting that $\binom{2n-k}{k}+\binom{2n-k-1}{k-1}=\frac{2n}{2n-k}\binom{2n-k}{k}$, Inclusion-Exclusion gives the probability that no couple sits together to be
$$
\begin{align}
&\frac1{(2n)!}\sum_{k=0}^n(-1)^k2^kk!\binom{n}{k}\binom{2n-k}{k}\frac{2n}{2n-k}(2n-2k)!\\
&=\sum_{k=0}^n(-2)^k\binom{n}{k}\frac{(2n-k-1)!}{(2n-1)!}\tag{2}
\end{align}
$$

Verification of the Limit
$$
\begin{align}
\sum_{k=0}^n(-2)^k\binom{n}{k}\frac{(2n-k-1)!}{(2n-1)!}
&=\sum_{k=0}^n\frac{(-1)^k}{k!}\frac{2^k\frac{n!}{(n-k)!}}{\frac{(2n-1)!}{(2n-k-1)!}}\\
&=\sum_{k=0}^n\frac{(-1)^k}{k!}\frac{2n(2n-2)\cdots(2n-2k+2)}{(2n-1)(2n-2)\cdots(2n-k)}\\
&\to\sum_{k=0}^\infty\frac{(-1)^k}{k!}\\[3pt]
&=\frac1e\tag{3}
\end{align}
$$
A: Start   with  the  inclusion-exclusion   argument.   We   compute  the
cardinality of a set of arrangements  where a given set of $q$ couples
or more  sit together. There  are $2^q$ possibilities of  ordering the
$q$ couples.  These are now fused and we have $2n-2q$ people left over
which are  joined by $q$ couples  in a circular  permutation, so there
are $(2n-2q+q)!/(2n-2q+q)  = (2n-q-1)!$ such  arangements. This yields
by inclusion-exclusion the closed formula
$$\sum_{q=0}^n {n\choose q} (-1)^q 2^q (2n-q-1)!$$
which  we  divide  by  $(2n-1)!$  to  get  the  probability.  For  the
asymptotics we use
$$(2n-q-1)! = \Gamma(2n-q) = \int_0^\infty x^{2n-q-1} e^{-x} dx.$$
This yields for the sum
$$\int_0^\infty x^{2n-1} e^{-x} 
\sum_{q=0}^n {n\choose q} (-1)^q 2^q x^{-q}
dx
\\ = \int_0^\infty x^{2n-1} e^{-x} 
\left(1-\frac{2}{x}\right)^n dx
= \int_0^\infty x^{n-1} (x-2)^n e^{-x} dx.$$
We     will     compute     the     asymptotic     with     Laplace's
method.  To  do  this
put $x=nz$ to obtain
$$n^{n-1} \times n^n \times n 
\times \int_0^\infty z^{n-1} (z-2/n)^n e^{-nz} dz
\\ = n^{2n}
\times \int_0^\infty z^{n-1} (z-2/n)^n e^{-nz} dz
.$$
The integrand is
$$\exp(n(\log(z)+\log(z-2/n)-\log(z)/n-z)) = \exp(n f(z)).$$
To find the saddle point we solve the saddle point equation $f'(z)=0$
which produces $$z_0 \approx 2 + \frac{1}{n^2}.$$
The asymptotic is then given by
$$\sqrt{\frac{2\pi}{n|f''(z_0)|}} e^{nf(z_0)}.$$
We have for $f(z_0)$ 
$$\log 2 + \log(1+1/n^2/2) 
+ \log 2 + \log(1-1/n+1/n^2/2)
\\ - 1/n \log 2 - 1/n \log (1 + 1/n^2/2) 
- 2 - 1/n^2.$$
The first two terns are
$$f(z_0) \approx 2\log 2 -2 - 1/n \log 2 - 1/n.$$
Furthermore we have
$$f''(z_0) \approx -\frac{1}{2}$$
so $f(z)$ does indeed have a maximum there.
Putting it all together we obtain
$$\sqrt{\frac{4\pi}{n}}
e^{2n\log 2 - 2n - \log 2 - 1}
= 2\sqrt{\frac{\pi}{n}}
\frac{2^{2n}}{2e^{2n+1}}
= \sqrt{\frac{\pi}{n}}
\frac{2^{2n}}{e^{2n+1}}.$$
This yields for our sum
$$n^{2n-1/2} \sqrt{\pi} \frac{2^{2n}}{e^{2n+1}}.$$
Now Stirling's formula yields
$$(2n-1)! \sim \sqrt{2\pi(2n-1)} (2n-1)^{2n-1} e^{-(2n-1)}.$$
This gives for the probability
$$\frac{1}{e^2} \frac{2\sqrt{n}}{\sqrt{2(2n-1)}}
\frac{2^{2n-1} n^{2n-1}}{(2n-1)^{2n-1}}
\sim \frac{1}{e^2}
\left(1+\frac{1}{2n-1}\right)^{2n-1}
\sim \frac{1}{e}$$
and the conjecture is proved.
Remark.  These  are  the  details  of  the  inclusion-exclusion
argument.  The  nodes of  the underlying poset  represent arrangements
where some set of $q$  couples sit together, possibly more. The Mobius
function of the poset is $(-1)^q.$ A specific arrangement with exactly
$p$ couples sitting  together is included in the  ${p\choose q}$ nodes
where $q$  couples from the $p$  couples sit together, or  more. It is
not included in  any nodes representing more than  $p$ couples sitting
together  because  this is  impossible  when  the  count is  precisely
$p$. Therefore the total weight with the given Mobius function is
$$\sum_{q=0}^p {p\choose q} (-1)^q.$$
This  is zero  when  $p\ge 1$  and one  when  $p =  0.$ These  weights
correctly represent  the requirement  of arrangements with  no couples
sitting together being counted once  and all others being counted zero
times.
Addendum.  An approach  to this  problem using  the  methods of
experimental mathematics might begin by computing the values for small
$n$  using  total  enumeration  and  hence  verify  the  formula  from
inclusion-exclusion, and  then use  that to compute  the probabilities
for  larger $n$,  leading  to a  conjecture  that we  have $1/e.$  The
enumeration  task  is practicable  to  about  $n=5$  and confirms  the
correctness  of the  formula. This  is the  code, which  consumes very
little memory. We get starting with one couple
$$0, 2, 32, 1488, 112512,\ldots$$
Observe  that we  got  the correct  value  for $n=1$  unlike from  the
formula,  which does  not work  here because  when there  is  only one
couple the  two orderings  of its constituents  are the same,  but the
factor $2^q$ produces a value of two instead of a value of one.

with(combinat);

X :=
proc(n)
option remember;
local perm, circ, count, pos;

    count := 0;

    perm := firstperm(2*n-1);

    while type(perm, `list`) do
        circ := [op(perm), 2*n, perm[1]];

        for pos to 2*n do
            if abs(circ[pos]-circ[pos+1]) = 1
            and type(min(circ[pos], circ[pos+1]), `odd`)
            then
                break;
            fi;
        od;

        if pos = 2*n + 1 then
            count := count + 1;
        fi;

        perm := nextperm(perm);
    od;

    count;
end;


Q := n -> add(binomial(n,q)*(-1)^q*2^q*(2*n-q-1)!, q=0..n);

