Structure Preserving Discretizations of the Liouville Equation and their
Numerical Tests

Structure Preserving Discretizations
of the Liouville Equation and their Numerical Tests\star\star$\star$This paper
is a contribution to the Special Issue on Exact Solvability and Symmetry
Avatars in honour of Luc Vinet. The full collection is available at
http://www.emis.de/journals/SIGMA/ESSA2014.html

Decio LEVI †, Luigi MARTINA ‡ and Pavel WINTERNITZ †§

D. Levi, L. Martina and P. Winternitz

† Mathematics and Physics Department, Roma Tre University and Sezione INFN of
Roma Tre,
† Via della Vasca Navale 84, I-00146 Roma, Italy<EMAIL_ADDRESS>
‡ Dipartimento di Matematica e Fisica - Università del Salento and Sezione
INFN of Lecce,
‡ Via per Arnesano, C.P. 193 I-73100 Lecce, Italy<EMAIL_ADDRESS>
§ Département de mathématiques et de statistique and Centre de recherches
mathématiques,
§ Université de Montréal, C.P. 6128, succ. Centre-ville, Montréal (QC) H3C
3J7,
§ Canada (permanent address<EMAIL_ADDRESS>
Received April 01, 2015, in final form September 22, 2015; Published online
October 02, 2015

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.

Lie algebras of Lie groups; integrable systems; partial differential
equations; discretization procedures for PDEs

17B80; 22E60; 39A14; 65Mxx

Dedicated to Luc Vinet on the occasion of his 60th birthday.

##  1 Introduction

This article is part of a general program the aim of which is to make full use
of the theory of Lie groups to study the solution space of discrete equations
and in particular to solve difference equations [7, 8, 21, 24, 28, 29, 30,
42]. This is one of the areas to which Luc Vinet made important contributions
[14, 15, 16, 27].

S. Lie introduced what is now called Lie groups as groups of transformations
of the independent and dependent variables figuring in a system of
differential equations [31, 36]. Of special importance are symmetry groups,
transforming solutions into solutions. These may be point transformations,
where new variables depend only on the old ones. They may be contact
transformations, where the new variables depend also on the first derivatives
of the dependent variables. They may also be generalized symmetries where the
new variables can also depend on all derivatives of the old ones.

Lie’s method is particularly powerful for ordinary differential equations
(ODEs). A one-dimensional point symmetry group can be used to lower the order
of the ODE by one. An $n$-dimensional (solvable) Lie point symmetry group can
be used to decrease the order by $n$. Thus, if the ODE is of order
$n^{\prime}\leq n$ Lie group theory can provide the general solution (i.e.,
one that satisfies arbitrary initial conditions) in explicit or implicit
analytic form. For partial differential equations (PDEs) the Lie point
symmetry group is used to decrease the number of independent variables in the
equation and to provide special solutions (group invariant solutions),
satisfying particularly symmetrical boundary conditions.

The aim of this general program is to extend the use of Lie symmetry groups to
difference systems ($\Delta$S), i.e., to difference equations together with
the lattice they are written on.

The program has two complementary aspects, an analytical and a numerical one.

The aim of the analytical aspect is to determine the maximal symmetry group of
the $\Delta$S, i.e., the group of transformations that takes solutions into
solutions, and then to use it to obtain exact analytic solutions, at least
special ones, if possible general ones. The $\Delta$S to which the approach is
applied can come from the study of discrete physical, chemical, biological or
other systems, for which symmetries play an important role. Among them we
mention phenomena in crystals, or in atomic or molecular chains.

On the other hand $\Delta$S can be obtained by discretizing ODEs, or PDEs,
that have nontrivial symmetry groups reflecting fundamental physical laws such
as Galilei, Lorentz, or conformal invariance. At the scale of the Planck
length space-time may very well be discrete. In this case continuous equations
are approximations (continuous limits) of discrete ones. From the physical
point of view the symmetries are very important and should be preserved, e.g.,
when studying quantum field theories on lattices.

One way of preserving symmetries in a discretization of continuous equations
(the one used in this article) is to use symmetry adapted lattices that
themselves transform under the group action. This greatly enlarges the set of
equations for which symmetry preserving discretization is possible. We will
however see that in some cases only a subgroup of the Lie point symmetry group
can be preserved as point symmetries.

The numerical aspect of our program is the following. When solving an ODE or
PDE numerically it is always necessary to replace the continuous equation by a
difference system. This can be done in a standard manner, applicable to all
equations, simply by replacing derivatives by discrete derivatives. The other
possibility takes us directly into the field of geometric integration [20, 22,
33, 34]. The idea is to focus on some important feature of the underlying
problem and to preserve it in the discretization. Such a feature may be, for
instance linearizability, hamiltonian structure, integrability in the sense of
the existence of a Lax pairs and generalized symmetries or point and contact
symmetries. We are concentrating on point symmetries and exploring the
possibility and usefulness of including them in numerical calculations.

Earlier work has shown that for first-order ODEs preserving a one-dimensional
symmetry group provides an exact discretization [40]. For second-order ODEs
preserving a 3-dimensional symmetry group often provides analytically solvable
schemes (either via a Lagrangian [11, 12] or via the adjoint equation method
[9]). For third- and higher-order ODEs symmetry preserving discretization
provides numerical solutions that are, usually, closer to exact ones then
those obtained by other methods, specially near to the singularities [5, 39].
For previous work on PDEs see [2, 3, 4, 6, 10, 19, 25, 26, 37, 38, 41].

Several recent articles [1, 23, 38] were devoted to discretizations of the
Liouville equation [32]

$\displaystyle z_{xy}=e^{z},$ (1.1)

or its algebraic version

$\displaystyle uu_{xy}-u_{x}u_{y}=u^{3},\qquad 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 two-dimensional space of constant
curvature [13]. 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 [32]

$\displaystyle u=2\frac{\phi_{x}\phi_{y}}{\phi^{2}},\qquad\phi_{xy}=0.$ (1.3)

In Lie theory this is probably the simplest PDE that has an infinite-
dimensional Lie point symmetry group [35]. The symmetry algebra of the
algebraic Liouville equation (1.2) is given by the vector fields

$\displaystyle X(f(x))=f(x)\partial_{x}-f_{x}(x)u\partial_{u},\qquad
Y(g(y))=g(y)\partial_{y}-g_{y}(y)u\partial_{u},$ (1.4)

where $f(x)$ and $g(y)$ are arbitrary smooth functions.

Equation (1.4) is a standard realization of the direct product of two
centerless Virasoro algebras and we shall denote the corresponding Lie group
${\rm VIR}(x)\otimes{\rm VIR}(y)$. Restricting $f(x)$ and $g(y)$ to second-
order polynomials we obtain the maximal finite-dimensional subalgebra
${\mathfrak{sl}}_{x}(2,\mathbb{R})\bigoplus{\mathfrak{sl}}_{y}(2,\mathbb{R})$
and the corresponding finite-dimensional subgroup ${\rm
SL}_{x}(2,\mathbb{R})\otimes{\rm SL}_{y}(2,\mathbb{R})$ of the symmetry group.

The Liouville equation is also an excellent tool for testing numerical methods
for solving PDE’s, since equation (1.3) provides a very large class of exact
analytic solutions, obtained by putting

$\displaystyle\phi(x,y)=\phi_{1}(x)+\phi_{2}(y),$ (1.5)

where $\phi_{1}(x)$ and $\phi_{2}(x)$ are arbitrary $\mathbb{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 and exactly solvable. In [38]
Rebelo and Valiquette wrote a discrete Liouville equation that has the same
infinite-dimensional ${\rm VIR}(x)\otimes{\rm VIR}(y)$ symmetry group as the
continuous Liouville equation. The transformations are however generalized
symmetries, rather than point ones. In our article [23] we presented a
discretization on a four-point stencil that preserves the maximal finite-
dimensional subgroup of the ${\rm VIR}(x)\otimes{\rm VIR}(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 [23] we also compared numerical solutions obtained using standard (non
invariant) discretizations, the Rebelo–Valiquette invariant discretization
[38] 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
discretizations 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 [23] ${\rm
SL}_{x}(2,\mathbb{R})\otimes{\rm SL}_{y}(2,\mathbb{R})$ symmetry preserving
discretization using a 4-point stencil and show that after a slight
modification 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 $\epsilon^{2}$ precision, as opposed to the
$\epsilon$ precision of the 4-point discretization. We show that increasing
the number of points does not allow us to preserve the entire infinite-
dimensional symmetry algebra, nor to treat the lines of zeroes of solutions in
a satisfactory manner. Further, in Section 3 we take a specific exact solution
of the continuous algebraic Liouville equation (1.2) and approximate it on a
9-point lattice by a numerical solution. The Adler–Startsev discretization [1]
is reproduced in Section 4 in a form suitable for numerical calculations.
Section 5 is devoted to numerical tests of the invariant 4-point scheme. Five
different exact solution of the algebraic Liouville equation are presented and
then used to calculate boundary conditions on two lines parallel to the $x$
and $y$ coordinate axes, respectively. The solutions are then calculated
numerically using four different discretizations. We compare the validity of
the different methods and their qualitative features. Some general
conclusions, placing this article in the context of geometric integration, are
presented in the final Section 6.

##  2 Point symmetries on a four point lattice
and solutions with zeroes

In our previous article [23] we discretized the algebraic Liouville equation
(1.2) on a four point regular orthogonal lattice preserving the ${\rm
SL}_{x}(2,\mathbb{R})\otimes{\rm SL}_{y}(2,\mathbb{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\geq x_{0}$, $y\geq 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 $\phi_{1}(x)$, $\phi_{2}(y)$
(1.5). We take a region in which we have $\phi_{1}(x)+\phi_{2}(y)\neq 0$.
Zeroes of $u(x,y)$ occur if $\phi_{1,x}(x)$, or $\phi_{2,y}(y)$ are zero at
some point $x_{s}$, or $y_{s}$ (or both), respectively. We then have

$\displaystyle u(x_{s},y)=0,\quad\forall\,y,\qquad\mbox{or}\qquad
u(x,y_{s})=0,\quad\forall\,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 [23] we considered several different boundary value problems. Here we
restrict to the case of boundary conditions given on the lines $x\geq x_{0}$,
$y\geq y_{0}$ parallel to the coordinate axes. We can impose

$\displaystyle u(x_{s},0)=0\qquad\mbox{and/or}\qquad u(0,y_{s})=0$

in order to obtain a solution satisfying (2.1).

The ${\rm SL}_{x}(2,\mathbb{R})\otimes{\rm SL}_{y}(2,\mathbb{R})$ invariants
used in [23] to describe both the lattice and the discrete algebraic Liouville
equation on a four point stencil were

$\displaystyle\xi_{1}=\frac{\left(x_{m,n+1}-x_{m,n}\right)\left(x_{m+1,n+1}-x_{m+1,n}\right)}{\left(x_{m,n}-x_{m+1,n}\right)\left(x_{m,n+1}-x_{m+1,n+1}\right)},$
$\displaystyle\eta_{1}=\frac{\left(y_{m,n}-y_{m+1,n}\right)\left(y_{m,n+1}-y_{m+1,n+1}\right)}{\left(y_{m,n+1}-y_{m,n}\right)\left(y_{m+1,n+1}-y_{m+1,n}\right)},$
(2.2) $\displaystyle J_{1}=u_{m+1,n}u_{m,n+1}h^{2}k^{2},\qquad
J_{2}=u_{m,n}u_{m+1,n+1}h^{2}k^{2}.$ (2.3)

The lattice equations

$\displaystyle\xi_{1}=0,\qquad\eta_{1}=0$ (2.4)

are satisfied by the uniform orthogonal lattice

$\displaystyle x_{m,n}=hm+x_{0},\qquad y_{m,n}=kn+y_{0},$ (2.5)

where the scale factors $h$ and $k$ are the same as in (2.3). The continuous
limit corresponds to $h\rightarrow 0$, $k\rightarrow 0$. Two further
independent ${\rm SL}_{x}(2,\mathbb{R})\otimes{\rm SL}_{y}(2,\mathbb{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.4) (see [23]).

The Liouville equation (1.2) was approximated in [23] by the difference scheme

$\displaystyle
J_{2}-J_{1}=a{\,\operatorname{sign}(J_{1})}|J_{1}|^{3/2}+b{\,\operatorname{sign}(J_{2})}J_{1}|J_{2}|^{1/2}$
$\displaystyle\hphantom{J_{2}-J_{1}=}{}+c{\,\operatorname{sign}(J_{1})}|J_{1}|^{1/2}J_{2}+d{\,\operatorname{sign}(J_{2})}|J_{2}|^{3/2},$
(2.6) $\displaystyle\xi_{1}=\eta_{1}=0,\qquad a+b+c+d=1.$

The symbols ${\,\operatorname{sign}(J_{1})}$ and
${\,\operatorname{sign}(J_{2})}$ were omitted in [23] and were not necessary
as we restricted our formulation to strictly positive solutions. Equation
(2.6) 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.6) in terms of $u_{m,n}$, choose
$b=d=0$, $c=1-a$, $a\in\mathbb{R}$ (in order to have an explicit scheme) and
solve for $u_{m+1,n+1}$. We have

$\displaystyle u_{m+1,n+1}=\frac{u_{m,n+1}u_{m+1,n}}{u_{m,n}}A_{m,n+1;m+1,n},$
(2.7) $\displaystyle
A_{m,n+1;m+1,n}=\frac{1+ahk{\,\operatorname{sign}(u_{m+1,n}u_{m,n+1})}\sqrt{|u_{m,n+1}u_{m+1,n}|}}{1+(a-1)hk{\,\operatorname{sign}(u_{m+1,n}u_{m,n+1})}\sqrt{|u_{m,n+1}u_{m+1,n}|}}.$
(2.8)

The expression ${\,\operatorname{sign}(u_{m+1,n}u_{m,n+1})}$ follows from
(2.6). Here the sign before the square root is important since it will change
when the sign of $u$ changes in the recurrence relation.

We shall use (2.7), (2.8) to investigate the behaviour of the numerical
schemes for solutions that have rows (horizontal lines) or columns (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 $\mu$ and $\nu$ on the coordinate axes that will 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” [17, 18] used as an integrability
criterion for difference equations.

We first note that we have

$\displaystyle\lim_{u_{m,n+1}\rightarrow
0}A_{m,n+1;m+1,n}=\lim_{u_{m+1,n}\rightarrow 0}A_{m,n+1;m+1,n}=1.$ (2.9)

Three cases will be considered separately:

1\. A column of zeroes. The boundary conditions are

$\displaystyle u_{m_{0},0}=\mu,\qquad u_{m,0}\neq 0\qquad\mbox{for}\quad m\neq
m_{0}.$

Using (2.7), (2.9) we obtain expressions for $u_{m_{0},n}$, namely

$\displaystyle u_{m_{0},n}=\frac{u_{m_{0}-1,n}}{u_{m_{0}-1,0}}\mu,\qquad n\geq
0.$

In the column to the right of the zeroes we obtain two equivalent expressions:

$\displaystyle
u_{m_{0}+1,n}=u_{m_{0}+1,n-1}\frac{u_{m_{0}-1,n}}{u_{m_{0}-1,n-1}},$ (2.10)
$\displaystyle
u_{m_{0}+1,n}=u_{m_{0}-1,n}\frac{u_{m_{0}+1,0}}{u_{m_{0}-1,0}}.$ (2.11)

Thus the zero quantity $\mu$ cancels out and $u_{m_{0}+1,n}$ is finite and
nonzero for all $n\geq 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

$\displaystyle u_{0,n_{0}}=\nu,\qquad u_{0,n}\neq 0\qquad\mbox{for}\quad n\neq
n_{0},$

and we obtain

$\displaystyle u_{m,n_{0}}=\frac{u_{m,n_{0}-1}}{u_{0,n_{0}-1}}\nu,\qquad m\geq
0,$

i.e., a row of zeroes for $\nu=0$. The row above the zeroes satisfies

$\displaystyle
u_{m,n_{0}+1}=u_{m-1,n_{0}+1}\frac{u_{m,n_{0}-1}}{u_{m-1,n_{0}-1}},$ (2.12)
$\displaystyle
u_{m,n_{0}+1}=u_{m,n_{0}-1}\frac{u_{0,n_{0}+1}}{u_{0,n_{0}-1}}.$ (2.13)

3\. Two intersecting lines of zeroes. The boundary conditions are

$\displaystyle u_{0,n_{0}}=\nu,\qquad u_{m_{0},0}=\mu;\qquad u_{m,0}\neq
0\qquad\mbox{for}\quad m\neq m_{0},\qquad u_{0,n}\neq 0\qquad\mbox{for}\quad
n\neq n_{0}.$

Using the same considerations as above we find a column and a row of zeroes
satisfying

$\displaystyle u_{m_{0},n}=\frac{u_{m_{0}-1,n}}{u_{m_{0}-1,0}}\mu,\qquad n\neq
n_{0},\qquad u_{m,n_{0}}=\frac{u_{m,n_{0}-1}}{u_{0,n_{0}-1}}\nu,\qquad m\neq
m_{0},$ $\displaystyle
u_{m_{0},n_{0}}=\frac{u_{m_{0}-1,n_{0}-1}}{u_{0,n_{0}-1}u_{m_{0}-1,0}}\mu\nu.$

Thus, for $\mu=0$, $\nu=0$ the solutions $u_{m,n}$ have zeroes precisely where
they should. Now let us use (2.7), (2.9) 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.10) is valid for all $n\neq
n_{0},n_{0}+1$ and (2.12) for all $m\neq m_{0},m_{0}+1$ with

$\displaystyle
u_{m_{0}+1,n_{0}+1}=u_{m_{0}-1,n_{0}+1}\frac{u_{m_{0}+1,n_{0}-1}}{u_{m_{0}-1,n_{0}-1}},$

while (2.11) is valid for all $n\neq n_{0}$ and (2.13) for all $m\neq m_{0}$.

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.7),
(2.8) and the boundary conditions. In other words the rows and columns of
zeroes do not interfere with the integration algorithm. This will be confirmed
by numerical calculations in Section 5.

##  3 Invariant discretization of the algebraic Liouville equation
using a larger number of points

There are several reasons to increase the number of points on the stencil that
we use.

1\. To determine whether the entire ${\rm VIR}(x)\otimes{\rm VIR}(y)$ symmetry
group can be preserved on a larger lattice.

2\. To determine whether the only other ${\rm SL}_{x}(2,\mathbb{R})\otimes{\rm
SL}_{y}(2,\mathbb{R})$ differential invariant [23], namely

$\displaystyle
I_{2}=\frac{1}{u^{6}}\big{(}2uu_{xx}-3u_{x}^{2}\big{)}\big{(}2uu_{yy}-3u_{y}^{2}\big{)}$
(3.1)

can be invariantly discretized on a larger lattice. Four points are clearly
not sufficient to approximate two first- and three second-order derivatives.

3\. To approximate the algebraic Liouville equation with a higher degree of
accuracy in $h$ and $k$ and thus possibly improve the numerical calculations.

In Section 2 and in [23] 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.

Figure 1: Points on a general lattice, e.g., $x_{00}=x$, $x_{10}=x+h_{10}$,
$x_{01}=x+\epsilon_{01}$, $x_{11}=x+h_{10}+\epsilon_{11}$,
$x_{20}=x+h_{10}+h_{20}$, $x_{02}=x+\epsilon_{01}+\epsilon_{02}$,
$x_{12}=x+h_{10}+\epsilon_{11}+\epsilon_{12}$,
$x_{21}=x+h_{10}+h_{20}+\epsilon_{21}$,
$x_{22}=x+h_{10}+h_{20}+\epsilon_{21}+\epsilon_{22}$, $y_{00}=y$,
$y_{01}=y+k_{01}$, $y_{10}=y+\delta_{10}$, $y_{11}=y+k_{01}+\delta_{11}$,
$y_{02}=y+k_{01}+k_{02}$, $y_{20}=y+\delta_{10}+\delta_{20}$,
$y_{12}=y+k_{01}+k_{02}+\delta_{12}$,
$y_{21}=y+k_{01}+\delta_{10}+\delta_{20}$,
$y_{22}=y+k_{01}+k_{02}+\delta_{12}+\delta_{22}$.  Figure 2: A stencil for the
9-points scheme.

Equation (2.6) satisfies

$\displaystyle\lim_{h,k\rightarrow
0}\frac{1}{h^{3}k^{3}}\big{\\{}J_{2}-J_{1}-\big{[}a|J_{1}|^{3/2}+bJ_{1}|J_{2}|^{1/2}+c|J_{1}|^{1/2}J_{2}+d|J_{2}|^{3/2}\big{]}\big{\\}}$
$\displaystyle\qquad{}=\big{[}uu_{xy}-u_{x}u_{y}-u^{3}\big{]}[1+\mathcal{O}(h,k)],$

and thus provides a first-order approximation of the algebraic Liouville
equation. In this section we will explore a second-order approximation (order
$\mathcal{O}(h^{2},k^{2},hk)$) of the equation (1.2). To do this we shall use
a 9-point stencil as shown on Figs. 1 and 2. The 4 well-behaved invariants
(2.2), (2.3) make use of the four vertices of rectangle I on Fig. 2. Instead
of the vertices 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 $\xi_{a}$ and $\eta_{a}$ ($a=1,\dots,4$) all vanish
on the orthogonal lattice (2.5). The invariants depending on the dependent
variables $u_{ij}$ that are finite and nonzero on this lattice are

$\displaystyle J_{1}=u_{01}u_{10}h^{2}k^{2},\qquad$ $\displaystyle
J_{2}=u_{00}u_{11}h^{2}k^{2},$ $\displaystyle
J_{3}=u_{11}u_{20}h^{2}k^{2},\qquad$ $\displaystyle
J_{4}=u_{10}u_{21}h^{2}k^{2},$ $\displaystyle
J_{5}=u_{11}u_{02}h^{2}k^{2},\qquad$ $\displaystyle
J_{6}=u_{01}u_{12}h^{2}k^{2},$ $\displaystyle
J_{7}=u_{12}u_{21}h^{2}k^{2},\qquad$ $\displaystyle
J_{8}=u_{11}u_{22}h^{2}k^{2}.$

The quantities $J_{1},\dots,J_{8}$ are linearly independent but one polynomial
relation exists between them, namely

$\displaystyle J_{4}J_{6}=J_{1}J_{7}.$

The continuous limit is obtained by expanding the invariants into Taylor
series and then taking $h\rightarrow 0$, $k\rightarrow 0$. We shall assume
that they tend to zero at the same rate, i.e., $k=\alpha h$, $\alpha\sim 1$.
We have

$\displaystyle
J_{1}=h^{2}k^{2}\big{[}u^{2}+huu_{x}+kuu_{y}+hku_{x}u_{y}+\tfrac{1}{2}h^{2}uu_{xx}+\tfrac{1}{2}k^{2}uu_{yy}+\cdots\big{]},$
$\displaystyle
J_{2}=h^{2}k^{2}\big{[}u^{2}+huu_{x}+kuu_{y}+hkuu_{xy}+\tfrac{1}{2}h^{2}uu_{xx}+\tfrac{1}{2}k^{2}uu_{yy}+\cdots\big{]},$
$\displaystyle
J_{3}=h^{2}k^{2}\big{[}u^{2}+3huu_{x}+kuu_{y}+\tfrac{1}{2}k^{2}uu_{yy}+hk(uu_{xy}+2u_{y}u_{x})+h^{2}\big{(}\tfrac{5}{2}uu_{xx}+2u_{x}^{2}\big{)}+\cdots\big{]},$
$\displaystyle
J_{4}=h^{2}k^{2}\big{[}u^{2}+3huu_{x}+kuu_{y}+\tfrac{1}{2}k^{2}uu_{yy}+hk(2uu_{xy}+u_{y}u_{x})+h^{2}\big{(}\tfrac{5}{2}uu_{xx}+2u_{x}^{2}\big{)}+\cdots\big{]},$
$\displaystyle
J_{5}=h^{2}k^{2}\big{[}u^{2}+huu_{x}+3kuu_{y}+k^{2}\big{(}\tfrac{5}{2}uu_{yy}+2u_{y}^{2}\big{)}+hk(uu_{xy}+2u_{y}u_{x})+\tfrac{1}{2}h^{2}uu_{xx}+\cdots\big{]},$
$\displaystyle
J_{6}=h^{2}k^{2}\big{[}u^{2}+huu_{x}+3kuu_{y}+k^{2}\big{(}\tfrac{5}{2}uu_{yy}+2u_{y}^{2}\big{)}+hk(2uu_{xy}+u_{y}u_{x})+\tfrac{1}{2}h^{2}uu_{xx}+\cdots\big{]},$
$\displaystyle
J_{7}=h^{2}k^{2}\big{[}u^{2}+3huu_{x}+3kuu_{y}+k^{2}\big{(}2u_{y}^{2}+\tfrac{5}{2}uu_{yy}\big{)}+hk(4uu_{xy}+5u_{y}u_{x})$
$\displaystyle\hphantom{J_{7}=}{}+h^{2}\big{(}2u_{x}^{2}+\tfrac{5}{2}uu_{xx}\big{)}+\cdots\big{]},$
$\displaystyle
J_{8}=h^{2}k^{2}\big{[}u^{2}+3huu_{x}+3kuu_{y}+k^{2}\big{(}\tfrac{5}{2}uu_{yy}+2u_{y}^{2}\big{)}+hk(5uu_{xy}+4u_{y}u_{x})$
$\displaystyle\hphantom{J_{8}=}{}+h^{2}\big{(}\tfrac{5}{2}uu_{xx}+2u_{x}^{2}\big{)}+\cdots].$
(3.2)

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

$\displaystyle
J_{2}-J_{1}=J_{4}-J_{3}=J_{6}-J_{5}=J_{8}-J_{7}=h^{3}k^{3}(uu_{xy}-u_{x}u_{y})[1+\mathcal{O}(h,k)].$

To obtain the left hand side of the algebraic Liouville equation (1.2) up to
order $\mathcal{O}(h^{2},hk,k^{2})$ we need the differences $J_{2a}-J_{2a-1}$
to a higher-order than in (3.2), namely

$\displaystyle
J_{2}-J_{1}=h^{3}k^{3}\left\\{(uu_{xy}-u_{x}u_{y})+\frac{h}{2}(uu_{xxy}-u_{y}u_{xx})+\frac{k}{2}(uu_{xyy}-u_{x}u_{yy})\right\\},$
$\displaystyle
J_{4}-J_{3}=h^{3}k^{3}\left\\{(uu_{xy}-u_{x}u_{y})+\frac{3h}{2}(uu_{xxy}-u_{y}u_{xx})+\frac{k}{2}(uu_{xyy}-u_{x}u_{yy})\right\\},$
$\displaystyle
J_{6}-J_{5}=h^{3}k^{3}\left\\{(uu_{xy}-u_{x}u_{y})+\frac{h}{2}(uu_{xxy}-u_{y}u_{xx})+\frac{3k}{2}(uu_{xyy}-u_{x}u_{yy})\right\\},$
$\displaystyle
J_{8}-J_{7}=h^{3}k^{3}\left\\{(uu_{xy}-u_{x}u_{y})+\frac{3h}{2}(uu_{xxy}-u_{y}u_{xx})+\frac{3k}{2}(uu_{xyy}-u_{x}u_{yy})\right\\}.$
(3.3)

In [23] equation (1.2) was approximated to order $\mathcal{O}(h,k)$. To
approximate it to $\mathcal{O}(h^{2},hk,k^{2})$ we must get rid of the terms
of order $\mathcal{O}(h,k)$ in (3.3).

The left hand side is approximated to the needed order by

$\displaystyle\alpha[4(J_{2}-J_{1})-(J_{6}-J_{5}+J_{4}-J_{3})]+\beta[4(J_{8}-J_{7})-3(J_{6}-J_{5}+J_{4}-J_{3})]$
$\displaystyle\qquad{}=2h^{3}k^{3}(uu_{xy}-u_{x}u_{y})(\alpha-\beta)\big{[}1+\mathcal{O}\big{(}h^{2},hk,k^{2}\big{)}\big{]},$
(3.4)

where $\alpha$ and $\beta$ are arbitrary real constants.

To express the right-hand side of (1.2) we use the basis

$\displaystyle B_{1}=\frac{1}{2}(3J_{1}-J_{8})=h^{2}k^{2}u^{2}(1+R_{1}),\qquad
B_{2}=J_{4}+J_{6}-J_{8}-J_{1}=h^{2}k^{2}R_{2},$ $\displaystyle
B_{3}=J_{2}-J_{1}=h^{2}k^{2}R_{3},\qquad
B_{4}=J_{4}-J_{3}=h^{2}k^{2}R_{4},\qquad B_{5}=J_{6}-J_{5}=h^{2}k^{2}R_{5},$
$\displaystyle B_{6}=J_{8}-J_{7}=h^{2}k^{2}R_{6},$ (3.5)

where $R_{1},\dots,R_{6}$ are all of the order $\mathcal{O}(h^{2},hk,k^{2})$.

The left hand side of (1.2) is already expressed in this basis (see (3.4)
using $B_{3},\dots,B_{6}$).

From the basis elements (3.5) we can calculate $u^{2}$ as

$\displaystyle
B_{1}+\sum_{i=2}^{6}c_{i}B_{i}=h^{2}k^{2}u^{2}\big{[}1+\mathcal{O}\big{(}h^{2},hk,k^{2}\big{)}\big{]},$

with 5 free real parameters $c_{i}$. To obtain $u^{3}$ we have several
possibilities. One is to take

$\displaystyle\left(B_{1}+\sum_{i=2}^{6}c_{i}B_{i}\right)^{3/2}=h^{3}k^{3}u^{3}\big{[}1+\mathcal{O}\big{(}h^{2},hk,k^{2}\big{)}\big{]}.$
(3.6)

The corresponding discrete Liouville equation is then

$\displaystyle\alpha[4B_{3}-(B_{4}+B_{5})]+\beta[4B_{6}-3(B_{4}+B_{5})]=\left(B_{1}+\sum_{i=2}^{6}c_{i}B_{i}\right)\\!\left\\{\left|B_{1}+\sum_{i=2}^{6}c_{i}B_{i}\right|\right\\}^{1/2}\\!\\!\\!,\\!\\!\\!$
(3.7)

with

$\displaystyle 2(\alpha-\beta)=1.$ (3.8)

Another possibility is to replace the basis (3.5) by

$\displaystyle A_{1}=B_{1},\qquad A_{a}=B_{1}+\frac{1}{2}B_{a},\qquad
a=2,\dots,6.$

We can then approximate the right-hand side of the discrete algebraic
Liouville equation by

$\displaystyle\sum_{a,b=1}^{6}\gamma_{a,b}A_{a}\sqrt{|A_{b}|}=h^{3}k^{3}u^{3}\sum_{a,b=1}^{6}\gamma_{a,b}\big{[}1+\mathcal{O}\big{(}h^{2},hk,k^{2}\big{)}\big{]}.$

Then the discrete algebraic Liouville equation reads

$\displaystyle
2\alpha[4A_{3}-2A_{1}-(A_{4}+A_{5})]+2\beta[4A_{6}+2A_{1}-3(A_{4}+A_{5})]=\sum_{a,b=1}^{6}\gamma_{a,b}A_{a}\sqrt{|A_{b}|},$
(3.9) $\displaystyle\sum_{a,b=1}^{6}\gamma_{a,b}=2(\alpha-\beta).$

The general invariant equations (3.9), (3.6) use all 9 points on the stencil
and contain a lot of free parameters. The parameters can be chosen to simplify
calculations, though the choice of $(\alpha,\beta,c_{1},\dots,c_{6})$ or
$(\alpha,\beta,\gamma_{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 $\alpha=-3\beta$ and $c_{2}=c_{3}=c_{4}=c_{5}=0$,
$c_{6}=1/2$ in (3.7), (3.8). Then $\beta=-1/8$ and the invariant Liouville
equation reduces to

$\displaystyle
J_{8}=J_{7}+3(J_{2}-J_{1})-\frac{1}{\sqrt{2}}(3J_{1}-J_{7})\sqrt{|3J_{1}-J_{7}|}.$
(3.10)

In terms of the field $u_{ij}$ (3.10) reads

$\displaystyle
u_{22}=\frac{1}{u_{11}}[u_{12}u_{21}+3(u_{11}u_{00}-u_{10}u_{01})$
$\displaystyle\hphantom{u_{22}=}{}-\frac{1}{\sqrt{2}}hk(3u_{01}u_{10}-u_{12}u_{21})\sqrt{|3u_{01}u_{10}-u_{12}u_{21}|}].$
(3.11)

so that $u_{22}$ is expressed in terms of $u_{00}$, $u_{01}$, $u_{10}$,
$u_{11}$, $u_{12}$ and $u_{21}$, i.e., only 7 points are involved.

Another simple possibility is to choose $\alpha=-3\beta$ and
$\gamma_{ab}=\delta_{a1}\delta_{b6}$ in (3.9). Then we have $\beta=-1/8$ and
we obtain

$\displaystyle
J_{8}=\frac{J_{7}+3(J_{2}-J_{1})-\frac{3}{\sqrt{2}}J_{1}{\,\operatorname{sign}(3J_{1}-J_{7})}\sqrt{|3J_{1}-J_{7}|}}{1-\frac{1}{\sqrt{2}}{\,\operatorname{sign}(3J_{1}-J_{7})}\sqrt{|3J_{1}-J_{7}|}}$
(3.12)

In terms of the field $u_{ij}$ (3.12) reads

$\displaystyle u_{22}=\bigg{\\{}u_{12}u_{21}+3(u_{11}u_{00}-u_{10}u_{01})$
$\displaystyle\hphantom{u_{22}=}{}-\frac{3}{\sqrt{2}}hku_{01}u_{10}{\,\operatorname{sign}(3u_{01}u_{10}-u_{12}u_{21})}\sqrt{|3u_{01}u_{10}-u_{12}u_{21}|}\bigg{\\}}$
$\displaystyle\hphantom{u_{22}=}{}\times\bigg{\\{}u_{11}\bigg{[}1-\frac{hk}{\sqrt{2}}{\,\operatorname{sign}(3u_{01}u_{10}-u_{12}u_{21})}\sqrt{|3u_{01}u_{10}-u_{12}u_{21}|}\bigg{]}\bigg{\\}}^{-1}.$
(3.13)

Again only 7 of the 9 points on a stencil are used.

Equations (3.11) and (3.13) are to be viewed as recursion relations,
expressing $u_{2,2}$ in terms of 6 points on a rectangle of which the point
(2,2) is the top right vertex (see Fig. 2).

By construction (3.11) and (3.13) are better approximations of the equation
(1.1) than is (2.6). 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. In our
numerical solutions we calculated the conditions using an exact solution so
$u(x,0)$, $u(0,y)$ and $u_{x}(0,y)$, $u_{y}(x,0)$ are consistent.

2\. Contrary to the case of a 4-point lattice, instabilities close to zero
lines of solutions cannot be avoided on 7- or 9-point lattices. Indeed let us
give initial conditions on the first square satisfying $u_{00}\neq 0$,
$u_{01}\neq 0$, $u_{10}\neq 0$, $u_{11}\neq 0$, $u_{12}\neq 0$,
$u_{20}=\epsilon_{1}$, $u_{21}=\epsilon_{2}$. From the known solution of the
PDE (1.1) we expect the solution to satisfy $u_{2,m}=0$ for $m\geq 2$.
Equation (3.11) implies

$\displaystyle
u_{22}=\frac{1}{u_{11}}\bigg{[}u_{12}\epsilon_{2}+3(u_{11}u_{00}-u_{10}u_{01})-\frac{1}{\sqrt{2}}hk(3u_{01}u_{10}-u_{12}\epsilon_{2})\sqrt{|3u_{01}u_{10}-u_{12}\epsilon_{2}|}\bigg{]}.$

Thus $u_{22}$ is not strictly zero for $\epsilon_{1}=\epsilon_{2}=0$, it does
however satisfy $u_{22}\sim\mathcal{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}\neq 0$,
$u_{31}\neq 0$. What we obtain from (3.11) is

$\displaystyle
u_{32}=\frac{1}{\epsilon_{2}}\bigg{[}\mathcal{O}\big{(}h^{2},k^{2},hk\big{)}u_{31}+3(\epsilon_{2}u_{10}-\epsilon_{1}u_{11})$
$\displaystyle\hphantom{u_{32}=}{}-\frac{1}{\sqrt{2}}hk\big{(}3\epsilon_{1}u_{11}-\mathcal{O}\big{(}h^{2},k^{2},hk\big{)}u_{31}\big{)}\sqrt{\big{|}3\epsilon_{1}u_{11}-\mathcal{O}\big{(}h^{2},k^{2},hk\big{)}u_{31}\big{|}}\bigg{]}.$

Thus, $u_{32}$ is singular for $\epsilon_{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.13) leads to the same kind of
problems.

Sadly (for the 9 points scheme) the answers to all three questions posed in
the beginning of this section are negative.

1\. The only ${\rm VIR}(x)\otimes{\rm VIR}(y)$ invariant lattice is given by
requiring conditions (2.2) in all four quadrangles of Fig. 2. Indeed we have

$\displaystyle\operatorname{pr}X\big{(}x^{3}\big{)}\xi_{a}=(x_{11}-x_{00})(x_{10}-x_{01})\xi_{a},$
$\displaystyle\operatorname{pr}Y\big{(}y^{3}\big{)}\eta_{a}=(y_{11}-y_{00})(y_{10}-y_{01})\eta_{a},\qquad
a=1,2,3,4$

so (weak) invariance requires $\xi_{a}=\eta_{a}=0$ and adding further points
does not help. Moreover, no function of $J_{1},\dots,J_{8}$ has the correct
continuous limit and is invariant under ${\rm VIR}(x)\otimes{\rm VIR}(y)$.

2\. We have determined that the invariant (3.1) cannot be discretized in an
${\rm SL}_{x}(2,\mathbb{R})\otimes{\rm SL}_{y}(2,\mathbb{R})$ invariant manner
using $J_{1},\dots,J_{8}$ but we do not present the details here since this
question is not related to the Liouville equation.

3\. Numerical results for several exact solutions showed that serious
instabilities occur for the 7-point scheme. As one can see in Fig. 3
representing the solution $f_{1}$ given in (5.9), instabilities in the 7-point
case occur almost immediately. The same is true for other solutions of
equation (1.2).

Our conclusion is that the 9-point (or 7-point) scheme is too unstable to be
useful. We present it here because we think that this negative result is not a
priori obvious and that this discussion may be useful.

Figure 3: Comparison between the 4-point and 7-point approximation for the
solution $f_{1}$ given in (5.9). The dark grey graph represents the results of
the application of (3.12), (3.13), while the light grey one is obtained by
using the 4-point formula (2.7), (2.8). The integration is carried out
starting from the the point $(-1.0,-1.0)$, the recursion formula is applied
for steps $h=k=0.02$ on a grid $20\times 20$.

##  4 The Adler–Startsev linearizable discrete Liouville equation

Adler and Startsev [1] have presented a discretization of the algebraic
Liouville equation (1.2) on a four-point lattice, namely

$\displaystyle
a_{m+1,n+1}\left(1+\frac{1}{a_{m+1,n}}\right)\left(1+\frac{1}{a_{m,n+1}}\right)a_{m,n}=1.$
(4.1)

This equation is linearizable by the substitution

$\displaystyle
a_{m,n}=-\frac{(b_{m+1,n}-b_{m,n})(b_{m,n+1}-b_{m,n})}{b_{m+1,n}b_{m,n+1}},$

where $b_{m,n}$ satisfies the linear equation

$\displaystyle b_{m+1,n+1}-b_{m+1,n}-b_{m,n+1}+b_{m,n}=0.$

Hence the exact general solution of (4.1) is

$\displaystyle
a_{m,n}=-\frac{(c_{m+1}-c_{m})(k_{n+1}-k_{n})}{(c_{m+1}+k_{n})(c_{m}+k_{n+1})},$

where $c_{m}$, $k_{n}$ are arbitrary functions of one index each.

In [23] we showed, following [1] that the continuous limit of (4.1), for
$a_{m,n}=-\frac{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}=\phi_{1}(x_{m,n})$,
$k_{n}=\phi_{2}(y_{m,n})$ with $x_{m,n}$ and $y_{m,n}$ defined in (2.5) we
have

$\displaystyle
c_{m+1}=\phi_{1}(x)+h\frac{d\phi_{1}}{dx}+\mathcal{O}\big{(}h^{2}\big{)},\qquad
k_{n+1}=\phi_{2}(y)+k\frac{d\phi_{2}}{dy}+\mathcal{O}\big{(}k^{2}\big{)},$

and thus
$a_{m,n}=-hk\frac{\phi_{1,x}\phi_{2,y}}{(\phi_{1}+\phi_{2})^{2}}+\mathcal{O}(h^{3},h^{2}k,hk^{2},k^{3})$
a first-order approximation of the general solution of (1.2) given by (1.3) .

##  5 Numerical tests of the 4-point scheme

In this section we shall apply the invariant recursion formula (2.7), (2.8) 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, respectively, 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
$\phi(x,y)=\phi_{1}(x)+\phi_{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_{\cal D}^{2}$. We compute the quantity

$\displaystyle\chi_{\alpha}(F)=\sqrt{\frac{\sum_{ij}(F_{ij}^{\alpha}-F_{ij})^{2}}{\sum_{ij}F_{ij}^{2}}},$
(5.1)

where $F_{ij}$ are the values of the exact solution $F$ on the lattice sites
and $F_{ij}^{\alpha}$, with $\alpha={\rm Inv},{\rm AS},{\rm RV}$, or ${\rm
stand}$ are the values computed numerically for the invariant, Adler–Startsev,
Rebelo–Valiquette or standard discretization, respectively. A similar analysis
is performed for the other recursion formulae. The summation will be over all
points of the lattice for which the calculation was performed.

The $\chi_{\alpha}$ quantity (5.1) provides information about the overall
averaged behavior of numerical solutions, rather then about their point-by-
point behavior. Geometrical features of solutions are better reflected by
plots of individual solutions and by a relative error function such as

$\displaystyle R_{ij}=\left|\frac{F_{ij}^{\alpha}-F_{ij}}{F_{ij}}\right|.$
(5.2)

A characteristic property of all solutions of the algebraic Liouville equation
concern the zeroes. They either have no zero in any finite domain
$\mathcal{D}$, or the zeroes are not isolated, but occur along continuous
lines parallel to the $x$ or $y$ axes.

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

$\displaystyle u_{11}u_{00}-u_{10}u_{01}=h^{2}k^{2}(uu_{xy}-u_{x}u_{y}),$

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 and their continuous limits up to
one order beyond the leading one are:

1\. The invariant method (2.7), (2.8) (preserving the ${\rm
SL}_{x}(2,\mathbb{R})\otimes{\rm SL}_{y}(2,\mathbb{R})$ symmetry group as
point symmetries)

$\displaystyle
u_{11}u_{00}-u_{10}u_{01}=hk[au_{01}u_{10}+(1-a)u_{00}u_{11}]{\,\operatorname{sign}(u_{01}u_{10})}\sqrt{|u_{01}u_{10}|},$
(5.3) $\displaystyle
uu_{xy}-u_{x}u_{y}={u}^{3}+\big{(}3{u}^{2}u_{y}+u_{yy}u_{x}+uu_{xyy}\big{)}\frac{k}{2}+\big{(}3{u}^{2}u_{x}+u_{y}u_{xx}-uu_{xxy}\big{)}\frac{h}{2}.$
(5.4)

2\. The Rebelo and Valiquette method (preserving the entire infinite-
dimensional symmetry algebra as generalized symmetries)

$\displaystyle u_{11}u_{00}-u_{10}u_{01}=hku_{00}u_{01}u_{10},$ (5.5)
$\displaystyle
uu_{xy}-u_{x}u_{y}={u}^{3}+\big{(}2{u}^{2}u_{y}+u_{x}u_{yy}-uu_{xyy}\big{)}\frac{k}{2}+\big{(}2{u}^{2}u_{x}+u_{xx}u_{y}-uu_{xxy}\big{)}\frac{h}{2}.$

3\. The Adler–Startsev method (preserving linearizability of the Liouville
equation)

$\displaystyle
u_{11}u_{00}-u_{10}u_{01}=hku_{00}u_{11}\left[\frac{u_{01}+u_{10}}{2}-hku_{01}u_{10}\right],$
(5.6) $\displaystyle
uu_{xy}-u_{x}u_{y}={u}^{3}+\big{(}2{u}^{2}u_{y}+u_{x}u_{yy}-uu_{xyy}\big{)}\frac{k}{2}+\big{(}2{u}^{2}u_{x}+u_{xx}u_{y}-uu_{xxy}\big{)}\frac{h}{2}.$
(5.7)

4\. The standard method (not preserving any specific structure) defined on the
4 points of a square lattice is

$\displaystyle u_{11}u_{00}-u_{10}u_{01}=hku_{00}^{3},$ (5.8) $\displaystyle
uu_{xy}-u_{x}u_{y}={u}^{3}+\left(u_{xx}u_{y}-uu_{xxy}\right)\frac{h}{2}+\left(u_{x}u_{yy}-uu_{xyy}\right)\frac{k}{2}.$

Each of these formulae gives a different explicit expression for $u_{11}$ in
terms of the already known values of $u_{00}$, $u_{01}$ and $u_{10}$.

Several comments are in order:

1\. In (5.4) there is no dependence on the parameter $a$. It will only appear
at the order $hk$ (not $h^{2}$ or $k^{2}$). That is the reason why the
dependence of the numerical results depends weakly on the choice of $a$ (see
Fig. 7 below to confirm this).

2\. Formulas (5.4) (invariant method) and (5.7) (Adler–Startsev method)
coincide. They differ in the higher-order terms. Table 1, as could be
expected, confirms that the results obtained by these two methods are similar.

We consider 5 different solutions of the continuous algebraic Liouville
equation (1.2), namely

$\displaystyle
f_{1}=\frac{2}{\left(x^{2}+1\right)\left(y^{2}+1\right)\left(\tan^{-1}(x)+\tan^{-1}(y)+6\right)^{2}},$
(5.9) $\displaystyle
f_{2}=\frac{8\left(1-4\left(x+\frac{1}{2}\right)\right)(1-4y)\exp\left(-2x(1+2x)-2y(y+2)\right)}{\left(e^{-2x(1+2x)}+e^{2y(1-2y)}+1\right)^{2}},$
$\displaystyle
f_{3}=-\frac{3.38\sin(1.3(x+10^{-2}))\cos(1.3(y+10^{-2}))}{(\cos(1.3(x+10^{-2}))+\sin(1.3(y+10^{-2}))+3)^{2}},$
$\displaystyle f_{4}=\frac{8xy}{\left(x^{2}+y^{2}+2\right)^{2}},$
$\displaystyle
f_{5}=\frac{383.1e^{3.862(2.5(x-0.5)+0.4y+2.5)}}{\left(e^{9.655(x+0.5)}+12.83e^{1.545y}\right)^{2}}.$

The functions $f_{1}$ and $f_{5}$ do not contain any zeroes in any finite
domain. The functions $f_{2}$ and $f_{4}$ have one row and one column of
zeroes each. The function $f_{3}$ contains infinitely many orthogonal lines of
zeroes, since it is a periodic function. Finally, $f_{5}$ is a wall like
function, with no zeroes. We mention that for $f_{1}$, $f_{3}$ and $f_{4}$ the
first-order corrections in (5.4) and (5.7) vanish. Plots of the exact
solutions $f_{1},\dots,f_{4}$ are given in Fig. 4, $f_{5}$ on Fig. 6a below.

a)  b)

c)  d)

Figure 4: Plots of the exact solutions (a) $f_{1}$, (b) $f_{2}$, (c) $f_{3}$
and (d) $f_{4}$ in the domain ${\cal D}_{1}$ used in the numerical
integrations.

The right-hand sides of (5.5), (5.6) and (5.8) are polynomials whereas the
invariant case (5.3) involves square roots.

The numerical computations were performed on the square domain ${\cal
D}_{0}=[-1.5,1.1]\times[-1.0,1.6]$, with steps of equal length $h=k=0.02$, for
a lattice of $130\times 130$ points. Somewhat arbitrarily we choose the
parameter $a$ in the symmetry invariant recursion formula to be $a=1.0$. The
boundary conditions are given on the bottom and left side of the square.

In Table 1 we give the $\chi$ quantity (5.1) for all five solutions using 4
different methods.

Table 1: Relative mean square distance (5.1) between the numerical solutions and the analytic one in the domain ${\cal D}_{0}$. | $\chi_{\rm Inv}$ | $\chi_{\rm AS}$ | $\chi_{\rm RV}$ | $\chi_{\rm stand}$
---|---|---|---|---
$f_{1}$ | $5.2\times 10^{-6}$ | $2.7\times 10^{-6}$ | $3.1\times 10^{-4}$ |  $9.2\times 10^{-4}$ 1pt
$f_{2}$ | $3.4\times 10^{-4}$ | $1.5\times 10^{-4}$ | $7.6\times 10^{-3}$ | $2.2\times 10^{-2}$
$f_{3}$ | $4.7\times 10^{-5}$ | $1.5\times 10^{-5}$ | $3.0\times 10^{-3}$ | $9.2\times 10^{-3}$
$f_{4}$ | $4.3\times 10^{-5}$ | $7.9\times 10^{-5}$ | $5.2\times 10^{-3}$ | $2.0\times 10^{-2}$
$f_{5}$ | $3.8\times 10^{-2}$ | $3.0\times 10^{-2}$ | $2.8\times 10^{-1}$ | $4.3\times 10^{-1}$

We see from Table 1 that $\chi_{\rm Inv}$ and $\chi_{\rm AS}$ are in general
of the same order, as are $\chi_{\rm RV}$ and $\chi_{\rm stand}$. The values
of $\chi_{\rm Inv}$ and $\chi_{\rm AS}$ are better than those of $\chi_{\rm
RV}$ and $\chi_{\rm stand}$ by at least one order of magnitude, usually by 2
orders, with $\chi_{\rm RV}$ always better than $\chi_{\rm stand}$. The faster
the solution changes, the worse is the result (for all methods), specially for
the solutions $f_{2}$ and $f_{5}$.

a)  b)

c)  d)

