There is no nice, closed-form expression for the expected maximum of IID geometric random variables. However, the expected maximum of the corresponding IID exponential random variables turns out to be a very good approximation. More specifically, we have the hard bounds
$$\frac{1}{\lambda} H_n \leq E_n \leq 1 + \frac{1}{\lambda} H_n,$$
and the close approximation
$$E_n \approx \frac{1}{2} + \frac{1}{\lambda} H_n,$$
where $H_n$ is the $n$th harmonic number $H_n = \sum_{k=1}^n \frac{1}{k}$, and $\lambda = -\log (1-p)$, the parameter for the corresponding exponential distribution.
Here's the derivation. Let $q = 1-p$. Use Did's expression with the fact that if $X$ is geometric with parameter $p$ then $P(X \leq k) = 1-q^k$ to get
$$E_n = \sum_{k=0}^{\infty} (1 - (1-q^k)^n).$$
By viewing this infinite sum as right- and left-hand Riemann sum approximations of the corresponding integral we obtain
$$\int_0^{\infty} (1 - (1 - q^x)^n) dx \leq E_n \leq 1 + \int_0^{\infty} (1 - (1 - q^x)^n) dx.$$
The analysis now comes down to understanding the behavior of the integral. With the variable switch $u = 1 - q^x$ we have
$$\int_0^{\infty} (1 - (1 - q^x)^n) dx = -\frac{1}{\log q} \int_0^1 \frac{1 - u^n}{1-u} du = -\frac{1}{\log q} \int_0^1 \left(1 + u + \cdots + u^{n-1}\right) du $$
$$= -\frac{1}{\log q} \left(1 + \frac{1}{2} + \cdots + \frac{1}{n}\right) = -\frac{1}{\log q} H_n,$$
which is exactly the expression the OP has above for the expected maximum of $n$ corresponding IID exponential random variables, with $\lambda = - \log q$.
This proves the hard bounds, but what about the more precise approximation? The easiest way to see that is probably to use the Euler-Maclaurin summation formula for approximating a sum by an integral. Up to a first-order error term, it says exactly that
$$E_n = \sum_{k=0}^{\infty} (1 - (1-q^k)^n) \approx \int_0^{\infty} (1 - (1 - q^x)^n) dx + \frac{1}{2},$$
yielding the approximation
$$E_n \approx -\frac{1}{\log q} H_n + \frac{1}{2},$$
with error term given by
$$\int_0^{\infty} n (\log q) q^x (1 - q^x)^{n-1} \left(x - \lfloor x \rfloor - \frac{1}{2}\right) dx.$$
One can verify that this is quite small unless $n$ is also small or $q$ is extreme.
All of these results, including a more rigorous justification of the approximation, the OP's recursive formula, and the additional expression
$$E_n = \sum_{i=1}^n \binom{n}{i} (-1)^{i+1} \frac{1}{1-q^i},$$
are in Bennett Eisenberg's paper "On the expectation of the maximum of IID geometric random variables" (Statistics and Probability Letters 78 (2008) 135-143).