If you solve the transient heat equation
$$\frac{\partial u}{\partial t} - k\nabla^2u = 0$$
for a very long time, and you then solve the Poisson equation
$$-k\nabla^2u = 0$$
with the same boundary conditions, you should get roughly the same result.
But we have to define what a "very long" time actually is.
You can write the solution of the time-dependent problem as
$$u(t) = \sum_ke^{-k\lambda_nt}\langle \phi_n, u(0)\rangle\phi_n + u(\infty)$$
where $\phi_n$, $\lambda_n$ are the eigenfunctions and eigenvalues of the Laplace operator on the domain, $u(0)$ is the initial value, and $u(\infty)$ is the steady-state value.
To 0th order, you only care about the smallest eigenmode because all others will decay much faster.
If $R$ is the radius of the largest circle that you can inscribe inside the domain, then
$$\lambda_1 \ge \pi^2 / R^2.$$
(See these notes, it's a consequence of the Cheeger inequality.)
So if $\pi^2 kt / R^2$ is greater than, say, 5 or 10 you should be pretty close to steady state.
Now, you mentioned the Jacobi, Gauss-Seidel, and SOR methods.
These iterative methods do converge to the actual solution of the linear system eventually.
But they are notoriously slow on their own, and they converge slowest for the smallest eigenmode $\phi_1$, which is the one that matters the most when you're looking at relaxation to steady state.
So if you find that your numerical solution for $u(t)$ is very different from $u(\infty)$ even for large $t$, it might actually be a consequence of using a slowly-converging iterative method.
In that case I'd look at using a sparse direct solver.