the singularity theorems [ 178 , 104 , 105 , 106 ] state that einstein s equations will not evolve regular initial data arbitrarily far into the future or the past . an obstruction such as infinite curvature or the termination of geodesics will always arise to stop the evolution somewhere . the simplest , physically relevant solutions representing for example a homogeneous , isotropic universe ( friedmann - robertson - walker ( frw ) ) or a spherically symmetric black hole ( schwarzschild ) contain space - like infinite curvature singularities . although , in principle , the presence of a singularity could lead to unpredictable measurements for a physically realistic observer , this does not happen for these two solutions . the surface of last scattering of the cosmic microwave background in the cosmological case and the event horizon in the black hole ( bh ) case effectively hide the singularity from present day , external observers . singularity is generic , and the types of singularities that appear in generic spacetimes , remain major open questions in general relativity . the questions arise quickly , since other exact solutions to einstein s equations have singularities which are quite different from those described above . for example , the charged bh ( reissner - nordstrom solution ) has a time - like singularity . it also contains a cauchy horizon ( ch ) marking the boundary of predictability of space - like initial data supplied outside the bh . a test observer can pass through the ch to another region of the extended spacetime . the big bang in frw is classified as an asymptotically velocity term dominated ( avtd ) singularity [ 69 , 119 ] , since any spatial curvature term in the hamiltonian constraint becomes negligible compared to the square of the expansion rate as the singularity is approached . however , some anisotropic , homogeneous models exhibit mixmaster dynamics ( md ) [ 15 , 139 ] and are not avtd  the influence of the spatial scalar curvature can never be neglected . for more rigorous discussions of the classification and properties of the types of singularities . however , such methods yield either detailed knowledge of unrealistic , simplified ( usually by symmetries ) spacetimes or powerful , general results that do not contain details . to overcome these limitations , one might consider numerical methods to evolve realistic spacetimes to the point where the properties of the singularity may be identified . of course , most of the effort in numerical relativity applied to bh collisions has addressed the avoidance of singularities . numerical calculations , even more than analytic ones , require finite values for all quantities . ideally then , one must describe the singularity by the asymptotic non - singular approach to it . a numerical method which can follow the evolution into this asymptotic regime since the numerical study must begin with a particular set of initial data , the results can never have the force of mathematical proof . one may hope , however , that such studies will provide an understanding of the  phenomenology  of singularities that will eventually guide and motivate rigorous results . some examples of the interplay between analytic and numerical results and methods will be given here . in the following , we shall consider examples of numerical study of singularities both for asymptotically flat ( af ) spacetimes and for cosmological models . these examples have been chosen to illustrate primarily numerical studies whose focus is the nature of the singularity itself . in the af context , we shall consider two questions . one approach has been to explore highly non - spherical collapse looking for spindle or pancake singularities . if the formation of an event horizon requires a limit on the aspect ratio of the matter , such configurations may yield a naked singularity . recent analytic results suggest that one must go beyond the failure to observe an apparent horizon to conclude that a naked singularity has formed . another approach is to probe the limits between initial configurations which lead to black holes and those which yield no singularity at all ( i.e. flat spacetime plus radiation ) to explore the singularity as the bh mass goes to zero . this quest led naturally to the discovery of critical behavior in the collapse of a scalar field . the critical ( choptuik ) solution is a zero mass naked singularity ( visible from null infinity ) . however , it is a non - generic one since ( in addition to the fine - tuning required for this critical solution ) christodoulou has shown   that for the spherically symmetric einstein - scalar field equations , there always exists a perturbation that will convert a solution with a naked singularity to one with a black hole . reviews of critical phenomena in gravitational collapse can be found in [ 33 , 93 , 96 ] . the other question which is now beginning to yield to numerical attack involves the stability of the cauchy horizon ( ch ) in charged or rotating black holes . it has been conjectured [ 178 , 56 ] that a real observer , as opposed to a test mass , can not pass through the ch since realistic perturbed spacetimes will convert the ch to a strong spacelike singularity . numerical studies [ 40 , 67 , 47 ] show that a weak , null singularity forms first , as had been predicted [ 157 , 151 ] . in cosmology , we shall consider both the behavior of the mixmaster model and the issue of whether or not its properties are applicable to generic cosmological singularities . although numerical evolution of the mixmaster equations has a long history , developments in the past decade were motivated by inconsistencies between the known sensitivity to initial conditions and standard measures of the chaos usually associated with such behavior [ 145 , 164 , 166 , 19 , 76 , 45 , 111 , 160 ] . most recently , a coordinate invariant characterization of mixmaster chaos has been formulated   and a new extremely fast and accurate algorithm for mixmaster simulations developed . belinskii , khalatnikov , and lifshitz ( bkl ) long ago claimed [ 11 , 12 , 13 , 15 , 14 ] that it is possible to formulate the generic cosmological solution to einstein s equations near the singularity as a mixmaster universe at every point . while others have questioned the validity of this claim , it is only very recently that evidence for oscillatory behavior in the approach to the singularity of spatially inhomogeneous cosmologies has been obtained [ 181 , 31 ] . the weak cosmic censorship conjecture   requires a singularity formed from regular , asymptotically flat initial data to be hidden from an external observer by an event horizon . an excellent review of the meaning and status of weak cosmic censorship has been given by wald . counter examples have been known for a long time but tend to be dismissed as unrealistic in some way . the strong form of the cosmic censorship conjecture   forbids timelike singularities , even within black holes . perhaps the first numerical approach to study the cosmic censorship conjecture consisted of attempts to create naked singularities . many of these studies were motivated by thorne s  hoop conjecture    that collapse will yield a black hole only if a mass m is compressed to a region with circumference c  4m in all directions . ( as is discussed by wald , one runs into difficulties in any attempt to formulate the conjecture precisely . for example , how does one define c and m , especially if the initial data are not at least axially symmetric ? ) if the hoop conjecture is true , naked singularities may form if collapse can yield c  4m in some direction . the existence of a naked singularity is inferred from the absence of an apparent horizon ( ah ) which can be identified locally by following null geodesics . although a definitive identification of a naked singularity requires the event horizon ( eh ) to be proven to be absent , to identify an eh requires knowledge of the entire spacetime . while one finds an ah within an eh [ 120 , 121 ] , it is possible to construct a spacetime slicing which avoids the ah even though an eh is present . methods to find an eh in a numerically determined spacetime have only recently become available and have not been applied to this issue [ 132 , 136 ] . in the best known attempt to produce naked singularities , shapiro and teukolsky ( st )   considered collapse of prolate spheroids of collisionless gas . ( nakamura and sato   had previously studied the collapse of non - rotating deformed stars with an initial large reduction of internal energy and apparently found spindle or pancake singularities in extreme cases . ) st solved the general relativistic vlasov equation for the particles along with einstein s equations for the gravitational field . they found that an ah ( and presumably a bh ) formed if c  4m < 1 everywhere but no ah ( and presumably a naked singularity ) in the opposite case . in the latter case , the evolution ( not surprisingly ) could not proceed past the moment of formation of the singularity . in a subsequent study , st   also showed that a small amount of rotation ( counter rotating particles with no net angular momentum ) does not prevent the formation of a naked spindle singularity . however , wald and iyer   have shown that the schwarzschild solution has a time slicing whose evolution approaches arbitrarily close to the singularity with no ah in any slice , ( but of course , with an eh in the spacetime ) . this may mean that there is a chance that the increasing prolateness found by st in effect changes the slicing to one with no apparent horizon just at the point required by the hoop conjecture . while , on the face of it , this seems unlikely , tod gives an example where the ah does not form on a chosen constant time slice  but rather different portions form at different times . he argues that a numerical simulation might be forced by the singularity to end before the formation of the ah is complete . such an ah would not be found by the simulations . in response to such a possibility , shapiro and teukolsky considered equilibrium sequences of prolate relativistic star clusters . the idea is to counter the possibility that an eh might form after the simulation must stop . if an equilibrium configuration is non - singular , it can not contain an eh , since singularity theorems say that an eh implies a singularity . however , a sequence of nonsingular equilibria with rising i ever closer to the spindle singularity would lend support to the existence of a naked spindle singularity since one can approach the singular state without formation of an eh . they constructed this sequence and found that the singular end points were very similar to their dynamical spindle singularity . wald believes , however , that it is likely that the st slicing is such that their singularities are not naked  the ah is present but has not yet appeared in their time slices . another numerical study of the hoop conjecture was made by chiba et al . rather than a dynamical collapse model , they searched for ah s in analytic initial data for discs , annuli , and rings . previous studies of this type were done by nakamura et al   with oblate and prolate spheroids and by wojtkiewicz   with axisymmetric singular lines and rings . the summary of their results is that an ah forms if c  4m  1.26 . ( analytic results due to barrabes et al [ 5 , 4 ] and tod   give similar quantitative results with different initial data classes and ( possibly ) definition of c. ) there is strong analytical evidence against the development of spindle singularities . it has been shown by chrusciel and moncrief that strong cosmic censorship holds in af electrovac solutions which admit a g2 symmetric cauchy surface . the evolutions of these highly nonlinear equations are in fact non - singular . motivated by st s results , echeverria   numerically studied the properties of the naked singularity that is known to form in the collapse of an infinite , cylindrical dust shell . while the asymptotic state can be found analytically , the approach to it must be followed numerically . the analytic asymptotic solution can be matched to the numerical one ( which can not be followed all the way to the collapse ) to show that the singularity is strong ( an observer experiences infinite stretching parallel to the symmetry axis and squeezing perpendicular to the symmetry axis ) . a burst of gravitational radiation emitted just prior to the formation of the singularity stretches and squeezes in opposite directions to the singularity . one wonders then if dust collapse gives any information about singularities of the gravitational field . nakamura et al ( nsn )   conjectured that even if naked spindle singularities could exist , they would either disappear or become black holes . this demise of the naked singularity would be caused by the back reaction of the gravitational waves emitted by it . while nsn proposed a numerical test of their conjecture , they believed it to be beyond the current generation of computer technology . the results of all these searches for naked spindle singularities are controversial , but they could be resolved if the presence or absence of the eh could be determined . one could demonstrate numerically whether or not wald s view of st s results is correct by using existing eh finders [ 132 , 136 ] in a relevant simulation . we now consider an effect that was originally found by choptuik   in a numerical study of the collapse of a spherically symmetric massless scalar field . for recent reviews see [ 93 , 96 ] . we note that this is the first completely new phenomenon in general relativity to be discovered by numerical simulation . in collapse of a scalar field , essentially two things can happen : either a black hole ( bh ) forms , or the scalar waves pass through each other and disperse . choptuik discovered that for any 1-parameter set of initial data labeled by p , there is a critical value p * such that p > p * yields a bh . he found 1\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${m_{bh } } \approx { c_f}{(p - p*)^\gamma}$$ 
 \end{document}mbhcf(pp*) where mbh is the mass of the eventual bh . the constant cf depends on the parameter of the initial data that is selected , but   .37 is the same for all choices . = ln(t0 *  t0 ) +  ( t0 is the proper time of an observer at r = 0 with t0 * the finite proper time at which the critical evolution concludes and  is a constant which scales r ) , the waveform x repeats ( echoes ) at intervals of  in  if  is rescaled to    , i.e. x(   , the scaling behavior ( 1 ) demonstrates that the minimum bh mass ( for bosons ) is zero . the critical solution itself is a counter - example to cosmic censorship , ( since the formation of the zero mass bh causes high curvature regions to become visible at r =  ) . the numerical demonstration of this feature of the critical solution was provided by hamad and stewart . this result caused hawking to pay off a bet to preskill and thorne [ 44 , 123 ] . soon after this discovery , scaling and critical phenomena were found in a variety of contexts . abrahams and evans   discovered the same phenomenon in axisymmetric gravitational wave collapse with a different value of  and , to within numerical error , the same value of . ( note that the rescaling of r with e  30 required choptuik to use adaptive mesh refinement ( amr ) to distinguish subsequent echoes . abrahams and evans smaller  ( e  1.8 ) allowed them to see echoing with their 2 + 1 code without amr . ) garfinkle   confirmed choptuik s results with a completely different algorithm that does not require amr . he used goldwirth and piran s   method of simulating christodoulou s   formulation of the spherically symmetric scalar field in null coordinates . this formulation allowed the grid to be automatically rescaled by choosing the edge of the grid to be the null ray that just hits the central observer at the end of the critical evolution . ( missing points of null rays that cross the central observer s world line are replaced by interpolation between those that remain . ) hamad and stewart   have also repeated choptuik s calculation using null coordinates and amr . evans and coleman   realized that self - similar rather than self - periodic collapse might be more tractable both numerically ( since ode s rather than pde s are involved ) and analytically . however , most of the action occurs in the center so that a match of the self - similar inner region to an outer af one should always be possible . ) in a series of papers , hirschmann and eardley [ 108 , 109 ] developed a ( numerical ) self - similar solution to the spherically symmetric complex scalar field equations . these are ode s with too many boundary conditions causing a solution to exist only for certain fixed values of . numerical solution of this eigenvalue problem allows very accurate determination of . the self - similarity also allows accurate calculation of  as follows : the critical p = p * solution is unstable to a small change in p. at any time t ( where t < 0 is increasing toward zero ) , the amplitude a of the perturbation exhibits power law growth : 2\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$a \propto ( p - p*){(- t)^{- \kappa}}$$ 
 \end{document}a(pp*)(t) where  > 0 . at any fixed t , larger a implies larger mbh . equivalently , any fixed amplitude a =  will be reached faster for larger eventual mbh . scaling arguments give the dependence of mbh on the time at which any fixed amplitude is reached : 3\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${m_{bh } } \propto ( - { t_1})$$ 
 \end{document}mbh(t1 ) where 4\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$(p - p*){(- { t_1})^\kappa } \propto \delta.$$ 
 \end{document}(pp*)(t1). thus 5\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${m_{bh } } \propto { ( p - { p^*})^{1/\kappa}}.$$ 
 \end{document}mbh(pp*)1/. therefore , one need only identify the growth rate of the unstable mode to obtain an accurate value of  = 1/. it is not necessary to undertake the entire dynamical evolution or probe the space of initial data . hirschmann and eardley obtain  = 0.387106 for the complex scalar field solution , while koike et al   obtain  = 0.35580192 for the evans - coleman solution . although the similarities among the critical exponents  in the collapse computations suggested a universal value , maison   used these same scaling - perturbation methods to show that  depends on the equation of state p = k of the fluid in the evans - coleman solution . gundlach reformulates the model as a nonlinear hyperbolic boundary value problem with eigenvalue  and finds  = 3.4439 . as with the self - similar solutions described above , the critical solution is found directly without the need to perform a dynamical evolution or explore the space of initial data . hara et al extended the renormalization group approach of   to the discretely - self - similar case . recently , gundlach   completed his eigenvalue analysis of the choptuik solution to find the growth rate of the unstable mode to be  = .374  .001 . self - similar critical behavior has been seen in string theory related axion - dilaton models [ 68 , 101 ] and in the nonlinear -model . garfinkle and duncan have shown that subcritical collapse of a spherically symmetric scalar field yields a scaling relation for the maximum curvature observed by the central observer with critical parameters that would be expected on the basis of those found for supercritical collapse . choptuik et al   have generalized the original einstein - scalar field calculations to the einstein - yang - mills ( eym ) ( for su(2 ) ) case . here the critical solution is a regular , unstable solution to the eym equations found previously by bartnik and mckinnon . in type ii collapse , the minimum bh mass is zero with the critical solution similar to that of choptuik ( with a different   0.20 ,   0 . the type i behavior arises when the collapsing system has a metastable static solution in addition to the choptuik critical one . brady , chambers and goncalves [ 54 , 37 ] conjectured that addition of a mass to the scalar field of the original model would break scale invariance and might yield a distinct critical behavior . phases  seen in the eym case . the type ii solution can be understood as perturbations of choptuik s original model with a small scalar field mass . here small means that   1 where  is the spatial extent of the original nonzero field region . ( the scalar field is well within the compton wavelength corresponding to . ) on the other hand ,   1 yields type i behavior . the minimum mass critical solution is an unstable soliton of the type found by seidel and suen . the massive scalar field can be treated as collapsing dust to yield a criterion for bh formation . the choptuik solution has also been found to be the critical solution for charged scalar fields [ 98 , 113 ] . as p  p * , q / m  0 for the black hole . numerical study of the critical collapse of collisionless matter ( einstein - vlasov equations ) has yielded a non - zero minimum bh mass . an astrophysical application of bh critical phenomena has been considered by nimeyer and jedamzik   and yokoyama . they consider its implications for primordial bh formation and suggest that it could be important . the question is then why these critical phenomena should appear in so many collapsing gravitational systems . the discrete self - similarity of choptuik s solution may be understood as scaling of a limit cycle . ( the self - similarity of other systems may be understood as scaling of a limit point . ) garfinkle   has conjectured that the scale invariance of einstein s equations might provide an underlying explanation for the self - similarity and discrete - self - similarity found in collapse . until recently , he considered spherical and non - spherical perturbations of \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$p = { 1 \over 3}\rho$\end{document}p=13 perfect fluid collapse . only the original ( spherical ) growing mode survived . unlike the simple singularity structure of the schwarzschild solution , where the event horizon encloses a spacelike singularity at r = 0 , charged and/or rotating bh s have a much richer singularity structure . the extended spacetimes have an inner cauchy horizon ( ch ) which is the boundary of predictability . poisson and israel [ 157 , 158 ] began an analytic study of the effect of perturbations on the ch . their goal was to check conjectures that the blue - shifted infalling radiation during collapse would convert the ch into a true singularity and thus prevent an observer s passage into the rest of the extended regions . by including both ingoing and back - scattered outgoing radiation , they find for the reissner - nordstrm ( rn ) solution that the mass function ( qualitatively r however , ori showed both for rn and kerr [ 151 , 152 ] that the metric perturbations are finite ( even though r 
 r diverges ) so that an observer would not be destroyed by tidal forces ( the tidal distortion would be finite ) and could survive passage through the ch . a numerical solution of the einstein - maxwell - scalar field equations could test these perturbative results . gnedin and gnedin   have numerically evolved the spherically symmetric einstein - maxwell with massless scalar field equations in a 2 + 2 formulation . the initial conditions place a scalar field on part of the rn event horizon ( with zero field on the rest ) . an asymptotically null or spacelike singularity whose shape depends on the strength of the initial perturbation replaces the ch . for a sufficiently strong perturbation , although they claim to have found that the ch evolved to become a spacelike singularity , the diagrams in their paper show at least part of the final singularity to be null or asymptotically null in most cases . brady and smith   used the goldwirth - piran formulation   to study the same problem . they follow the evolution of the ch into a null singularity , demonstrate mass inflation , and support ( with observed exponential decay of the metric component g ) the validity of previous analytic results [ 157 , 158 , 151 , 152 ] including the  weak  nature of the singularity that forms . they find that the observer hits the null ch singularity before falling into the curvature singularity at r = 0 . whether or not these results are in conflict with gnedin and gnedin   is unclear . the collapse of a scalar field in a charged , spherically symmetric spacetime causes an initial rn ch to become a null singularity , except at r = 0 , where it is spacelike . the observer falling into the bh experiences ( and potentially survives ) the weak , null singularity [ 151 , 152 , 36 ] before the spacelike singularity forms . this has been confirmed by droz   using a plane wave model of the interior and by burko   using a collapsing scalar field . recently , numerical studies of the interiors of non - abelian black holes have been carried out by breitenlohner et al [ 41 , 42 ] and by galtsov et al [ 66 , 78 , 80 , 79 ] ( see also ) . although there appear to be conflicts between the two groups , the differences can be attributed to misunderstandings of each other s notation . the main results include an interesting oscilliatory behavior of the metric . replacing the af boundary conditions with schwarzschild - de sitter and rn - de sitter bh s was long believed to yield a counterexample to strong cosmic censorship . ( see [ 137 , 138 , 156 , 53 ] and references therein for background and extended discussions . ) the stability of the ch is related to the decay tails of the radiating scalar field . numerical studies recently determined these to be exponential [ 38 , 53 , 55 ] rather than power law as in af spacetimes . the decay tails of the radiation are necessary initial data for numerical study of ch stability   and are crucial to the development of the null singularity . analytic studies had indicated that the ch is stable in rn - de sitter bh s for a restricted range of parameters near extremality . however , brady et al   have shown ( using linear perturbation theory ) that , if one includes the backscattering of outgoing modes which originate near the event horizon , the ch is always unstable for all ranges of parameters numerical studies are needed to demonstrate the existence of a null singularity at the ch in nonlinear evolution . as a potentially useful approach to the numerical study of singularities , we consider hbner s [ 115 , 116 , 117 ] numerical scheme to evolve on a conformal compactified grid using friedrich s formalism . so far , this code has been used to locate singularities and to identify choptuik s scaling . the first , called asymptotically velocity term dominated ( avtd ) [ 69 , 119 ] , refers to a cosmology that approaches the kasner ( vacuum , bianchi i ) solution   as   . ( spatially homogeneous universes can be described as a sequence of homogeneous spaces labeled by . here we shall choose  so that  =  coincides with the singularity . ) an example of such a solution is the vacuum bianchi ii model   which begins with a fixed set of kasner - like anisotropic expansion rates , and , possibly , makes one change of the rates in a prescribed way ( mixmaster - like bounce ) and then continues to  =  as a fixed kasner solution . in contrast are the homogeneous cosmologies , which display mixmaster dynamics such as vacuum bianchi viii and ix [ 15 , 139 , 100 ] and bianchi vi0 and bianchi i with a magnetic field [ 131 , 20 , 130 ] . mixmaster dynamics describes an approach to the singularity which is a sequence of kasner epochs with a prescription , originally due to belinskii , khalatnikov , and lifshitz ( bkl ) , for relating one kasner epoch to the next . some of the mixmaster bounces ( era changes ) display sensitivity to initial conditions one usually associates with chaos , and in fact mixmaster dynamics is chaotic . the vacuum bianchi i ( kasner ) solution is distinguished from the other bianchi types in that the spatial scalar curvature r , ( proportional to ) the minisuperspace ( mss ) potential [ 139 , 167 ] , vanishes identically . but r arises in other bianchi types due to spatial dependence of the metric in a coordinate basis . thus an avtd singularity is also characterized as a regime in which terms containing or arising from spatial derivatives no longer influence the dynamics . this means that the mixmaster models do not have an avtd singularity , since the influence of the spatial derivatives ( through the mss potential ) never disappears  there is no last bounce . symplectic numerical methods have proven useful in studies of the approach to the singularity in cosmological models . symplectic ode and pde [ 75 , 141 ] methods comprise a type of operator splitting . details of the application to the gowdy model ( pde s in one space and one time direction ) are given elsewhere . for a field q(t ) and its conjugate momentum p(t ) split the hamiltonian operator into kinetic and potential energy subhamiltonians . thus , 6\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$h = \int\nolimits { dx } \{{1 \over 2}{p^2 } + v(q)\ } = { h_1}(p ) + { h_2}(q).$$ 
 \end{document}h=dx{12p2+v(q)}=h1(p)+h2(q ) . if the vector x = ( p , q ) defines the variables at time t , then the time evolution is given by 7\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${{dx } \over { dt } } = { \{h , x\ } _ { pb } } \equiv ax$$ 
 \end{document}dxdt={h , x}pbax where { } pb is the poisson bracket . the usual exponentiation yields an evolution operator 8\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${e^{a\delta t } } = { e^{{a_1}}}^{(\delta t/2)}{e^{{a_2}\delta t}}{e^{{a_1}}}^{(\delta t/2 ) } + o(\delta { t^3})$$ 
 \end{document}eat = ea1(t/2)ea2tea1(t/2)+o(t3 ) for a = a1 + a2 the generator of the time evolution . higher order accuracy may be obtained by a better approximation to the evolution operator [ 172 , 173 ] . this method is useful when exact solutions for the subhamiltonians are known . for the given h , variation of h1 yields the solution 9\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } \begin{document}$$q = { q_0 } + { p_0}\,\delta t\quad , \quad p = { p_0},$$ 
 \end{document}q = q0+p0t , p = p0 , while that of h2 yields 10\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$q = { q_0}\quad , \quad p = { p_0 } - { \left . { { { dv } \over { dq } } } \right|_{{q_0}}}\delta t \ , \ , .$$ 
 \end{document}q = q0,p = p0dvdq|q0t . note that h2 is exactly solvable for any potential v no matter how complicated , although the required differenced form of the potential gradient may be non - trivial . one evolves from t to t + t using the exact solutions to the sub - hamiltonians according to the prescription given by the approximate evolution operator ( 8) . extension to more degrees of freedom and to fields is straightforward [ 30 , 21 ] . belinskii , khalatnikov , and lifshitz   ( bkl ) described the singularity approach of vacuum bianchi ix cosmologies as an infinite sequence of kasner   epochs whose indices change when the scalar curvature terms in einstein s equations become important . they were able to describe the dynamics approximately by a map evolving a discrete set of parameters from one kasner epoch to the next [ 15 , 57 ] . for example , the kasner indices for the power law dependence of the anisotropic scale factors can be parametrized by a single variable u  1 . bkl determined that 11\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${u_{n + 1 } } = \left\{{\begin{array}{*{20}c } { { u_n } - 1\quad , \quad 2 \le { u_n}}\quad \quad \quad \\ { { { ( { u_n } - 1)}^{- 1}}\quad , \quad 1 \le { u_n } \le 2 } \\ the subtraction in the denominator for 1  un  2 yields the sensitivity to initial conditions associated with mixmaster dynamics ( md ) . misner   described the same behavior in terms of the model s volume and anisotropic shears . a multiple of the scalar curvature acts as an outward moving potential in the anisotropy plane . kasner epochs become straight line trajectories moving outward along a potential corner while bouncing from one side to the other . a change of corner ends a bkl era when u  ( u  1 ) . numerical evolution of einstein s equations was used to explore the accuracy of the bkl map as a descriptor of the dynamics as well as the implications of the map [ 145 , 164 , 166 , 19 ] . rendall has studied analytically the validity of the bkl map as an approximation to the true trajectories . later , the bkl sensitivity to initial conditions was discussed in the language of chaos [ 6 , 125 ] . an extended application of bernoulli shifts and farey trees was given by rugh   and repeated by cornish and levin . however , the chaotic nature of mixmaster dynamics was questioned when numerical evolution of the mixmaster equations yielded zero lyapunov exponents ( le s ) [ 76 , 45 , 111 ] . only an exponential divergence , characteristic of a chaotic system , will yield positive exponent . ) this issue was resolved when the le was shown numerically and analytically to depend on the choice of time variable [ 164 , 18 , 73 ] . although md itself is well - understood , its characterization as chaotic or not had been quite controversial . leblanc et al   have shown ( analytically and numerically ) that md can arise in bianchi vi0 models with magnetic fields ( see also ) . in essence , the magnetic field provides the wall needed to close the potential in a way that yields the bkl map for u . a similar study of magnetic bianchi i has been given by leblanc . recently , cornish and levin ( cl )   identified a coordinate invarient way to characterize md . sensitivity to initial conditions can lead to qualitatively distinct outcomes from initially nearby trajectories . while the le measures the exponential divergence of the trajectories , one can also  color code  the regions of initial data space corresponding to particular outcomes . cl defined the following set of discrete outcomes : during a numerical evolution , the bkl parameter u is evaluated from the trajectories . the first time u > 7 appears in an approximately kasner epoch , the trajectory is examined to see which metric scale factor has the largest time derivative . this defines three outcomes and thus three colors for initial data space . to study the cl fractal and ergodic properties of mixmaster evolution , one could take advantage of a new numerical algorithm due to berger , garfinkle , and strasser . symplectic methods are used to allow the known exact solution for a single mixmaster bounce   to be used in the ode solver . standard ode solvers   can take large time steps in the kasner segments but must slow down at the bounce . tens of orders of magnitude improvement in speed are obtained while the accuracy of machine precision is maintained . in , this required an improvement of the bkl map ( for parameters other than u ) to take into account details of the exponential potential . one can easily invent prescriptions other than that given by cornish and levin   which would yield discrete outcomes . it is not clear how the value of the fractal dimension as measured by cornish and levin would be affected . for example , carretero - gonzalez et al   find evidence of chaotic behavior in bianchi ix - brans - dicke solutions while cotsakis et al   have shown that bianchi ix models in 4th order gravity theories have stable non - chaotic solutions . barrow and levin find evidence of chaos in bianchi ix einstein - yang - mills ( eym ) cosmologies . finally , we remark on a successful application of numerical regge calculus in 3 + 1 dimensions . bkl have conjectured that one should expect a generic singularity in spatially inhomogeneous cosmologies to be locally of the mixmaster type . for a review of homogeneous cosmologies , inhomogeneous cosmologies , and the relation between them , see . the main difficulty with the acceptance of this conjecture has been the controversy over whether the required time slicing can be constructed globally [ 8 , 89 ] . montani , belinskii , and kirillov and kochnev [ 128 , 127 ] have pointed out that if the bkl conjecture is correct , the spatial structure of the singularity could become extremely complicated as bounces occur at different locations at different times . a class of cosmological models which might have local md are vacuum universes on t  r with a u(1 ) symmetry . the non - commuting killing vectors of local md can be constructed since only one killing vector is already present . the two commuting killing vectors of the even simpler plane symmetric gowdy cosmologies [ 88 , 17 ] preclude their use to test the conjecture . however , these models are interesting in their own right since they have been conjectured to possess an avtd singularity . polarized plane symmetric cosmologies have been evolved numerically using standard techniques by anninos , centrella , and matzner [ 2 , 3 ] . a long - term project involving berger , garfinkle , and moncrief and our students to study the generic cosmological singularity numerically has been reviewed in detail elsewhere   but will be discussed briefly here . the gowdy model on t  r is described by gravitational wave amplitudes p( ,  ) and q( ,  ) which propagate in a spatially inhomogeneous background universe described by ( ,  ) . ( we note that the physical behavior of a gowdy spacetime can be computed from the effect of the metric evolution on a test cylinder . ) the time variable  measures the area in the symmetry plane with  =  a curvature singularity . the first is nonlinearly coupled wave equations for dynamical variables p and q ( where , a = /a ) obtained from the variation of   12\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\begin{array}{*{20}c } { h \quad = \quad { 1 \over 2}\int\limits_0^{2\pi } { d\theta } [ \pi _ p^2 + { e^{- 2p}}\pi _ q^2]\quad \quad \quad \quad \quad \quad \quad \quad } \\ { + \quad { 1 \over 2}\int\limits_0^{2\pi } { d\theta } [ { e^{- 2\tau}}(p,_\theta ^2 + { e^{2p}}q,_{\theta}^2 ) ] = { h_1 } + { h_2}. } \\ if the model is , in fact , avtd , the approximation in the symplectic numerical scheme should become more accurate as the singularity is approached . these can be expressed as first order equations for  in terms of p and q. this break into dynamical and constraint equations removes two of the most problematical areas of numerical relativity from this model  the initial value problem and numerical preservation of the constraints . for the special case of the polarized gowdy model ( q = 0 ) , p satisfies a linear wave equation whose exact solution is well - known . avtd behavior is defined in   as follows : solve the gowdy wave equations neglecting all terms containing spatial derivatives . if the approach to the singularity is avtd , the full solution comes arbitrarily close to an avtd solution at each spatial point as   . as    , the avtd solution becomes 13\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$p(\theta , \ , \tau ) \to v ( \theta)\tau \quad , \quad q(\theta , \tau ) \to { q_0}(\theta)$$ 
 \end{document}p(,)v(),q(,)q0( ) where  > 0 . substitution in the wave equations shows that this behavior is consistent with asymptotic exponential decay of all terms containing spatial derivatives only if 0    1 . we have shown that , except at isolated spatial points , the nonlinear terms in the wave equation for p drive  into this range [ 25 , 27 ] . the exceptional points occur when coefficients of the nonlinear terms vanish and are responsible for the growth of spiky features seen in the wave forms [ 30 , 25 ] . we conclude that generic gowdy cosmologies have an avtd singularity except at isolated spatial points [ 25 , 27 ] . recently , it has been proven analytically that gowdy solutions with 0 <  < 1 and avtd behavior almost everywhere are generic . addition of a magnetic field to the vacuum gowdy models ( plus a topology change ) which yields the inhomogeneous generalization of magnetic bianchi vi0 models provides an additional potential which grows exponentially if 0 <  < 1 . local mixmaster behavior has recently been observed in these magnetic gowdy models . given the success of the symplectic method in studying the singularity behavior of the gowdy model , we can consider its extension to the case of u(1 ) symmetric cosmologies . moncrief has shown   that cosmological models on t  r with a spatial u(1 ) symmetry can be described by five degrees of freedom { x , z ,  ,  ,  } and their respective conjugate momenta { px , pz , p , p , r } all variables are functions of spatial variables u , v and time , . einstein s equations can be obtained by variation of 14\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\begin{array}{*{20}c } { h\quad = \quad \oint { \oint { du\,dv\,{\mathcal h } } } \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \,\ , } \\ { \quad = \quad \oint { \oint { du\,dv\,\left({{1 \over 8}p_z^2 + { 1 \over 2}{e^{4z}}p_x^2 + { 1 \over 8}{p^2 } + { 1 \over 2}{e^{4\varphi}}{r^2 } - { 1 \over 2}p_\lambda ^2 + 2{p_\lambda } } \right ) } } } \\ { + { e^{- 2\tau } } = \oint { \oint { du\,dv } \left\{{\left({{e^\lambda}{e^{ab } } } \right){,_{ab } } - \left({{e^\lambda}{e^{ab } } } \right){,_a}\lambda { , _ b}\quad \quad } \right . } \quad } \\ { + { e^\lambda}\left [ { \left({{e^{- 2z } } } \right){,_u}x{,_v } - \left({{e^{- 2z } } } \right){,_v}x{,_u } } \right]\quad \quad \quad \quad \quad \quad \quad \quad } \\ { \left . { + 2{e^\lambda}{e^{ab}}\varphi { , _ a}\varphi { , _ b } + { 1 \over 2}{e^\lambda}{e^{- 4\varphi}}{e^{ab}}\omega { , _ a}\omega { , _ b } } \right\}\quad \quad \quad \quad \quad \quad } \\ { = \quad { h_1 } + { h_2}.\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad } \\ \end{array}$$ 
 \end{document}h=dudv=dudv(18pz2 + 12e4zpx2 + 18p2 + 12e4r212p2 + 2p)+e2=dudv{(eeab),ab(eeab),a,b+e[(e2z),ux , v(e2z),vx , u]+2eeab,a,b+12ee4eab,a,b}=h1+h2 . here and  are analogous to p and q while e is a conformal factor for the metric eab(x , z ) in the u - v plane perpendicular to the symmetry direction . note particularly that h1 contains two copies of the gowdy h1 plus a free particle term and is thus exactly solvable . however , it still contains no momenta , so its equations are trivially exactly solvable . a spectral evaluation of derivatives   which has been shown to work in gowdy simulations   does not appear to be helpful in the u(1 ) case . ) a particular solution to the initial value problem is used , since the general solution is not available . currently , except as discussed below , the constraints are monitored but not explicitly enforced . despite the current limitations of the u(1 ) code , conclusions can be confidently drawn for polarized u(1 ) models in our restricted class of initial data . methods similar to those used in , indicate that the term containing gradients of  in ( 14 ) acts as a mixmaster - like potential to drive the system away from avtd behavior in generic u(1 ) models . numerical simulations provide support for this suggestion [ 27 , 31 ] . whether this potential term grows or decays depends on a function of the field momenta . however , failure to enforce the constraints can cause an erroneous relationship among the momenta to yield qualitatively wrong behavior . there is numerical evidence that this error tends to suppress mixmaster - like behavior leading to apparent avtd behavior in extended spatial regions [ 21 , 22 ] . in fact , it has been found very recently   that when the hamiltonian constraint is enforced at every time step , the predicted local oscillatory behavior of the approach to the singularity is observed . mixmaster simulations with the new algorithm   can easily evolve more than 250 bounces reaching ||  10 . the larger number of bounces quickly reveals that it is necessary to enforce the hamiltonian constraint . a similar conclusion for gravitational wave equations although this particular analysis may be incorrect , it is still likely that constraint enforcement will be essential for sufficiently long simulations . an explicitly constraint enforcing u(1 ) code was developed some years ago by ove ( see   and references therein ) . numerical investigation of singularities provides an arena for the close coupling of analytic and numerical techniques . the references provided here contain many examples of analytic results guided by numerical results and numerical studies to demonstrate whether or not analytic conjectures are valid . while the search for naked singularities in the collapse of highly prolate systems has yielded controversial results , a naked singularity was discovered in the collapse of spherically symmetric scalar fields . the numerical exploration of cosmological singularities has yielded strong evidence that the asymptotic behavior is local  each spatial point evolves toward the singularity as a separate universe . this area of research then should begin to overlap with the studies of black hole interiors ( see for example ) .
this review updates a previous review article [ 22 ] 
 . numerical exploration of the properties of singularities could , in principle , yield detailed understanding of their nature in physically realistic cases . 
 examples of numerical investigations into the formation of naked singularities , critical behavior in collapse , passage through the cauchy horizon , chaos of the mixmaster singularity , and singularities in spatially inhomogeneous cosmologies are discussed .