Figure 5: Plots of the relative error $R_{ij}$ defined in (5.2) for the
numerical approximation of the solution $f_{2}$ in the domain ${\cal D}_{1}$
using the individual discretizations: (a) Invariant, (b) Adler–Startsev, (c)
Rebelo–Valiquette and (d) standard. The graphs are not in the same scale. The
maximal value of the error is approximately $2\times 10^{-3}$, $2\times
10^{-4}$, $1.5\times 10^{-2}$, $4\times 10^{-2}$ for the discretizations (a),
(b), (c ) and (d), respectively.

a)  b)

c)  d)

e)

Figure 6: Pictures of the function $f_{5}$ in the domain ${\cal D}_{1}$ for
the exact solution (a) and for the individual discretizations: (b) Invariant,
(c) Adler–Startsev, (d) Rebelo–Valiquette and (e) standard. It is clear by
inspection that in this case all proposed discretizations are numerically
unstable, leading to decaying or blow up of the computed solutions by
recursive formulae.

In order to test the stability of the algorithms with respect to the size of
the adopted meshes, we made another series of calculations involving the above
test functions over a fixed domain ${\cal
D}_{1}=[-1.905,1.895]\times[-1.905,1.895]$, larger than ${\cal D}_{0}$, and
spanned it using different lattice scales with $h=k$.

