I came across this question while preparing a talk on this paper. There we gave the explicit formula for $p_n$ (in the GOE) that the other answers anticipated could be found by random matrix methods.

The formula is given as follows:

First, let $n\geq 1$ be any integer, and define $n':=2 \lceil n/2 \rceil$. When real symmetric matrices are chosen according to the $n$-dimensional Gaussian Orthogonal Ensemble, the probability of positive definiteness is given by: \begin{equation}p_n=\frac{{\mathrm{ Pf}}(A)}{2^{n(n+3)/4}\prod_{m=1}^n\Gamma(\frac m 2)},\label{theorem3} \end{equation} where $A$ is the $n'\times n'$ skew-symmetric matrix whose $(i,j)$-entry $a_{ij}$ is given for $i<j$ by: \begin{align}\label{pfaffianformula} a_{ij} = \begin{cases} \,2^{i+j-2}\Gamma(\frac{i+j}{2})\left(\beta_{\frac{1}{2}}(\frac{i}{2},\frac{j}{2})- \beta_{\frac{1}{2}}(\frac{j}{2},\frac{i}{2})\right) & \mbox{if }~ i<j\leq n, \\[.075in] \,2^{i-1}\Gamma(\frac{i}{2}) & \mbox{if }~ i<j=n+1.\\[.075in] \end{cases} \end{align} $($Note that the second case arises only when $n$ is odd.$)$

Here $\mathrm{Pf}$ is the Pfaffian, $\Gamma$ denotes the usual gamma function $\Gamma(s)=\int_{0}^\infty x^{s-1}e^{-x} \mathrm{d} x$, and $\beta_{t}$ is the usual incomplete beta function $\beta_{t}(i,j) = \int_0^{t} x^{i-1}(1-x)^{j-1} \mathrm{d} x.$ This formula follows from de Bruijn's identity for evaluating determinantal integrals.

As for the asymptotics, the leading term of $\log p_n$ is indeed $-\frac{\log(3)}{4} n^2$ as noted in Mikael's answer. If you are interested in the leading order asymptotic for $p_n$ it can be found from the results of Borot, Eynard, Majumdar and Nadal. Specialising to the GOE they give: \begin{align} p_n \sim \exp\left(-\frac{\log3}{4} n^2 - \frac{1}{2}\log\left(1+\frac{2}{\sqrt{3}}\right)n - \frac{1}{24}\log n + \alpha_1 \right)(1+o(1)).\label{asympexp} \end{align} The constant $\alpha_1$ is given by: \begin{align} \alpha_1= \frac{\zeta'(-1)}{2} - \frac{\log(2)}{12} -\frac{\log(3)}{16}+ \frac{1}{4}\log\left(1+\frac{2}{\sqrt{3}}\right)\simeq-0.0172, \end{align} where $\zeta(z)$ is the Riemann zeta function.

I came across this question while preparing a talk on this paper. There we gave the explicit formula for $p_n$ (in the GOE) that the other answers anticipated could be found by random matrix methods.

The formula is given as follows:

First, let $n\geq 1$ be any integer, and define $n':=2 \lceil n/2 \rceil$. When real symmetric matrices are chosen according to the $n$-dimensional Gaussian Orthogonal Ensemble, the probability of positive definiteness is given by: \begin{equation}p_n=\frac{{\mathrm{ Pf}}(A)}{2^{n(n+3)/4}\prod_{m=1}^n\Gamma(\frac m 2)},\label{theorem3} \end{equation} where $A$ is the $n'\times n'$ skew-symmetric matrix whose $(i,j)$-entry $a_{ij}$ is given for $i<j$ by: \begin{align}\label{pfaffianformula} a_{ij} = \begin{cases} \,2^{i+j-2}\Gamma(\frac{i+j}{2})\left(\beta_{\frac{1}{2}}(\frac{i}{2},\frac{j}{2})- \beta_{\frac{1}{2}}(\frac{j}{2},\frac{i}{2})\right) & \mbox{if }~ i<j\leq n, \\[.075in] \,2^{i-1}\Gamma(\frac{i}{2}) & \mbox{if }~ i<j=n+1.\\[.075in] \end{cases} \end{align} $($Note that the second case arises only when $n$ is odd.$)$

Here $\mathrm{Pf}$ is the Pfaffian, $\Gamma$ denotes the usual gamma function $\Gamma(s)=\int_{0}^\infty x^{s-1}e^{-x} \mathrm{d} x$, and $\beta_{t}$ is the usual incomplete beta function $\beta_{t}(i,j) = \int_0^{t} x^{i-1}(1-x)^{j-1} \mathrm{d} x.$ This formula follows from de Bruijn's identity for evaluating determinantal integrals.

As for the asymptotics, the leading term of $\log p_n$ is indeed $-\frac{\log(3)}{4} n^2$ as noted in Mikael's answer. If you are interested in the leading order asymptotic for $p_n$ it can be found from the results of Borot, Eynard, Majumdar and Nadal. Specialising to the GOE they give: \begin{align} p_n \sim \exp\left(-\frac{\log3}{4} n^2 - \frac{1}{2}\log\left(1+\frac{2}{\sqrt{3}}\right)n - \frac{1}{24}\log n + \alpha_1 \right)(1+o(1)).\label{asympexp} \end{align} The constant $\alpha_1$ is given by: \begin{align} \alpha_1= \frac{\zeta'(-1)}{2} - \frac{\log(2)}{12} -\frac{\log(3)}{16}+ \frac{1}{4}\log\left(1+\frac{2}{\sqrt{3}}\right)\simeq-0.0172, \end{align} where $\zeta(z)$ is the Riemann zeta function.

Edit: The exact formula lends itself easily to numerics. I have added a plot of $\log(p_n)$ for $1\leq n \leq 400.$ Fitting the data for $201\leq n \leq 400$ to $a n^2 + b n + c\log n +d$ we find $$−0.27465n^2 - 0.38382n -0.041504 \log n -0.018061.$$ A numerical approximation for the asymptotic formula is $$−0.27465n^2 - 0.38382n -0.041667 \log n -0.017223,$$ so we see good agreement.

$\log p_n$ for for $1\leq n \leq 400.$" />