Structure Preserving Discretizations of the Liouville Equation and their Numerical Tests

The main purpose of this article is to show how symmetry structures in partial differential equations can be preserved in a discrete world and reflected in difference schemes. Three different structure preserving discretizations of the Liouville equation are presented and then used to solve specific boundary value problems. The results are compared with exact solutions satisfying the same boundary conditions. All three discretizations are on four point lattices. One preserves linearizability of the equation, another the infinite-dimensional symmetry group as higher symmetries, the third one preserves the maximal finite-dimensional subgroup of the symmetry group as point symmetries. A 9-point invariant scheme that gives a better approximation of the equation, but significantly worse numerical results for solutions is presented and discussed.


Introduction
This article is part of a general program the aim of which is to make full use of the theory of continuous groups (Lie groups) to study the solution space of discrete equations and in particular to solve difference equations [6,13,17,23]. This is one of the areas to which Luc Vinet made important original contributions [8][9][10]15].
Several recent articles [1,14,22] were devoted to discretizations of the Liouville equation [19] z xy = e z , (1.1) or its algebraic version u u xy − u x u y = u 3 , u = e z . (1. 2) The Liouville equation is of interest for many reasons. In differential geometry it is the equation satisfied by the conformal factor z(x, y) of the metric ds 2 = z 2 (dx 2 + dy 2 ) of a twodimensional space of constant curvature [7]. In the theory of infinite dimensional nonlinear integrable systems it is the prototype of a nonlinear partial differential equation (PDE) linearizable by a transformation of variables, involving the dependent variables (and their first derivatives) alone [19] In Lie theory this is probably the simplest PDE that has an infinite-dimensional Lie point symmetry group [21]. The symmetry algebra of the algebraic Liouville equation (1.2) is given by the vector fields Y (g(y)) = g(y)∂ y − g y (y)u∂ u , (1.4) where f (x) and g(y) are arbitrary smooth functions. Eq. (1.4) is a standard realization of the direct product of two centerless Virasoro algebras and we shall denote the corresponding Lie group V IR (x) ⊗ V IR (y). Restricting f (x) and g(y) to second order polynomials we obtain the maximal finite-dimensional subalgebra sl x (2, R) sl y (2, R) and the corresponding finite dimensional subgroup SL x (2, R) ⊗ SL y (2, R) of the symmetry group.
The Liouville equation is also an excellent tool for testing numerical methods for solving PDE's, since eq. (1.3) provides a very large class of exact analytic solutions, obtained by putting φ(x, y) = φ 1 (x) + φ 2 (y), (1.5) where φ 1 (x) and φ 2 (x) are arbitrary C (2) (I) functions on some interval I.
In [1] Adler and Startsev presented a discrete Liouville equation that preserves the property of being linearizable. In [22] Rebelo and Valiquette wrote a discrete Liouville equation that has the same infinite dimensional V IR (x) ⊗ V IR (y) symmetry group as the continuous Liouville equation. The transformations are however generalized symmetries, rather than point ones. In our article [14] we presented a discretization on a four-point stencil that preserves the maximal finite-dimensional subgroup of the V IR (x) ⊗ V IR (y) group as point symmetries. It was also shown that it is not possible to conserve the entire infinite dimensional Lie group of the Liouville equation as point symmetries. In [14] we also compared numerical solutions obtained using standard (non invariant) discretizations, the Rebelo-Valiquette invariant discretization [22] and our discretization with exact solutions (for 3 different specific solutions). It turned out that the discretization based on preserving the maximal subgroup of point transformations always gave the most accurate results for the considered solutions (all of them strictly positive in the area of integration).
The purpose of this article is to further explore and compare the different discratizations of the Liouville equation from two points of view. One is a theoretical one, namely to investigate the degree to which different discretizations preserve the qualitative feature of the equation: its exact linearizability, its infinite dimensional Lie point symmetry algebra, the behavior of the zeroes of the solutions. The other point of view is that of geometric integration: what are the advantages and disadvantages of the different discretizations as tools for obtaining numerical solutions.
In Section 2 we reproduce our previous [14] SL x (2, R) ⊗ SL y (2, R) symmetry preserving discretization using a 4-point stencil and show that in theory it can reproduce solutions that have horizontal or vertical lines of zeroes (or both). In Section 3 we propose an alternative discretization, using a 9-point stencil, instead of the 4-point one. It approximates the continuous Liouville equation with 2 precision, as opposed to the precision of the 4-point discretization. We show that increasing the number of points does not allow us to preserve the entire infinitedimensional symmetry algebra, nor to treat the lines of zeroes of solutions in a satisfactory manner. The Adler-Startsev discretization [1] is reproduced in Section 4 in a form suitable for numerical calculations. Section 5 is devoted to numerical experiments. Four different exact solutions of the algebraic Liouville equation are presented and then used to calculate boundary conditions on two lines parallel to the coordinate axes. The solutions are then calculated numerically using four different discretizations. We compare the goodness of the different methods and their qualitative features.