Table 2: ${\rm Log}_{10}\chi$ in the domain ${\cal D}_{1}$ for each discretization procedure we considered changing the size of the lattice meshes. We report only the results for $f_{1}$ with the parameter $a=1$. $h=k$ | $\chi_{\rm Inv}$ | $\chi_{\rm AS}$ | $\chi_{\rm RV}$ | $\chi_{\rm stand}$
---|---|---|---|---
$8.\times 10^{-1}$ | $-3.98219$ | $-4.30917$ | $-2.57467$ |  $-2.13488$ 1pt
$4.\times 10^{-1}$ | $-4.5782$ | $-4.91072$ | $-2.87464$ | $-2.41647$
$2.\times 10^{-1}$ | $-5.17707$ | $-5.51254$ | $-3.17517$ | $-2.7076$
$1.\times 10^{-1}$ | $-5.77759$ | $-6.11449$ | $-3.47595$ | $-3.00363$

A general flavor of such calculations can be extracted from Table 2 for the
function $f_{1}$ For the solution $f_{1}$ (with no zeroes) the value of
$\chi_{\rm Inv}$ and $\chi_{\rm AS}$ are comparable and at least two orders of
magnitude lower than the other two. The values of $\chi_{\rm RV}$ are always
lower than $\chi_{\rm stand}$ but of the same order. Generally speaking,
decreasing the scale of the mesh by a factor of $0.3$ implies decreasing the
value of $\chi$ by a factor of $0.6$ for the invariant and Adler–Startsev
discretization and by a factor of $0.3$ for the other two discretization
(quadratic as opposed to linear convergence).

