When one uses complex variables in this way one never multiplies two variables because the whole system is linear: if $z$ is the oscillating variable and you choose to represent it by a complex number, then things like $z\times z$ don't arise in a linear equation, so you don't get the kind of contradiction you astutely and clearly pointed out above. If the variable $z = r \exp(-i\,\omega\,t + i\,\phi)$ is multiplied by a complex coefficient $\alpha = a e^{i\,\theta}$ then this is done to impart (1) a scaling of amplitude by the real proportionality constant $a$ and (2) a phase shift of $\theta$ radians, i.e. a delay given by a fraction $\theta / (2\pi)$ of the oscillation period: there is no physical meaning to taking the real part before the scaling and delaying operation is imparted and the two complex numbers $\alpha$ and $z$ have been fully multiplied. So the procedure in this technique is that complex quantities are transformed to the corresponding real ones by taking the real part at the very end of the calculation, never before.
What does all this work? It's simply linearity: any multiple of the solution is also a solution and any sum of two solutions is also a solution. The general steady-state solution to a constant coefficient linear equation with forcing term ("input") $\cos(\omega\,t)$ is always a function of the form $a\cos(\omega\,t + \phi) = a \cos\phi \cos(\omega\,t) - a \sin\phi \sin(\omega\,t)$, so we can see that our "solution space" is a vector space spanned by the two functions $\cos(\omega\,t)$ and $\sin(\omega t)$. One is just simply changing basis and using the functions $\exp(i\,\omega\,t) = \cos(\omega\,t) + i\,\sin(\omega t)$ and $\exp(-i\,\omega\,t) = \cos(\omega\,t) - i\,\sin(\omega t)$ instead: by linearity is a perfectly valid. Moreover, if the linear system has real coefficients, then the solution corresponding to an input of the "negative frequency" wave $\exp(i\,\omega\,t)$ is simply the complex conjugate of the solution corresponding to the input of the "positive frequency" wave $\exp(-i\,\omega\,t)$. So we can deduce the system's behaviour in response to any forcing term of the form $a\,\cos(\omega\,t) + b\, \sin(\omega t)$ by simply deducing the behaviour in response to the prototypical $\exp(-i\,\omega\,t)$. We take the real part at the end of a calculation what we are doing in detail is this: we're averaging our solution with its complex conjugate, i.e. in one step we are inferring the system's behaviour for an input of the form $\exp(+i\,\omega\,t)$ and averaging with our solution, so we never really see the $\exp(+i\,\omega\,t)$ terms.
Why do we do this? It's easier. It's important to take heed that it's not essential. We could do everything in terms of sines and cosines, and the foregoing paragraph can be made to show that we get the same result by using complex ecponentials, and the kind of automatic "double handling" of positive and negative exponentials I described is likely the succintest description of how this convenience arises.
