Realization of complex high information density LCDs and systematic optimization of their electro-optical and ergonomic performance would not be possible in the required time-frame without reliable numerical modeling of the electro-optical performance of such display devices. In this paper we outline the history of numerical LDC modeling starting with Berreman and van Doorn, finally arriving at modern state-of-the-art LCD-modeling in two and three dimensions. Numerical modeling of LCDs is carried out in two steps: first, the effect of the electrical field on the orientation of the liquid crystalline alignment has to be evaluated before the corresponding optical properties can be computed. Starting from LC-elasticity theory we present suitable numerical methods for computing various states of LC-deformation (stable, metastable, bistable, etc.) in one- dimensional problems Light propagation in layered anisotropic absorbing media is evaluated with methods that are based on Maxwell's equations (Berreman 4 X 4-matrix approach). This approach can be simplified to yield methods with reduced computing time and sufficient accuracy for many problems (e.g. extended Jones 2 X 2-matrix formalism). A finite element method with automatic mesh generation and refinement for computing accurate solutions in two- dimensional problems is presented and its application illustrated with examples (e.g. IPS-effect, VAN-cells, etc.). In two- and three-dimensional problems, i.e. in cells with lateral dimensions comparable to the cell thickness, a variety of different director configurations are possible for a given geometry and electrical driving and addressing, making the modeling more complicated. Moreover, local defects can occur, which should also be considered in the simulation. Suitable approaches for the director field calculation, i.e. the vector and the tensor approach, are discussed. The complexity of the problem increases considerably when a third dimension is added, e.g. the geometry of the problem has to be defined in three dimensions together with the respective boundary conditions (anchoring geometry and elasticity) and electrodes. If strong deformations or even distortions are present in the orientation of the LC-layer, the applicability of known one- dimensional approaches for computing the optical properties must be checked and new approaches eventually have to be developed. The third dimension prohibits the use of some standard methods (e.g. FDTD), solely because of the enormous memory requirements and the long calculation times. Other approaches are presented and discussed.