Table 3: Value of $f_{2}$ at the four lattice points nearest to the saddle point ($-0.25,+0.25$) for the exact solution and the four numerical approximations. point | Exact | Inv | AS | RV | stand
---|---|---|---|---|---
($-0.255,0.255$) | $-4.14391\cdot 10^{-4}$ | $-4.1419\cdot 10^{-4}$ | $-4.14369\cdot 10^{-4}$ | $-4.1281\cdot 10^{-4}$ |  $-4.11579\cdot 10^{-4}$ 1pt
($-0.245,0.255$) | $4.14391\cdot 10^{-4}$ | $4.14257\cdot 10^{-4}$ | $4.14369\cdot 10^{-4}$ | $4.12877\cdot 10^{-4}$ | $4.11513\cdot 10^{-4}$
($-0.255,0.245$) | $4.14391\cdot 10^{-4}$ | $4.14123\cdot 10^{-4}$ | $4.14369\cdot 10^{-4}$ | $4.12744\cdot 10^{-4}$ | $4.11645\cdot 10^{-4}$
($-0.245,0.245$) | $-4.14391\cdot 10^{-4}$ | $-4.1419\cdot 10^{-4}$ | $-4.14369\cdot 10^{-4}$ | $-4.1281\cdot 10^{-4}$ | $-4.11579\cdot 10^{-4}$

