Your time equations of motion look good. Most of what you have done is right, there are some subtleties to the meaning of the characteristic equation $ J_1J_2\omega^4 - k_t(J_1 + J_2)\omega^2 = 0$ (which you've gotten right) and you need to look at the physics and its relation to the initial conditions.
Assuming the real, trigonometric form is a headache, because you need to assume two different phase angles for the two different $\theta$s, i.e. you need to assume $\theta_i(t) = \hat{\theta}_i\cos(\omega\,t+\phi_i)$. So you don't get the cosine functions dropping out of your equations quite like you imagine. What you need is to assume functions of the form $\theta_i(t) = \hat{\theta}_i\exp(-i\,\omega\,t)$, so that the phase angle can be encoded in the phase of the complex scaling constants $\hat{\theta}_i$. When you use this form, you then find two of your natural frequencies, exactly as you have derived. There is indeed both a positive and negative natural frequency; both must be present so that we can build real valued trigonometric functions out of the complex exponentials: any sum of solutions is also a solution for these linear equations. But there are also two solutions $\omega^2 = 0$. Thus there must be a constant term in the solution as well (corresponding to the zero frequency solution). Actually, the repeated root means that a linear function of time will fulfill the equation (try solutions of the form $\theta_j(t)=A_j\,t+B_j$; these are solutions as long as $A_1=A_2$ and $B_1 = B_2$). Because we know the solution must be real for all $t$, the solutions must be of the form:
$$\theta_1(t) = \hat{\theta}_1\,\exp(-i\,\omega_0\,t) + \hat{\theta}_1^*\,\exp(+i\,\omega_0\,t) + A\,t+B$$
$$\theta_2(t) = \hat{\theta}_2\,\exp(-i\,\omega_0\,t) + \hat{\theta}_2^*\,\exp(+i\,\omega_0\,t) + A\,t+B$$
where $\omega_0 = \sqrt{\frac{k_t\,(J_1+J_2)}{J_1\,J_2}}$ as you found and $A\,B$ are real. Notice that $B$ corresponds to a constant angular offset of the system, and $A$ a constant angular speed of the system about its axis: you can take any general solution you get and set the system in uniform rotational motion on top of that solution and the total motion will still be a general solution.
Since $\dot{\theta}_j(0) = 0$ we get
$$i\,\omega_0(\hat{\theta}_j^*-\hat{\theta}_j) + A=0$$
so that:
$$-2\,\mathrm{Im}(\hat{\theta}_j)=-2\,\mathrm{Im}(\hat{\theta}_1)=-2\,\mathrm{Im}(\hat{\theta}_2)=A$$
To find the value of $A$, we take heed that the system begins with an angular momentum of nought ($\dot{\theta}_1(0)=\dot{\theta}_2(0)=0$). Angular momentum must be conserved, therefore at all times, we must have:
$$\begin{array}{lcl}J_1\,\dot{\theta}_1(t) + J_2\,\dot{\theta}_2(t) &=& J_1\,\left(i\,\omega_0(\hat{\theta}_1^*\,\exp(i\,\omega_0\,t)-\hat{\theta}_1\,\exp(-i\,\omega_0\,t) + A\right)+J_2\,\left(i\,\omega_0(\hat{\theta}_2^*\,\exp(i\,\omega_0\,t)-\hat{\theta}_2\,\exp(-i\,\omega_0\,t) + A\right)\\
&=&J_1\,\left(A-2\,\omega_0\,|\hat{\theta}_1|\,\sin(\omega_0\,t+\arg\hat{\theta}_1)\right)+J_2\,\left(A-2\,\omega_0\,|\hat{\theta}_2|\,\sin(\omega_0\,t+\arg\hat{\theta}_2)\right)\\
&=&0\;\forall\,t>0\end{array}$$
and so we see that
$$A=0$$
$$|\hat{\theta}_1|\,J_1+|\hat{\theta}_2|\,J_2=0$$
$$\arg\hat{\theta}_2 = \pi + \arg\hat{\theta}_1$$
and since we have already found that $-2\,\mathrm{Im}(\hat{\theta}_j)=A$ we now know $\arg\hat{\theta}_1 = 0;\,\arg\hat{\theta}_1 = \pi$ and so
$$\theta_1(t) = \alpha\,J_2\,\cos(\omega_0\,t) + B$$
$$\theta_2(t) = -\alpha\,J_1\,\cos(\omega_0\,t) + B$$
where it remains to find the common real scaling constant $\alpha$ and the offset $B$. From our initial conditions, we get from the above equations:
$$B=\frac{J_1\,\theta_1(0)+J_2\,\theta_2(0)}{J_1+J_2}$$
$$\alpha = \frac{\theta_1(0)-\theta_2(0)}{J_1+J_2}$$
Phew! We're here at last!