Point symmetries on a four point lattice and solutions with zeroes
In our previous article [14] we discretized the algebraic Liouville equation (1.2) on a four point regular orthogonal lattice preserving the SL x (2, R) ⊗ SL y (2, R) subgroup of its Lie point symmetries. The discretization was shown to provide good numerical results for solutions that were strictly positive in the entire integration region (a quadrant to the right and above a chosen point (x 0 , y 0 ), i.e. for x ≥ x 0 , y ≥ y 0 ). A particular property of the Liouville equation is that the zeroes of its solutions are not isolated. They occur on lines parallel to the x or y axes. Indeed, consider the infinite family of solutions of (1.2) parametrized by two arbitrary smooth functions of one variable φ 1 (x), φ 2 (y) (1.5). We take a region in which we have φ 1 (x) + φ 2 (y) = 0. Zeroes of u(x, y) occur if φ 1,x (x), or φ 2,y (y) are zero at some point x s , or y s (or both), respectively. We then have u(x s , y) = 0, ∀y, or u(x, y s ) = 0, ∀x . (2.1) This must be reflected in any computational scheme and the value u(x, y) = 0 will also occur on the intersection with the corresponding coordinate axis.
In [14] we considered several different boundary value problems. Here we restrict to the case of boundary conditions given on the lines x ≥ x 0 , y ≥ y 0 parallel to the coordinate axis. We can impose u(x s , 0) = 0 and/or u(0, y s ) = 0 (2.2) in order to obtain a solution satisfying (2.1).
The SL x (2, R)⊗SL y (2, R) invariants used in [14] to describe both the lattice and the discrete algebraic Liouville equation on a four point stencil were The lattice equations are satisfied by the uniform orthogonal lattice x m,n = hm + x 0 , y m,n = kn + y 0 (2.6) where the scale factors h and k are the same as in (2.4). The continuous limit corresponds to h → 0, k → 0. Two further independent SL x (2, R) ⊗ SL y (2, R) invariants exist on the four point stencil but any combination of them will either vanish, or be infinite on the lattice given by (2.5) (see [14]). The Liouville equation (1.2) was approximated in [14] by the difference scheme Eq. (2.7) can be solved for u m+1,n+1 in terms of u n,m , u n+1,m and u n,m+1 . On the first stencil we have m = n = 0. The boundary conditions are u m,0 = f (m) and u 0,n = g(n) with f and g given.
Let us now rewrite the recurrence relation (2.7) in terms of u m,n , choose b = d = 0, c = 1 − a, a ∈ R (in order to have an explicit scheme) and solve for u m+1,n+1 . We have

8)
A m,n+1;m+1,n = 1 + ahk |u m,n+1 u m+1,n | 1 + (a − 1)hk |u m,n+1 u m+1,n | . (2.9) We shall use (2.8, 2.9) to investigate the behavior of the numerical schemes for solutions that have rows (horizontal lines) or column (vertical lines) of zeroes. We impose boundary conditions on the lines x = x s and y = y s . To see the influence of the boundary conditions we introduce small quantities µ and ν on the coordinate axes that shall later be set to zero. We shall see that these small values do not propagate elsewhere but are confined to the columns and rows where they were introduced. This procedure is analogous to "singularity confinement" [11,12] used as an integrability criterion for difference equations.
We first note that we have Using (2.8, 2.10) we obtain expressions for u m 0 ,n , namely In the column to the right of the zeroes we obtain two equivalent expressions: (2.14) Thus the zero quantity µ cancels out and u m 0 +1,n is finite and nonzero for all n ≥ 0. Moreover u m 0 +1,n is expressed in terms of the given initial values and values calculated at previous nonzero values.