Let us now turn to the point-by-point behavior of the solutions. As discussed
in Section 2 a characteristic property of all solutions of the Liouville
equation is that zeroes are not isolated but occur in straight lines parallel
to the axes. To see how well this is reflected in numerical solutions, let us
first concentrate on solution $f_{2}$ which has zeroes on horizontal and
vertical lines passing through the saddle point $x=-0.25$, $y=+0.25$. These
points are not on the lattice due to the definition of the domain
$\mathcal{D}_{1}$ with $h=k=0.01$. On Table 3 we give the values of the
solution $f_{2}$ at the four lattice points nearest to the saddle point. The
AS solution has the first 4 digits coinciding with the exact one, the
invariant one has 3. The RV and the standard method have 2. Similar results
are also valid for the other solutions with zeroes ($f_{3}$ and $f_{4}$).

The behavior of the solution $f_{2}$ is plotted on Fig. 5 for the entire
region $\mathcal{D}_{1}$ where we show the values of the error function
$R_{i,j}$ of (5.2). The maximal value of $R_{i,j}$ is $2\times 10^{-4}$ in
Fig. 5b for the AS method, 10 times higher for the invariant method, 75 higher
for RV and 200 times higher for the standard one. For the invariant method the
error is concentrated at the saddle point (see Fig. 5a) with a tail behind it.
The AS method has maximal error at the four extrema (see Fig. 5b) with no
tail. The other two methods have maximal errors on the maxima of the solutions
(not the minima) with tails in both directions.

