The probability is $e^{-2n (\log \log n + 1 - o(1))}$. $\DeclareMathOperator{\kernel}{kernel} \DeclareMathOperator{\range}{range}$
Lower Bound
The above bound is achieved by random $f$ and $g$ with $\range(f) ⊆ \operatorname{kernel}(g)$ and $\range(g) ⊆ \kernel(f)$ (i.e. $∀x \, f(g(x))=g(f(x))=0$); typically for such $f$ and $g$, $|\kernel(f)| ≈ |\kernel(g)| ≈ n/\log n$ (up to $1±o(1)$ factor).
Commutativity is extremely unlikely absent special patterns, and (as we show below) the most likely is (approximately) the above annihilation pattern.
Besides choosing the value of "zero", the one change I have to improve the probability is to have $≈\log^2 n$ exceptions for which commutativity happens accidently: For a typical input $x$, there are $≈n^2/\log^2 n$ ordinary choices of $(f(x),g(x))$, and $≈n$ choices for accidental commutativity (i.e. $f(g(x))≠0$). I suspect that this model captures $1-o(1)$ fraction of commuting $f$ and $g$.
With this model, we can compute the asymptotics quite precisely, though it is easy to make a mistake about a constant or otherwise small factor. I came up with the following probability, where $h(p) = -p \log p - (1-p) \log(1-p)$ (entropy for probability $p$):
$(1±o(1)) \max\limits_{n^{-1/3}≤p<1} p^3 n e^{ -2n(p \log n - (1-p) \log p - h(p)) \, + \, (1-p)^4 p^{-2} \, - \, 2(1/p-1)}$
Intuition for $p$: $p ≈ |\kernel(f)|/n ≈ |\kernel(g)|/n ≈ \log^{-1} n$.
Explanation of all factors:
$n$ - number of choices for "zero"
$p = (p^{1/2})^2$ - stems from the error tolerance (stdev) for $p$ (computed using the second derivatives) being $≈p^{-1/2}$ times smaller than $\sqrt{p/n}$; the square is because the correction is for both $f$ and $g$ (and the covariance of the errors is small enough).
$p^2$ - both kernels include 0
$e^{-2(1/p-1)} = n^{-2±o(1)}$ - probability that $f$ and $g$ are never 0 outside of what we designated as their kernels
$e^{-2n p \log n}$ - kernel
$e^{-2n(1-p) \log p}$ - range restriction
$e^{2h(p)n}$ - effective number of choices for the kernels (before the $p^3$ correction)
$e^{(1-p)^4 p^{-2}} ≈ (1 + \frac{(n(1-p))^2/n}{(pn)^2})^{n(1-p)^2} = n^{(1±o(1))\log n}$ is the correction for the exceptions with accidental commutativity; $≈(1-p)^2$ is the probability of being outside both kernels; $(n(1-p))^2/n$ is accidental commutativity ($n(1-p)$ choices each times $n^{-1}$ collision chance); $n^{-1/3}≤p$ (chosen arbitrarily) more than suffices for the accuracy of the approximation by the exponential.
Numerical results: For $n≤20$, the exact values (times $n^{2n}$) are published as an OEIS sequence A181162 and are in an order-of-magnitude agreement with the predicted values (the above formula underestimates in about 4-5 times there, but (for $n>10$) improving with $n$). This testing is limited because $p$ goes to zero very slowly.
Upper Bound
For a mapping $f$, the probability of commutativity with a random mapping $g$ is at most $\min_{S⊆\operatorname{domain}(f)} n^{-|f[S] \setminus S|}$ since for commuting $f$ and $g$, $g\restriction (f[S] \setminus S)$ is uniquely determined by $f$ and $g \restriction S$.
Thus, for random commuting mappings, with high probability, $|\range(f)|$ and $|\range(g)|$ are both $O(n / \log n ⋅ \log \log n)$ (we will use this below). This range restriction has probability $e^{-(2±o(1)) n \log \log n}$. If $|\range(f)|$ were too large here, then unless $f$ is identity on too many points to be likely, we get a large $|f[S] \setminus S|$ and thus commute probabilities well below the lower bound.
We now get a more precise upper bound. For $x∈\range(g)$ define the equivalence class $[x]_f = \{y∈\range(g): f(x)=f(y)\}$, and analogously define $[x]_g$ (for $x∈\range(f)$). If $f$ and $g$ commute, we get a corresponding bijection $P$ between the set of all $[x]_f$ and the set of all $[x]_g$: $P([g(x)]_f) = [f(x)]_g$. Let $S = \{x∈\range(g): |[x]_f|+|P([x]_f)|≥ne^{-\sqrt{\log n}}\}$, and set $T = P[S] = ∪\{P([x]_f):x∈S\}$ (here, $e^{-\sqrt{\log n}}$ can be arbitrary within bounds). Let $m = |\range(f)∪\range(g)|$.
The probability that random $f$ and $g$ commute is at most $\sup_{|S|≤m, \, |T|≤m, \, m=O(n / \log n ⋅ \log \log n)} e^{o(n)} n^{-|S|-|T|} (\frac{|S||T|}{n^2} + e^{-\sqrt{\log n}})^{n-m} = $ $e^{-2n (\log \log n + 1 - o(1))}$ as desired. Here, $n^{-|S|-|T|}$ is the probability that $f$ and $g$ both have few enough values on $S$ and $T$ respectively (ignoring $e^{o(n)}$ factors). Next, after setting $f$ and $g$ on inputs in $\range(f)∪\range(g)$, we test $f(g(x))=g(f(x))$ for $x$ outside of $\range(f)∪\range(g)$ (hence the $n-m$ exponent) with $f$ and $g$ random on $x$. We either hit a large equivalence class (the $|S||T|$ factor) or deal with the low likelihood that $f$ and $g$ hit matching classes (hence the $e^{-\sqrt{\log n}}$ above).
Furthermore, the $|S||T|$ above is an overestimate unless there is a single dominant equivalence class: The $|S||T|$ can be refined to $\sum_i |S_i| |P(S_i)|$ where $S_i⊆S$ is an equivalence class (and with a further reduction if $|S_i||P(S_i)|$ does not reflect the probability of $f(x)∈S_i$). To match the $e^{-2n (\log \log n + 1 - o(1))}$ bound, $|S|$, $|T|$, and for $n \! - \! o(n)$ $x$, $|[g(x)]_f|$ and $|[f(x)]_g|$ must all be $(1±o(1))n/\log n$.
Thus, with high probability, for random commuting $f$ and $g$, the composition $fg$ is constant, except on $o(n)$ inputs.
Related Problems
As you note, the probability that a random permutation $f$ commutes with a random mapping $g$ is $Θ(n^{-n-1/2} ρ^{-n})$ for $ρ≈0.29224$. This probability should be dominated by $f$ for which the number of $k$-cycles decreases roughly exponentially with $k$. Here is a derivation of the weaker $e^{Θ(n)}/n!$ bound.
Briefly, for commutativity, a $k$-cycle of $f$ must be mapped to a cycle whose length divides $k$, and the total number of commuting mappings $g$ is $\prod_k (\sum_{q|k} qc_q)^{c_k}$ where $c_k$ is the number of $k$-cycles of $f$ (so $\sum_k kc_k=n$). We get the $e^{Θ(n)}/n!$ lower bound by considering cycles of length 1 and 2. For the $e^{Θ(n)}/n!$ upper bound, if $f$ has $m$ cycles, the number of choices of $f$ is (at most) $e^{O(n)} n^{n-m}$, while the commute probability with a random $g$ is $≤n^{-(n-m)}$.
One can also ask about the likelihood of commutativity in the absence of special patterns:
- For almost all mappings $f$, I think the number of commuting mappings $g$ is $e^{Θ(n)}$: Once we set $g(x)$ for a given $x$, we can go back through all values of $f^{-1}(x)$, $f^{-2}(x)$, ..., and typically have $O(1)$ choices for $g$ per point.
- For almost all mappings $f$, the number of commuting permutations $g$ is $e^{Θ(n)}$. For the lower bound, for every $x∉f[\range(f)]$, one can arbitrarily permute $f^{-1}(x)$. The upper bound follows from the bound on the probability that a random permutation and a random mapping commute.
- For almost all permutations $f$, the number of commuting mappings $g$ should be $n^{(\frac{1}{2}±o(1)) \log n}$: A random permutation $f$ has an expected $≈\log n$ cycles, with an average of $1/k$ cycles of length $k$; and a large cycle of $f$ should be typically mapped to itself by $g$. Furthermore, for every constant $ε>0$, for $1-ε$ fraction of permutations $f$, a $Θ(1)$ fraction of commuting $g$ should be permutations.