2.
A row of zeroes can be treated completely analogously.
The boundary conditions are replaced by u 0,n 0 = ν, u 0,n = 0 for n = n 0 , and we obtain i.e. a row of zeroes for ν = 0. The row above the zeroes satisfies 3. Two intersecting lines of zeroes.
The boundary conditions are Using the same considerations as above we find a column and a row of zeroes satisfying Thus, for µ = 0, ν = 0 the solutions u m,n have zeroes precisely where they should. Now let us use (2.8, 2.10) to calculate the values of u m 0 +1,n and u m,n 0 +1 , i.e. the column at the right and the row above the zeroes. The final result is that (2.13) is valid for all n = n 0 , n 0 + 1 and (2.17) for all m = m 0 , m 0 + 1 with Finally we see that the zeroes are confined to the rows and columns determined by a zero in the boundary condition and that the values of u m,n everywhere else are finite, non zero and determined by the equations (2.8, 2.9) and the boundary conditions.

Invariant discretization of the algebraic Liouville equation on a 9-point stencil
In Section 2 and in [14] we have shown that the Liouville equation can be approximated on 4 points. To approximate an arbitrary second order PDE for a function u(x, y) we need at least 6 points. An invariant discretization may need more than six.
and thus provides a first order approximation of the algebraic Liouville equation. In this section we will explore a second order approximation (order O(h 2 , k 2 , hk)) of the equation (1.2). To do this we shall use a 9-point stencil as shown on Fig. 1 , y 00 = y, y 01 = y + k 01 , y 10 = y + δ 10 , y 11 = y + k 01 + δ 11 , y 02 = y + k 01 + k 02 , y 20 = y + δ 10 + δ 20 , y 12 = y + k 01 + k 02 + δ 12 , y 21 = y + k 01 + δ 10 + δ 20 , of rectangle I we could use any other 4 points, and we shall use the vertices of the rectangles II, III and IV . The invariants involving the independent variables ξ a and η a (a = 1, · · · , 4) all vanish on the orthogonal lattice (2.6). The invariants depending on the dependent variables u ij that are finite and nonzero on this lattice are x ) + · · · ] We see that u 22 , u 02 , u 20 and u 00 figure only once each in the invariants, namely in J 8 , J 5 , J 3 and J 2 , respectively. On the other hand u 01 , u 10 , u 12 , and u 21 figure twice each, respectively in (J 1 , J 6 ), (J 1 , J 4 ), (J 6 , J 7 ) and (J 4 , J 7 ). The value u 11 figures in all four of J 2 , J 3 , J 5 and J 8 .
To lowest order we have To obtain the left hand side of the algebraic Liouville equation (1.2) up to order O(h 2 , hk, k 2 ) we need the differences J 2a − J 2a−1 to a higher order than in (3.3), namely In [14] equation (1.2) was approximated to order O(h, k). To approximate it to O(h 2 , hk, k 2 ) we must get rid of the terms of order O(h, k) in (3.6). The left hand side is approximated to the needed order by where α and β are arbitrary real constants.
To express the right hand side of (1.2) we use the basis: where R 1 , · · · , R 6 are all of the order O(h 2 , hk, k 2 ). The left hand side of (1.2) is already expressed in this basis (see (3.7) using B 3 , · · · , B 6 ). From the basis elements (3.8) we can calculate u 2 as with 5 free real parameters c i . To obtain u 3 we have several possibilities. One is to take The corresponding discrete Liouville equation is then Another possibility is to replace the basis (3.8) by We can then approximate the RHS of the discrete algebraic Liouville equation by In any case we have a large number of free parameters that can be chosen a priori to simplify calculations. The choice of the parameters (α, β, c 1 , · · · , c 6 ) or (α, β, γ a,b ) is restricted by the type of boundary conditions we wish to impose.
The quantity u 22 figures in J 8 only. An explicit scheme is obtained if J 8 figures linearly in the corresponding invariant discrete Liouville equation.
One possibility is to choose α = −3β and c 2 = c 3 = c 4 = c 5 = 0, c 6 = 1/2 in (3.11, 3.12). Then β = −1/8 and the invariant Liouville equation reduces to: In term of the field u ij (3.16) reads: Another simple possibility is to choose α = −3β and γ ab = δ a1 δ b6 in (3.15). Then we have β = −1/8 and we obtain In term of the field u ij (3.18) reads: Equations (3.17) and (3.19) are to be viewed as recursion relations, expressing u 2,2 in terms of the closest 8 points on a rectangle of which the point (22) is the top right vertex (see Fig.2).
By construction (3.17) and (3.19) are better approximations of the equation (1.1) than is (2.7). This does not mean that they will provide better numerical results and some comments are in order 1. Boundary conditions for a numerical solution on a 4-point lattice require the knowledge of u(x, y) on two lines, e.g. u m,0 and u 0,n , i.e. u(x, 0) and u(0, y). On the 9-point lattice we must start with 2 sets of parallel lines, e.g. u m,0 , u m,1 and u 0,n , u 1,n . This amounts to giving u(x, 0), u(0, y) and the first term of u y (x, 0), u x (0, y). This is more information than is needed in standard (non invariant) discretizations and indeed more information than is needed in theory to determine a solution completely. Hence once u(x, 0) and u(0, y) are given u y (x, 0) and u x (0, y) cannot be chosen arbitrarily. Thus u 22 is not strictly zero for 1 = 2 = 0, it does however satisfy u 22 ∼ O(h 2 , k 2 , hk). This is acceptable, however the problem arises when we shift the stencil and calculate u 32 which is supposed to be finite and nonzero if we assume u 30 = 0, u 31 = 0. What we obtain from (3.17) is Thus, u 32 is singular for 2 = 0 and becomes finite only in the continuous limit h = k = 0. This will quite obviously create numerical instabilities. They are avoided only for very special initial conditions, such that u 22 = 0 for all h and k. Using (3.19) leads to the same kind of problems.
Numerical results for several exact solutions showed that serious instabilities occur for the 9point scheme and thus we abandon it in favour of the 4-point one.