On Fig. 6 we analyze the solution $f_{5}$ in detail in a point-by-point manner
in the domain $\mathcal{D}_{1}$. Fig. 6a represents the exact solution (a wall
of constant height). It has no zeroes anywhere in the finite real plane. The
height of the wall gradually decreases for solutions $b$, $d$ and $e$, but
much more slowly for the invariant method $b$. For the AS method the height
increases so we present the height of the solution on Fig. 6c in a different
scale. The increase in Fig. 6c is comparable with the decrease in Fig. 6b (a
factor of about 2.5).

Finally we analyze the role of the parameter $a$ in the formula (2.7), (2.8).
The continuous limit (5.4) shows that $a$ appears for the first time in terms
of order $O(hk)$. Numerical calculations of $\chi_{\rm Inv}$ for the function
$f_{2}$ shows a variation of about a factor 2, when $a\in[-0.5,1.5]$ (see Fig.
7). Furthermore, this function takes a minimum for $a=0.17$. However, this
value is strongly dependent on the test function considered.

Figure 7: Plot of $\chi_{\rm Inv}$ for the function $f_{2}$ in the parameter
range $a\in[-0.5,1.5]$.

##  6 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, both using
the global $\chi^{2}$ criterion and comparing local point by point convergence
using the $R_{ij}$ criterion. 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 [38] 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 point
of view of we have shown that that the precision of the RV solutions is
systematically better than that of those obtained by the standard method
(though of the same order of magnitude). The measure of the validity is given
by the quantity $\chi$ of (5.1).

