There is a closed form for this probability, but it is not particularly concise. To obtain the formula, start with the probability of the event of interest, conditional on the event $X_1+X_2=r$. It can easily be shown that:
$$X_1 | X_1 + X_2 = r, n, \mathbf{p} \sim \text{Bin} \Big( r, \frac{p_1}{p_1+p_2} \Big).$$
Under the stated conditioning event, the event $X_1 \geqslant X_2$ is equivalent to the event $X_1 \geqslant r/2$. We therefore have the conditional probability:
$$\begin{equation} \begin{aligned}
\mathbb{P}(X_1 \geqslant X_2 | X_1 + X_2 = r, n, \mathbf{p})
&= \mathbb{P}(X_1 \geqslant r-X_1 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt]
&= \mathbb{P}(X_1 \geqslant r/2 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt]
&= \mathbb{P}(X_1 \geqslant \lceil r/2 \rceil | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt]
&= \frac{1}{(p_1+p_2)^r} \sum_{x = \lceil r/2 \rceil}^r {r \choose x} \cdot p_1^x \cdot p_2^{r-x}. \\[6pt]
\end{aligned} \end{equation}$$
In practice, the upper tail of the binomial distribution is usually computed using the regularised incomplete beta function, rather than by summing the mass function. In this latter form, we have:
$$\begin{equation} \begin{aligned}
\mathbb{P}(X_1 \geqslant X_2 | X_1 + X_2 = r, n, \mathbf{p})
&= 1-\mathbb{P}(X_1 \leqslant \lceil r/2 \rceil - 1 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt]
&= 1-I_{p_2/(p_1+p_2)}(r - \lceil r/2 \rceil + 1, \lceil r/2 \rceil). \\[6pt]
\end{aligned} \end{equation}$$
Applying the law of total probability gives the marginal probability of interest:
$$\begin{equation} \begin{aligned}
\mathbb{P}(X_1 \geqslant X_2 | n, \mathbf{p})
&= \sum_{r=0}^n \mathbb{P}(X_1 + X_2 = r | n) \cdot \mathbb{P}(X_1 \geqslant X_2 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt]
&= \sum_{r=0}^n \text{Bin}(r|n, p_1+p_2) \cdot \mathbb{P}(X_1 \geqslant X_2 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt]
&= \sum_{r=0}^n {n \choose r} \cdot (1-p_1+p_2)^{n-r} \sum_{x = \lceil r/2 \rceil}^r {r \choose x} \cdot p_1^x \cdot p_2^{r-x} \\[6pt]
&= \sum_{r=0}^n {n \choose r} \cdot (p_1+p_2)^{r} \cdot (1-p_1+p_2)^{n-r} \\[6pt]
&\quad \quad \quad \times (1-I_{p_2/(p_1+p_2)}(r - \lceil r/2 \rceil + 1, \lceil r/2 \rceil)). \\[6pt]
\end{aligned} \end{equation}$$
Either of the last two lines give a closed form solution to the problem. (The first is a sum over relevant values of the mass function of the binomial. The second uses the regularised incomplete beta function for the distribution. The latter is closer to what you would use when implementing the function in statistical computing.) Either expression is somewhat cumbersome for large $n$.
Computing the probability: We can construct a function in R
to compute this probability using the standard statistical distribution functions. (Since this computation involves adding very small probabilities we work in log-space for the computation, and then convert back to a standard probability at the end.)
PROB <- function(n, p, log = FALSE) {
#Normalise probability vector
p <- p/sum(p);
#Computeprobability
LOGTERMS <- rep(0, n+1);
LOGTERMS[1] <- dbinom(0, size = n, prob = p[1]+p[2], log = TRUE);
CLOGPROB <- rep(0, n+1);
for (r in 1:n) {
CLOGPROB[r+1] <- pbinom(ceiling(r/2)-1, size = r, prob = p[1]/(p[1]+p[2]),
lower.tail = FALSE, log.p = TRUE);
LOGTERMS[r+1] <- dbinom(r, size = n, prob = p[1]+p[2], log = TRUE) + CLOGPROB[r+1]; }
LOGPROB <- matrixStats::logSumExp(LOGTERMS);
if (log) { LOGPROB } else { exp(LOGPROB) } }
This function is sufficiently efficient to compute the probability of interest up to quite large values of $n$. Basic tests show that it will compute for $n = 10^7$ in a few seconds. Below are some examples of use of this function for various values of $n$. As can be seen, even when the probability gets too small to see on the standard scale, we can still compute it on the log-probability scale.
p <- c(0.20, 0.25, 0.55);
PROB(1, p);
[1] 0.75
PROB(2, p);
[1] 0.6625
PROB(10, p);
[1] 0.4987803
PROB(100, p);
[1] 0.2503806
PROB(10^5, p);
[1] 1.174474e-123
PROB(10^6, p);
[1] 0
PROB(10^6, p, log = TRUE);
[1] -2795.471
PROB(10^7, p, log = TRUE);
[1] -27909.27
This function can be used to obtain the probability of interest up to large values of $n$. For values that are so large that they are computationally infeasible with this algorithm, we can then fall back on the normal approximation to the multinomial distribution.