The Adler-Startsev linearizable discrete Liouville equation
Adler where b m,n satisfies the linear equation Hence the general solution of (4.1) is where c m , k n are arbitrary functions of one index each.
In [14] we showed, following [1] that the continuous limit of (4.1), for a m,n = − hk 2 u m,n when h and k go to zero, gives (1.2) and that it has no continuous point symmetries but must have generalized symmetries. Moreover by defining c m = φ 1 (x m,n ), k n = φ 2 (y m,n ) with x m,n and y m,n defined in (2.6) we have and thus a m,n = −hk a first order approximation of the general solution of (1.2) given by (1.3) .

Numerical experiments
In this section we shall apply the invariant recursion formula (2.8, 2.9) to solve a set of boundary value problems on a quadrant in the xy-plane. Boundary conditions will be given on two orthogonal lines parallel to the x and y axes and numerical solutions will be constructed above and to the right of these lines. The numerical solutions will be compared with exact solutions of the continuous equation for the same boundary conditions. In practice we will start from exact solutions given by choosing φ(x, y) = φ 1 (x) + φ 2 (y) in (1.3) and calculate the values of these functions on the boundaries. The global estimator which we use is the discrete analog of relative distance in L 2 D . We compute the quantity where F ij are the values of the exact solution F on the lattice sites and F α ij , with α = Inv, AS, RV, or stand are the values computed numerically for the invariant, Adler-Startsev, Rebello-Valiquette or standard discretization, respectively. The summation will be over all points of the lattice for which the calculation was performed.
We will compare results using four different discretization methods and thus four different recursion formulae, expressing u m+1,n+1 in terms of u m,n , u m+1,n and u m,n+1 . For comparison we present the four formulae for the first position of the stencil, i.e. m = n = 0. In all cases the left hand side of (1.2) is approximated by where h and k are the lengths of the steps in the x and y directions, respectively. The right hand side of (1.2) is approximated differently in each case. The corresponding recursion formulae are: 1. The invariant method (2.8, 2.9) (preserving the SL x (2, R) ⊗ SL y (2, R) symmetry group as point symmetries): We mention that the right hand sides of (5.4-5.6) are polynomials whereas the invariant case (5.3) involves square roots. These square roots may cause certain numerical instabilities for x and y close to the lines of zeroes of u(x, y). To address this problem we have made use of the results of Section 2. Namely if u m 0 ,n and u m 0 −1,n have different signs (a column of zeroes is between them), we replace the calculated u m 0 +1,n by the expression (2.14). Similarly we use We summarize our calculations in Table 1. One sees that the quality of the approximation Table 1: Relative mean square distance (5.1) between the numerical solutions and the analytic one.   We see from Table 2 that for all four methods the value of χ decreases faster then linearly as h = k decreases linearly. For this solution (with no zeroes) the value of χ Inv and χ AS are comparable and at least two order of magnitude lower than the other two. The values of χ RV are always lower than χ stand but of the same order.
The results presented in Table 1 clearly show that the performance of the different discretization procedures as numerical methods greatly depends on the presence or absence of lines of zeroes in the integration region. In order to investigate this aspect we performed another set of calculations, in which the boundary values are assigned on different pairs of orthogonal semi-lines. In the calculations of u m+1,n+1 we may. or may not cross the lines of zeroes.
In Table 3 Table 3 is divided by horizontal lines into three sections. In the first section no line of zeroes is crossed during the integration. In the second section one line is crossed. In the third section both lines are crossed, so that the saddle point is included. We see from Table 3 that in the first section we have χ Inv ≈ χ AS , with χ AS slightly better for all three solutions. Similarly we have χ RV ≈ χ stand with χ RV Table 3: The relative mean square distance (5.1) between the numerical solutions and the analytic one for f 2 with (x s = −.25, y s = .25), f 3 with (x s = 2.407, y s = 1.218) andf 4 with (x s = 0.0, y s = 0.0), for the different boundary value problems, defined on two orthogonal semilines parallel to axes x, y. In the first column we first identify the function f a and then we give the values (x 0 , y 0 ) indicating the left bottom corner where the numerical calculation starts. The mesh constants are h = k = .02 and the integration is applied for 129 × 129 steps. The parameter in the invariant recursive formula is a = 1.
always better. Remarkably χ Inv ≈ χ AS is always about two order of magnitude lower than χ RV ≈ χ stand (in the first section). For all three solutions in all sections χ AS is excellent and is not much influenced by the zeroes. Similarly χ RV is always lower than χ stand in all sections, sometimes significantly so (for f 3 ). The method most influenced by the presence of zeroes in the integration domain is the invariant one. This is particularly evident for the rapidly varying solution f 2 where χ Inv in section 2 and 3 becomes comparable to χ stand (or worse). For f 3 and f 4 χ Inv becomes comparable to χ stand only after passing past the saddle point. Another observation is that for all functions f a the invariant method seems to accumulate errors while approaching a line of zeroes.The further from the zero the calculation starts the larger is the error. This is specially visible in Table 1 where the distance between the starting point (x 0 , y 0 ) and the saddle point (x s , y s ) is much larger than in the case of Table 3.