Finally, the method proposed in [23] 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 very well for all solutions in some case
even better than the Adler–Startsev case.

In future work we plan to study symmetry preserving discretizations of other
equations with infinite-dimensional symmetry groups, such as the
Kadomtsev–Petviashvili equation, and the three-wave interaction equation. A
symmetry preserving discretization of the Korteweg–de Vries equation has
provided encouraging results [3].

### Acknowledgments

DL has been partly supported by the Italian Ministry of Education and
Research, 2010 PRIN Continuous and discrete nonlinear integrable evolutions:
from water waves to symplectic maps.

LM has been partly supported by the Italian Ministry of Education and
Research, 2011 PRIN Teorie geometriche e analitiche dei sistemi Hamiltoniani
in dimensioni finite e infinite. DL and LM are supported also by INFN IS-CSN4
Mathematical Methods of Nonlinear Physics.

The research of PW is partially supported by a research grant from NSERC of
Canada. PW thanks the European Union Research Executive Agency for the award
of a Marie Curie International Incoming Research Fellowship making his stay at
University Roma Tre possible. He thanks the Department of Mathematics and
Physics of Roma Tre for hospitality.

We thank the referees for many valuable comments which allowed us to greatly
improve the article.

## References

  * [1] Adler V.E., Startsev S.Ya., On discrete analogues of the Liouville equation, Theoret. and Math. Phys. 121 (1999), 1484–1495, arXiv:solv-int/9902016.
  * [2] Bihlo A., Invariant meshless discretization schemes, J. Phys. A: Math. Theor. 46 (2013), 062001, 12 pages, arXiv:1210.2762.
  * [3] Bihlo A., Coiteux-Roy X., Winternitz P., The Korteweg–de Vries equation and its symmetry-preserving discretization, J. Phys. A: Math. Theor. 48 (2015), 055201, 25 pages, arXiv:1409.4340.
  * [4] Bihlo A., Nave J.-C., Invariant discretization schemes using evolution-projection techniques, SIGMA 9 (2013), 052, 23 pages, arXiv:1209.5028.
  * [5] Bourlioux A., Cyr-Gagnon C., Winternitz P., Difference schemes with point symmetries and their numerical tests, J. Phys. A: Math. Gen. 39 (2006), 6877–6896, math-ph/0602057.
  * [6] Budd C., Dorodnitsyn V., Symmetry-adapted moving mesh schemes for the nonlinear Schrödinger equation, J. Phys. A: Math. Gen. 34 (2001), 10387–10400.
  * [7] Dorodnitsyn V., Transformation groups in difference spaces, J. Sov. Math. 55 (1991), 1490–1517.
  * [8] Dorodnitsyn V., Applications of Lie groups to difference equations, Differential and Integral Equations and Their Applications, Vol. 8, CRC Press, Boca Raton, FL, 2011.
  * [9] Dorodnitsyn V., Kaptsov E., Kozlov R., Winternitz P., The adjoint equation method for constructing first integrals of difference equations, J. Phys. A: Math. Theor. 48 (2015), 055202, 32 pages, arXiv:1311.1597.
  * [10] Dorodnitsyn V., Kozlov R., A heat transfer with a source: the complete set of invariant difference schemes, J. Nonlinear Math. Phys. 10 (2003), 16–50, math.AP/0309139.
  * [11] Dorodnitsyn V., Kozlov R., Winternitz P., Lie group classification of second-order ordinary difference equations, J. Math. Phys. 41 (2000), 480–504.
  * [12] Dorodnitsyn V., Kozlov R., Winternitz P., Continuous symmetries of Lagrangians and exact solutions of discrete equations, J. Math. Phys. 45 (2004), 336–359, nlin.SI/0307042.
  * [13] Dubrovin B.A., Fomenko A.T., Novikov S.P., Modern geometry – methods and applications. Part I. The geometry of surfaces, transformation groups, and fields, Graduate Texts in Mathematics, Vol. 93, 2nd ed., Springer-Verlag, New York, 1992.
  * [14] Floreanini R., Negro J., Nieto L.M., Vinet L., Symmetries of the heat equation on the lattice, Lett. Math. Phys. 36 (1996), 351–355.
  * [15] Floreanini R., Vinet L., Lie symmetries of finite-difference equations, J. Math. Phys. 36 (1995), 7024–7042.
  * [16] Floreanini R., Vinet L., Quantum symmetries of $q$-difference equations, J. Math. Phys. 36 (1995), 3134–3156.
  * [17] Grammaticos B., Ramani A., Painlevé equations, continuous, discrete and ultradiscrete, in Symmetries and Integrability of Difference Equations (Université de Montréal, Montréal, QC, June 8–21, 2008), London Mathematical Society Lecture Note Series, Vol. 381, Editors D. Levi, P.J. Olver, Z. Thomova, P. Winternitz, Cambridge University Press, Cambridge, 2011, 50–82.
  * [18] Grammaticos B., Ramani A., Papageorgiou V., Do integrable mappings have the Painlevé property?, Phys. Rev. Lett. 67 (1991), 1825–1828.
  * [19] Grant T.J., Bespoke finite difference schemes that preserve multiple conservation laws, LMS J. Comput. Math. 18 (2015), 372–403.
  * [20] Hairer E., Lubich C., Wanner G., Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations, Springer Series in Computational Mathematics, Vol. 31, 2nd ed., Springer-Verlag, Berlin, 2006.
  * [21] Hydon P.E., Difference equations by differential equation methods, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, 2014.
  * [22] Iserles A., A first course in the numerical analysis of differential equations, 2nd ed., Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 2009.
  * [23] Levi D., Martina L., Winternitz P., Lie-point symmetries of the discrete Liouville equation, J. Phys. A: Math. Theor. 48 (2015), 025204, 18 pages, arXiv:1407.4043.
  * [24] Levi D., Olver P.J., Thomova Z., Winternitz P. (Editors), Symmetries and integrability of difference equations, London Mathematical Society Lecture Note Series, Vol. 381, Cambridge University Press, Cambridge, 2011.
  * [25] Levi D., Rodríguez M.A., Construction of partial difference schemes: I. The Clairaut, Schwarz, Young theorem on the lattice, J. Phys. A: Math. Theor. 46 (2013), 295203, 19 pages.
  * [26] Levi D., Rodríguez M.A., On the construction of partial difference schemes: II. Discrete variables and invariant schemes, arXiv:1407.0838.
  * [27] Levi D., Vinet L., Winternitz P., Lie group formalism for difference equations, J. Phys. A: Math. Gen. 30 (1997), 633–649.
  * [28] Levi D., Winternitz P., Continuous symmetries of discrete equations, Phys. Lett. A 152 (1991), 335–338.
  * [29] Levi D., Winternitz P., Continuous symmetries of difference equations, J. Phys. A: Math. Gen. 39 (2006), R1–R63, nlin.SI/0502004.
  * [30] Levi D., Winternitz P., Yamilov R.I., Lie point symmetries of differential-difference equations, J. Phys. A: Math. Theor. 43 (2010), 292002, 14 pages, arXiv:1004.5311.
  * [31] Lie S., Theorie der Transformationsgruppen, Bds. 1–3, B.G. Teubner, Leipzig, 1888, 1890, 1893.
  * [32] Liouville J., Sur l’equation aux différences partielles $\frac{d^{2}\log\lambda}{dudv}\pm\frac{\lambda}{2a^{2}}=0$, J. Math. Pures Appl. 18 (1853), 71–72.
  * [33] Marsden J.E., West M., Discrete mechanics and variational integrators, Acta Numer. 10 (2001), 357–514.
  * [34] McLachlan R.I., Quispel G.R.W., Geometric integrators for ODEs, J. Phys. A: Math. Gen. 39 (2006), 5251–5285.
  * [35] Medolaghi P., Classificazione delle equazioni alle derivate parziali del secondo ordine, che ammettono un gruppo infinito di trasformazioni puntuali, Ann. Mat. Pura Appl. 1 (1898), 229–263.
  * [36] Olver P.J., Applications of Lie groups to differential equations, Graduate Texts in Mathematics, Vol. 107, 2nd ed., Springer-Verlag, New York, 1993.
  * [37] Rebelo R., Valiquette F., Symmetry preserving numerical schemes for partial differential equations and their numerical tests, J. Difference Equ. Appl. 19 (2013), 738–757, arXiv:1110.5921.
  * [38] Rebelo R., Valiquette F., Invariant discretization of partial differential equations admitting infinite-dimensional symmetry groups, J. Difference Equ. Appl. 21 (2015), 285–318, arXiv:1401.4380.
  * [39] Rebelo R., Winternitz P., Invariant difference schemes and their application to ${\rm sl}(2,{\mathbb{R}})$ invariant ordinary differential equations, J. Phys. A: Math. Theor. 42 (2009), 454016, 10 pages, arXiv:0906.2980.
  * [40] Rodríguez M.A., Winternitz P., Lie symmetries and exact solutions of first-order difference schemes, J. Phys. A: Math. Gen. 37 (2004), 6129–6142, nlin.SI/0402047.
  * [41] Valiquette F., Winternitz P., Discretization of partial differential equations preserving their physical symmetries, J. Phys. A: Math. Gen. 38 (2005), 9765–9783, math-ph/0507061.
  * [42] Winternitz P., Symmetry preserving discretization of differential equations and Lie point symmetries of differential-difference equations, in Symmetries and Integrability of Difference Equations (Université de Montréal, Montréal, QC, June 8–21, 2008), London Mathematical Society Lecture Note Series, Vol. 381, Editors D. Levi, P.J. Olver, Z. Thomova, P. Winternitz, Cambridge University Press, Cambridge, 2011, 292–341.