Numerical methods

3.2 Numerical methods

3.2.1 Symplectic methods

Symplectic numerical methods have proven useful in studies of the approach to the singularity in cosmological models [38 ]. Symplectic ODE and PDE [101 , 189 ] methods comprise a type of operator splitting. An outline of the method (for one degree of freedom) follows. Details of the application to the Gowdy model (PDE’s in one space and one time direction) are given elsewhere [42 ].

For a field $q(t)$ and its conjugate momentum $p(t)$ , the Hamiltonian operator splits into kinetic and potential energy sub-Hamiltonians. Thus, for an arbitrary potential $V (q)$ ,

If the vector $X = (p,q)$ defines the variables at time $t$ , then the time evolution is given by

where ${} P B$ is the Poisson bracket. The usual exponentiation yields an evolution operator

for $A = A1 + A2$ being the generator of the time evolution. Higher order accuracy may be obtained by a better approximation to the evolution operator [234 , 235 ]. This method is useful when exact solutions for the sub-Hamiltonians are known. For the given $H$ , variation of $H1$ yields the solution

while that of $H2$ yields

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 [42

, 27

3.2.2 Other methods

Symplectic methods can achieve convergence far beyond that of their formal accuracy if the full solution is very close to the exact solution from one of the sub-Hamiltonians. Examples where this is the case are given in [39 , 36 ]. On the other hand, because symplectic algorithms are a type of operator splitting, suboperators can be subject to instabilities that are suppressed by the full operator. An example of this may be found in [41 ]. Other types of PDE solvers are more effective for such spacetimes. One currently popular method is iterative Crank–Nicholson (see [237 ]) which is, in effect, an implicit method without matrix inversion. It was first applied to numerical cosmology by Garfinkle [111 ] and was recently used in this context to evolve $T2$ symmetric cosmologies [41 ].

As pointed out in [42 , 36 , 41 , 43 ], spiky features in collapsing inhomogeneous cosmologies will cause any fixed spatial resolution to become inadequate. Such evolutions are therefore candidates for adaptive mesh refinement such as that implemented by Hern and Stuart [143 , 142 ].