Nowadays the rotation of the Earth can be observed with an accuracy of about 0.01 milliarcseconds (mas ), while theoretical models are able to describe this motion at a level of 1 mas. This mismatch is partly due to the enormous complexity of the involved processes, operating on different time scales and driven by a large variety of physical effects. But al so partly due to the used models, which often use simplified and linearized equations to obtain the solution analytically. In this work we present our new numerical theory of the rotation of the Earth. The model underlying the theory is fully compatible with the post - Newtonian approximation of general relativity and is formulated using ordinary differential equations for the angles describing the orientation of the Earth (or its particular layers) in the GCRS. These equations are then solved numerically to describe the rotational motion with highest accuracy. Being initially developed for a rigid Earth our theory was extended towards a more realistic Earth model. In particular, we included 3 different layers (crust, fluid outer core and solid inner core) and all important coupling torques between them as well as all important effects of non - rigidity, such as elastic deformation, relative angular momenta due to atmosphere and ocean etc. In our presentation we will describe the details of our work and compare i t to the currently used models of Earth rotation. Further, we discuss possible applications of our numerical theory to obtain high - accuracy models of rotational motion of other celestial bodies such as Mercury.