Conclusions.
Both from the point of physics and from the point of view of geometric integration we see that for discretizing the Liouville equation we have to choose which characteristic feature of the equation we wish to preserve. Adler and Startsev [1] have shown how to preserve linearizability and the existence of a class of exact solutions depending on two arbitrary functions of one variable. We have shown that for a wide class of solutions a recurrence formula based on their method provides the most accurate results. On the other hand, linearizability, just like integrability is a property of a very restricted class of nonlinear PDEs.
The existence of a nontrivial Lie point symmetry group is a much more generic property, specially for PDEs coming from fundamental physical theories. From this point of view the Liouville equation is again special: its Lie point symmetry group is infinite dimensional. Rebelo and Valiquette [22] have presented a discretization that preserves this entire infinite-dimensional symmetry group as a special type of generalized symmetries. As opposed to more general higher symmetries, their symmetries have a global group action and are very interesting from the theoretical point of view. From the numerical one we have shown that the goodness of the solutions based on their recursion formula is systematically better (though of the same order of magnitude) than that of the standard method. The measure of goodness is the quantity χ of (5.1).
Finally, the method proposed in [14] and further developed in this article preserves point invariance under the maximal finite subgroup of the infinite dimensional symmetry group. Numerical methods based on this partial preservation of symmetries perform extremely well for solutions not vanishing on any line in the integration region, or for boundary conditions such that no line of zeroes is crossed in the calculation. On the other hand, if the integration region includes a line of zeroes, or specially two orthogonal lines of zeroes, the quality of the numerical solutions may be worse than for the standard method.
In future work we plan to study symmetry preserving discretizations of other equations with infinite-dimensional symmetry groups, such as the Kadomtsev-Petviashvili equation [3,4], the three-wave interaction equation [20] and specific Liouville type equations [25]. A symmetry preserving discretization of the Korteweg-de Vries equation has provided encouraging results [2].