The Duffing equation in its full form is
$$\ddot{x} + \delta \dot{x} -\alpha x + \beta x^3 = \gamma \cos(\omega t)$$
Now for specific values of the parameters several attractors exist (or not). Let's assume that $\alpha = \beta = \omega = 1$, $\delta = 0.15$, while $\gamma = 0.2445$. For these values of the parameters the system has two fixed-point attractors and a period-$3$ attractor. The period-$1$ attractors are located at about $(0.815, 0.242)$ and $(−0.933, 0.299)$. The period-$3$ attractor is located at about $(−1.412,−0.137)$, $(−0.354,−0.614)$, and $(0.645,−0.464)$.
My question is the following: How do we obtain the values of the positions of the attractors? Is there a theoretical way or is it done purely numerically? And if so, how?
The best analytical approach is in this case to use averaging techniques for the full nonlinear, weakly damped and weakly forced Duffing equation (i.e. $\beta$ is order 1, but both $\delta$ and $\gamma$ are reasonably small). The technique of averaging is well explained here for small $\beta$, which makes the situation a bit easier, as the 'basis functions' to average over are well-known trigonometric functions. If $\beta= 1$, those 'basis functions' or 'generating functions' for the unperturbed systems are Jacobi elliptic functions, see DLMF for more information. The averaging technique works exactly the same, only the integrals become a bit more cumbersome.
In I Kovacic and M J Brennan (eds), The Duffing Equation: Nonlinear Oscillators and their Behaviour, Wiley (2011), ISBN: 978-0-470-71549-9, this is mentioned on page 82. In particular, they give an extensive reference list to find results obtained using this technique and other techniques. To summarise: look at references [15-33] in Kovacic & Brennan, which are listed here.