2.4 Hyperbolicity and numerical simulations

Knowing that the field equations for GR can be cast in a symmetric hyperbolic form, we can now ask how this fact can be of help for numerical calculations, (besides the extremely important fact that the problem would be well posed and so tractable!). There are at least two reasons why one should use variables in which the system is hyperbolic when performing numerical simulations. The first reason is that having a strongly hyperbolic system allows for standard constructions of numerical codes which are stable. If the system is symmetric, hyperbolic codes with better properties can be constructed. In particular, schema for symmetric hyperbolic systems can be constructed with numerical dissipation of lower order than what is needed for generic strongly hyperbolic ones (see Chapter VI of [36]). In variables where the system becomes diagonal, one can also use methods which take advantage of that structure.

The second reason for using a hyperbolic formulation for numerical analysis is that well posedness of the system gives bounds on the growth of the solution and its derivatives, as long as the solution is smooth. This property, when used in conjuction with stable algorithms, implies that one can bound the errors made on the simulation. That is, one knows not only that the error goes to zero as some power of the step size, but also the proportionality factor of that power law. In simulating phenomena whose observation would require hundreds of millions of dollars, a tight control of the accuracy reached should be required. Nevertheless it should be noted that raw hyperbolicity estimates alone usually give exponential bounds with very large growth coefficients, and that they are not of much value for numerical work.

Many of the systems we shall analyze can be cast as flux conservative equations with sources. This is a direct consequence of the facts that the principal part of the equations depends only on the metric variable, and that the equation for the time derivative of it does not contain derivatives of any of the dynamical variables. This property is important when using codes with variable grid spacing, even more if one considers that there are many standard codes for fluids – which are truly flux conserving – with adaptive grid schema.

It has to be said that flux conservation is important when dealing with systems that develop shock waves, that is in convective or more precisely in genuinely non-linear systems (for a definition of this term and many of the results, see [4353]). One should be cautious about any expectation of improvement by using flux conservative properties in general relativity, since here the shocks would probably not develop – in particular the systems are not genuinely non-linear. Rather, when singularities appear, they would be much worse than mere discontinuities of some of the dynamical fields. Due to bad gauge choices, discontinuities resembling shocks have been observed in numerical simulations, see [5]. Perhaps, instead of trying to devise an algorithm which allows one to go through these discontinuities, one should concentrate on finding better gauges, where it could even happen that the system cannot be put in flux conserving form. Thus, it is not clear whether flux conservative forms are relevant for vacuum general relativity.

Going Further

For the reader wishing to delve further into this precious theory of hyperbolic systems, while keeping a physicist’s approach, I recommend [34]. For those wishing to see more of the machinery at work, I recommend the book of Kreiss and Lorentz [48]. Finally, for those who really want to get the latest on the technical aspects and the modern approach to the problem, I recommend Taylor’s book [54]. Considerations about numerical analysis and algorithms can be found in [36]. In particular, that book contains general stable algorithms for strongly and symmetric hyperbolic systems and numerical error bounds in terms of analytic bounds of the exact solutions applicable to non-linear systems.