∎

11institutetext: C. C. da Silva, C. Lessig 22institutetext: Otto-von-Guericke-
Universität Magdeburg, Institut für Simulation und Graphik
22email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Variational symplectic diagonally implicit Runge-Kutta methods for
isospectral systems††thanks: Funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – Project-ID 422037413 – TRR 287.

Clauson Carvalho da Silva  Christian Lessig

(Received: date / Accepted: date)

###### Abstract

Isospectral flows appear in a variety of applications, e.g. the Toda lattice
in solid state physics or in discrete models for two-dimensional
hydrodynamics, with the isospectral property often corresponding to
mathematically or physically important conservation laws. Their most prominent
feature, i.e. the conservation of the eigenvalues of the matrix state
variable, should therefore be retained when discretizing these systems.
Recently, it was shown how isospectral Runge-Kutta methods can, in the Lie-
Poisson case also considered in our work, be obtained through Hamiltonian
reduction of symplectic Runge-Kutta methods on the cotangent bundle of a Lie
group. We provide the Lagrangian analogue and, in the case of symplectic
diagonal implicit Runge-Kutta methods, derive the methods through a discrete
Euler-Poincaré reduction. Our derivation relies on a formulation of diagonally
implicit isospectral Runge-Kutta methods in terms of the Cayley transform,
generalizing earlier work that showed this for the implicit midpoint rule. Our
work is also a generalization of earlier variational Lie group integrators
that, interestingly, appear when these are interpreted as update equations for
intermediate time points. From a practical point of view, our results allow
for a simple implementation of higher order isospectral methods and we
demonstrate this with numerical experiments where both the isospectral
property and energy are conserved to high accuracy.

###### Keywords:

Runge Kutta methods Lie group integrators variational integrators

###### MSC:

MSC 65L06 MSC 65P10

††journal: BIT

##  1 Introduction

Isospectral systems, for which the time evolution preserves the spectrum of a
matrix state variable, arise in the mathematical description of various
systems of interest. They are generically given by

$\displaystyle\dot{A}=\big{[}B(A),A\big{]},\quad A(0)=A_{0}$ (1)

where $A$, $B(A)$ are matrices and the spectrum of $A$ is preserved
(hairer2006geometric, 8, IV.3.2). See calvo1997numerical (5) for a detailed
discussion of the properties of these systems.

A classical example for an isospectral system is the ideal rigid body in its
reduced formulation Marsden1999 (11), as described by Euler’s equations. In
this case the angular momentum (as an anti-symmetric matrix in
$\mathfrak{so}^{*}(3)$) is the isospectrally conserved quantity. The periodic
Toda lattice in its Lax pair formulation is also isospectral with the
isospectral state $L$ in this case being symmetric and the generator of the
dynamics $B(L)$ being anti-symmetric. Another example is the discrete 2D
barotropic vorticity equation on the torus and 2-sphere when the “matrix
model”, i.e. the representation theory of the discrete Heisenberg group
respectively $\mathrm{SO}(3)$, is used for spatial discretization Zeitlin1991
(20, 9, 21). The isospectral property is in this case a discrete analogue
(with $N$ conserved quantities for $N$ degrees of freedom) of the integrated
powers of vorticity that are conserved in the continuous system Khesin1998
(1).

Since conserved quantities are an important intrinsic characteristic of a
physical system, it is natural to preserve isospectrality under
discretization. An early work that accomplished this was those by Moser and
Veselov Moser1991 (16) for rigid body like systems. Moore, Mahony, and Helmke
Moore1994 (15) introduced isospectral time integrators for the solution of
problems in linear algebra, since many classical ones there, such as the QR or
eigen decomposition, can be cast as isospectral flows. Isospectral numerical
time integrators were also considered by Diele, Lopez and Politi who proposed
the use of the Cayley transform Diele1998 (6). Calvo, Iserles and Zanna
calvo1997numerical (5) showed that Runge-Kutta methods applied to Eq. 1 break
isospectrality. To overcome this, they introduce modified Gauss-Legendre
Runge-Kutta methods that provide isospectral time integration schemes of
arbitrary high order. Recently, Modin and Vivani modin2019lie (13) developed
an alternative route to overcome the obstruction from calvo1997numerical (5)
by applying a Runge-Kutta scheme to the cotangent bundle $T^{*}G$ of a
quadratic matrix Lie group $G$ and then reducing to the dual Lie algebra
$\mathfrak{g}^{*}$ using the momentum map, analogous to how Lie-Poisson
reduction works in the continuous case, cf. (Marsden1999, 11, Ch. 13). The
ansatz exploits that many isospectral flows are Lie-Poisson, i.e. their phase
space is a dual Lie algebra $\mathfrak{g}^{*}$, and the dynamics are reduced
ones from $T^{*}G$. Vivani viviani2019minimal (19) recently extended
modin2019lie (13) and showed that the implicit midpoint (Runge-Kutta) scheme
can be written and implemented efficiently with an intermediate time step (in
the spirit of the Verlet method). He also pointed out that the Cayley
transform then arises naturally and does not have to be introduced as a group
map as in Diele1998 (6).

Another line of research on discretizations that retain conserved quantities
is based on the Lagrangian formulation of mechanics and uses a discrete
Hamilton’s principle marsden2001discrete (12). This ansatz has also been
developed for reduced dynamical systems on Lie groups and algebras
BouRabee2007 (4, 3). Gawlik et al. gawlik2011geometric (7), in particular,
developed a formulation for semi-direct products that, as a special case, is
applicable to systems such as the rigid body.

In the present work, we will derive the Lagrangian analogue of the work of
Modin and Vivani modin2019lie (13) and show that isospectral Runge-Kutta
methods can be obtained from the reduced Hamilton’s principle of Bou-Rabee and
Marsden BouRabee2009 (3) and Gawlik et al. gawlik2011geometric (7). We do so
for the case of symplectic diagonally implicit Runge-Kutta methods (SDIRKs),
which are an important special case since they are the simplest and most
easily implemented symplectic Runge-Kutta schemes. Similar to modin2019lie
(13), we will consider isospectral Lie-Poisson systems for our derivation. We
also extend the results of Vivani viviani2019minimal (19) by showing that the
Cayley transform arises for an SDIRK applied at the group level when it is
written using intermediate time steps located between the Runge-Kutta steps.
Another result of the present work is that the midpoint rule by Gawlik et al.
gawlik2011geometric (7) contains, in disguise, an isospectral method; more
precisely, the scheme from gawlik2011geometric (7) updates the midpoints of
the isospectral midpoint rule of Vivani viviani2019minimal (19). We exemplify
our results by implementing our isospectral SDIRKs for the rigid body, the
Toda lattice, and the matrix model discretization of the 2D Euler fluid. We
show that with a $7$-stage $4^{\textrm{th}}$-order SDIRK, the oscillations in
the energy, which are standard for a symplectic integrator, are on the order
of machine precision. Within the accuracy afforded by standard double
precision, the integrator therefore preserves the invariants of the continuous
system. SDIRKs provide here the advantage that these can be implemented more
easily and can be computed more efficiently than general isospectral Runge-
Kutta methods.

The remained of the paper is structured as follows. In Sec. 2 we will recall
the necessary background on Lie groups, their algebras and isospectral Lie-
Poisson systems. In Sec. 3 we will show that the Cayley transform arises for
any SDIRK on a quadratic matrix Lie group by introducing intermediate time
steps. The central result of the paper, namely a variational derivation of
isospectral SDIRKs, will be presented in Sec. 4. Numerical results for the
rigid body, the Toda lattice, and the discrete 2D Euler fluid are discussed in
Sec. 5.

##  2 Preliminaries

###  2.1 Lie groups, Lie algebras and retraction maps

Let $G$ be a matrix Lie group with identity $e$ and $\mathfrak{g}$ its Lie
algebra. We will denote the conjugate transpose of an element $g\in G$ with
$g^{\dagger}$. If conjugation is defined as a left action $g\cdot
h=g\,h\,g^{-1}$, the adjoint action
$\mathrm{Ad}_{g}:G\times\mathfrak{g}\to\mathfrak{g}$ of the group $G$ on its
Lie algebra $\mathfrak{g}$ is given by $\text{Ad}_{g}\cdot\xi=g\,\xi\,g^{-1}$.
Infinitesimally it becomes $\mathrm{ad}_{\xi}(\eta)=[\xi,\eta]$. The dual of
$\mathrm{ad}$ is the infinitesimal coadjoint action
$\text{ad}^{*}:\mathfrak{g}\times\mathfrak{g}^{*}\to\mathfrak{g}^{*}$ defined
through
$\smash{\langle\mathrm{ad}_{\xi}(\eta),\alpha\rangle=\langle\eta,\mathrm{ad}_{\xi}^{*}(\alpha)\rangle}$,
which it plays an important role for Lie-Poisson dynamics. Unless mentioned
otherwise, we will work with the pairing between $\mathfrak{g}$ and
$\mathfrak{g}^{*}$ associated with the Frobenius inner-product on
$\mathfrak{g}$ given by
$\langle\xi,\alpha\rangle=\langle\xi,\alpha^{\sharp}\rangle=\mathrm{tr}(\xi^{\dagger}\alpha^{\sharp})$
which also allows us to identify $\mathfrak{g}$ and $\mathfrak{g}^{*}$. We
ommit the sharp symbol $\sharp$ throughout the paper. A local diffeomorphism
$\tau:\mathfrak{g}\to G$ satisfying $\tau(0)=e$, $\tau(\xi)^{-1}=\tau(-\xi)$
and whose tangent map fulfills $T_{0}\tau=\text{Id}$ is called a retraction
map (or local group diffeomorphism). Such a map provides a local chart for the
Lie group around the identity (with coordinates of the first kind) and induces
an atlas for $G$ through the translation action. The tangent $T_{\xi}\tau$ of
$\tau$ maps $T_{\xi}\mathfrak{g}\cong T_{e}G=\mathfrak{g}$ to
$T_{\tau(\xi)}G$. The same is accomplished with the tangent map
$T_{e}R_{\tau(\xi)}$ where $R_{g}$ is right translation $R_{g}(h)=h\cdot g$.
The maps $T_{\xi}\mathfrak{g}$ and $T_{e}R_{\tau(\xi)}$ are related by the
right trivialization $\mathrm{d}\tau_{\xi}:\mathfrak{g}\to\mathfrak{g}$ of
$\tau$,

$\displaystyle T_{\xi}\tau=T_{e}R_{\tau(\xi)}\circ\mathrm{d}\tau_{\xi}.$

Since $\tau$ is a local diffeomorphism, the map $T_{\xi}\tau$ is invertible
for $\xi$ close enough to $0\in\mathfrak{g}$ which implies that
$\mathrm{d}\tau_{\xi}$ is also invertible. We will write
$\mathrm{d}\tau^{-1}_{\xi}\equiv(\mathrm{d}\tau_{\xi})^{-1}:\mathfrak{g}\to\mathfrak{g}$.
A basic property of the right trivialization used throughout the paper is
$\mathrm{d}\tau_{\xi}=\text{Ad}_{\tau(\xi)}\circ\mathrm{d}\tau^{-1}_{-\xi}$,
cf. Owren2016 (17).

In the following, we will consider quadratic matrix Lie groups, i.e. we will
work with $G\subset\mathrm{GL}(N)$ so that there exists a $J$ such that
$J^{\dagger}\,g\,J=g,\forall g$. For $\xi\in\mathfrak{g}$ it then holds that
$J\,\xi+\xi^{\dagger}J=0$. Examples of quadratic Lie groups are
$\mathrm{O}(n)$ and $\mathrm{SO}(n)$ with $J=\mathrm{Id}$ and the symplectic
group $\mathrm{Sp}(2N,\mathbb{R})$ with $J$ the symplectic matrix. An
expedient retraction map for quadratic Lie groups is the Cayley transform
BouRabee2007 (4)

$\displaystyle\mathrm{cay}(\xi)=\left(\text{Id}-\frac{\xi}{2}\right)^{-1}\\!\left(\text{Id}+\frac{\xi}{2}\right).$
(2)

The two factors in the definition above commute. The right trivialization
$\mathrm{d}\tau$ for the Cayley transform is given by hairer2006geometric (8)

$\displaystyle\mathrm{d}\mathrm{cay}_{\xi}\cdot\eta$
$\displaystyle=\left(\text{Id}-\frac{\xi}{2}\right)^{-1}\eta\left(\text{Id}+\frac{\xi}{2}\right)^{-1}$
$\displaystyle\mathrm{d}\mathrm{cay}^{-1}_{\xi}\cdot\eta$
$\displaystyle=\left(\text{Id}-\frac{\xi}{2}\right)\eta\left(\text{Id}+\frac{\xi}{2}\right).$

The configuration space considered in the paper are quadratic Lie groups.
Therefore the Cayley transform plays a major role throughout the paper.

###  2.2 Lie-Poisson systems on reductive Lie algebras

We recall Lie-Poisson systems on the dual of a Lie algebra
$\mathfrak{g}\subseteq\mathfrak{gl}(n)$ and the associated reconstruction
process to the cotangent bundle of the Lie group $T^{*}G$. The flow of a Lie-
Poisson system on a dual Lie algebra $\mathfrak{g}^{*}$ is a solution to

$\displaystyle\dot{\mu}=\text{ad}^{*}_{\nabla
H(\mu)}\mu,\quad\mu\in\mathfrak{g}^{*}$ (3)

where $H:\mathfrak{g}^{*}\to\mathbb{R}$ is a Hamiltonian function. Using the
Frobenius inner product
$\langle\xi,\eta\rangle=\mathrm{tr}(\xi^{\dagger}\eta)$ in $\mathfrak{g}$, we
can identify $\mathfrak{g}$ with $\mathfrak{g}^{*}$ and the dual of the
infinitesimal adjoint action then becomes
$\mathrm{ad}^{*}_{\xi}(\eta)=\mathrm{\Pi}([\xi^{\dagger},\eta])$ where
$\mathrm{\Pi}$ is the projection from $\mathfrak{gl}(n)$ onto $\mathfrak{g}$.
For the case where $\mathfrak{g}$ is a reductive Lie algebra, it is known that
$[\mathfrak{g}^{\dagger},\mathfrak{g}]\subset\mathfrak{g}$ and therefore the
Lie-Poisson equation 3 becomes

$\displaystyle\dot{\mu}=\left[\nabla H(\mu)^{\dagger},\mu\right].$ (4)

An important property of such systems is isospectrality, i.e, the eigenvalues
of the matrices $\mu(t)$ are independent of $t$.

Every Lie-Poisson system can be reconstructed to a canonical Hamiltonian
system on $T^{*}G$. The main tool in Lie-Poisson reconstruction is the
momentum map $\mathcal{J}(g,p):T^{*}G\to\mathfrak{g}^{*}$ which for the right
action $(g,p)\cdot h=(gh,p(h^{-1})^{\dagger})$ of the group
$G\subseteq\mathrm{GL}(n)$ on the phase space $T^{*}G$, is given by
$\mu=\mathcal{J}(g,p)=g^{\dagger}p$. With $\mathcal{J}$ we can define a left
invariant Hamiltonian $\tilde{H}(g,p)=H(\mathcal{J}(g,p))=H(g^{\dagger}p)$ and
the canonical Hamilton equations on $T^{*}G$ are

$\displaystyle\begin{aligned} \dot{g}&=g\nabla
H(g^{\dagger}p)=gB(g^{\dagger}p)^{\dagger}\\\\[4.0pt] \dot{p}&=-p\nabla
H(g^{\dagger}p)^{\dagger}=-pB(g^{\dagger}p)\end{aligned}$ (5)

where we write $B(\cdot)=\nabla
H(\cdot)^{\dagger}:\mathfrak{g}\to\mathfrak{g}$ as in modin2019lie (13). For
the right invariant version of the Lie-Poisson reconstruction, see Appendix A.

###  2.3 Symplectic diagonally implicit Runge-Kutta methods

The Butcher tableau for an $s$-stage symplectic diagonally implicit Runge-
Kutta methods (SDIRK) is (hairer2006geometric, 8, Ch. V.2.1)

$\displaystyle\begin{array}[]{c|cccc}c_{1}&\frac{b_{1}}{2}&\\\\[3.0pt]
c_{2}&b_{1}&\frac{b_{2}}{2}&&\\\\[3.0pt]
\vdots&\vdots&\vdots&\ddots\\\\[3.0pt]
c_{s}&b_{1}&b_{2}&\cdots&\frac{b_{s}}{2}\\\\[3.0pt]
\hline\cr\rule{0.0pt}{12.11993pt}&b_{1}&b_{2}&\cdots&b_{s}\\\\[3.0pt]
\end{array}$ (11)

Given a fixed time step $h$, we will write $h_{i}=h\,b_{i}$ and
$r_{i}=\sum_{j=1}^{i}b_{j}$ for $i=0,\ldots,s$ so that, through the diagonal
entries in Eq. 11, $c_{i}$ and $r_{i-1}$ always differ by $b_{i}/2$. Hence
$\sum^{s}_{i}h_{i}=h$ and $r_{s}=1$ and also $b_{0}=0$. An important result
for us is that an $s$-stage SDIRK can be decomposed in $s$ implicit midpoint
rule steps with suitable step size (hairer2006geometric, 8, Theorem VI.4.4).

##  3 The Cayley transform and symplectic diagonally implicit Runge-Kutta
methods

As discussed in the last section, any $s$-stage symplectic diagonally implicit
Runge-Kutta method (SDIRK) can be decomposed into $s$ implicit midpoint rule
steps (hairer2006geometric, 8, Theorem VI.4.4). In this section we will show
that this carries over to isospectral SDIRKs and hence the simplicity of the
isospectral minimal midpoint algorithm presented in viviani2019minimal (19)
also holds for higher order isospectral SDIRKs. We will also show that the
Cayley transform appears naturally in this context. The key to this result is
to introduce not only the usual intermediate time steps that appear in a
Runge-Kutta method, i.e. $k_{n,i}^{g}$ and $k_{n,i}^{p}$ in the Hamiltonian
setting, but also intermediate half time steps at
$t_{i}^{n}=hn+h\sum_{j=1}b_{j}$. These intermediate half points, denoted
$g_{n+r_{i}}$ and $p_{n+r_{i}}$, can be seen as a generalization of the
intermediate step that one considers for the Verlet integrator and they are
defined through the strictly lower diagonal part of the Butcher tableau in Eq.
11.

Although the Lie groups we consider are nonlinear, since they are quadratic it
nonetheless holds that a symplectic Runge-Kutta method applied linearly to
Hamilton’s equations in Eq. 5 yields $g_{n}\in G,\forall n$. Exploiting this
for an $s$-stage SDIRK, the Runge-Kutta steps can be written in terms of the
intermediate Runge-Kutta points as

$\displaystyle k_{n,i}^{g}$
$\displaystyle=\left(g_{n,r_{i-1}}+\frac{h_{i}}{2}\,k_{n,i}^{g}\right)B({\mu}_{n,c_{i}})^{\dagger}$
(12a) $\displaystyle k_{n,i}^{p}$
$\displaystyle=-\left(p_{n,r_{i-1}}+\frac{h_{i}}{2}\,k_{n,i}^{p}\right)B({\mu}_{n,c_{i}}).$
(12b) where the reduced momentum is
$\displaystyle{\mu}_{n,c_{i}}=\left(g_{n,r_{i-1}}+\frac{h_{i}}{2}\,k_{n,i}^{g}\right)^{\\!\dagger}\\!\left(p_{n,r_{i-1}}+\frac{h_{i}}{2}\,k_{n,i}^{p}\right)$
(12c)

by the definition of the momentum map. The intermediate points $g_{n+r_{i}}$
and $p_{n+r_{i}}$ are defined as

$\displaystyle g_{n,r_{i}}=g_{n}+\sum_{j=0}^{i}h_{j}k_{n,j}^{g}\quad\quad
p_{n,r_{i}}=p_{n}+\sum_{j=0}^{i}h_{j}k_{n,j}^{p}$

and they are always a half step away from $k_{n,i+1}^{g}$ and $k_{n,i+1}^{p}$,
cf. Eq. 11. Here we have $g_{n,r_{0}}=g_{n}$ and $g_{n,r_{s}}=g_{n+1}$ (and
analogous for the variable $p$).

The first intermediate point at $t=hn+b_{1}$ is therefore

$\displaystyle
g_{n,r_{1}}=g_{n}+h_{1}\,k_{n,1}^{g}=g_{n}+h_{1}\,\left(g_{n}+\frac{h_{1}}{2}k_{n,1}^{g}\right)B({\mu}_{n,c_{1}})^{\dagger}.$

Fully expanding the right hand side of the last equation and using from the
left hand side that $k_{1}^{g}=(g_{n,r_{1}}-g_{n})/h_{1}$ we obtain

$\displaystyle
g_{n,r_{1}}\left(\text{Id}-\frac{h_{1}}{2}B({\mu}_{n,c_{1}})^{\dagger}\right)=g_{n}\left(\text{Id}+\frac{h_{1}}{2}B({\mu}_{n,c_{1}})^{\dagger}\right).$

The last equation states conceptually that taking a (linear) half step
backward from $g_{n,r_{1}}$ and a (linear) half step forward from $g_{n}$
agree. With the definition of the Cayley transform in Eq. 2 and taking into
account that the terms in the parenthesis commute, this is equivalently given
by

$\displaystyle
g_{n,r_{1}}=g_{n}\,\mathrm{cay}(h_{1}\,B({\mu}_{n,c_{1}})^{\dagger}).$

Following the same steps as above for $g_{n,r_{i}}$ one obtains
$g_{n,r_{i}}=g_{n,{r_{i-1}}}\,\mathrm{cay}(h_{i}\,B({\mu}_{n,c_{i}})^{\dagger})$.
and we conclude that

$\displaystyle
g_{n+1}=g_{n}\prod_{i=1}^{s}\mathrm{cay}(h_{i}\,B({\mu}_{n,c_{i}})^{\dagger}).$
(13)

An analogous calculation shows that

$\displaystyle
p_{n+1}=p_{n}\prod_{i=1}^{s}\mathrm{cay}(h_{i}\,B({\mu}_{n,c_{i}}))^{-1}.$
(14)

With the momentum map $\mathcal{J}$, the phase space points $(q_{n},p_{n})$
are related to a sequence of reduced momenta $\mu_{n}\in\mathfrak{g}^{*}$ by
$\mu_{n}=\mathcal{J}(g_{n},p_{n})=g_{n}^{\dagger}p_{n}$. Using Eq. 13 and Eq.
14 we thus have

$\displaystyle\mu_{n+1}=\left(\prod_{i=s}^{1}\mathrm{cay}(h_{i}\,B({\mu}_{n,c_{i}}))\right)\,\mu_{n}\left(\prod_{i=1}^{s}\mathrm{cay}(h_{i}\,B({\mu}_{n,c_{i}}))^{-1}\right).$
(15)

Using the momentum map to define the reduced momentum variable also at the
intermediate point we have

$\displaystyle\mu_{n,r_{i}}=g_{n,r_{i}}^{\dagger}p_{n,r_{i}}=\mathrm{cay}(h_{i}\,B({\mu}_{n,c_{i}}))\mu_{n,r_{i-1}}\mathrm{cay}(h_{i}\,B({\mu}_{n,c_{i}}))^{-1}.$
(16)

and a calculation shows that

$\displaystyle\mu_{n,c_{i}}=\mathrm{d}\mathrm{cay}_{h_{i}\,B({\mu}_{n,c_{i}}))}\,\mu_{n,r_{i-1}}$

Eq. 16 can be written as

$\displaystyle\mu_{n,r_{i}}=\mathrm{Ad}_{\mathrm{cay}(h_{i}\,B({\mu}_{n,c_{i}}))}\,\mu_{n,r_{i-1}}=\mathrm{d}\mathrm{cay}^{-1}_{-h_{i}\,B({\mu}_{n,c_{i}}))}\,\circ\mathrm{d}\mathrm{cay}_{h_{i}\,B({\mu}_{n,c_{i}}))}\mu_{n,r_{i-1}}.$

One thus has the following algorithm for updating the intermediate algebra
points

$\displaystyle\mu_{n,r_{i-1}}$
$\displaystyle=\mathrm{d}\mathrm{cay}^{-1}_{h_{i}\,B({\mu}_{n,c_{i}}))}\,\mu_{n,c_{i}}$
(17a) $\displaystyle\mu_{n,r_{i}}$
$\displaystyle=\mathrm{d}\mathrm{cay}^{-1}_{-h_{i}\,B({\mu}_{n,c_{i}}))}\,\mu_{n,c_{i}}$
(17b)

where one solves implicitly the first equation to find $\mu_{n,c_{i}}$ from
the known $\mu_{n,r_{i-1}}$ and uses the solution in the second equation to
update to $\mu_{n,r_{i}}$. The $\mu_{n,r_{i}}$ and $\mu_{n,c_{i}}$ thus form a
staggered grid that can be used in conjunction to advance in time to the full
steps. The existence of solution of the implicit step is proven for
sufficiently small $h$ in (viviani2019minimal, 19, Lemma 3). See appendix A
for the right invariant version of the above algorithm. We summarize this
section with the following theorem

###### Theorem 3.1

Consider the reduced isospectral Hamiltonian system in Eq. 4 and the Butcher
tableau for a symplectic diagonally implicit Runge-Kutta method in Eq. 11.
Then Eq. 17 provides an isospectral integrator for the system.

The theorem follows directly from Eq. 15 and this provides an alternative
proof of Theorem 4 in the recent work by Modin and Vivani modin2019lie (13)
for the case of a symplectic diagonally implicit Runge-Kutta method.
Previously, Vivani viviani2019minimal (19) already obtained an analogous
result for the case of the implicit midpoint rule. Our extension to SDIRKs
partially verifies a conjecture that can be found in this work
(viviani2019minimal, 19, Remark 4). Furthermore, we showed that
(hairer2006geometric, 8, Theorem VI.4.4), which states that an $s$-stage SDIRK
can be decomposed in $s$ implicit midpoint rule steps with suitable step size,
also holds in the isospectral reduced case.

###### Remark 1

In the work by Calvo, Iserles, and Zanna calvo1997numerical (5), one has that
$\mu_{n,1/2}=(\mu_{n}+\mu_{n+1})/2$, i.e, the midpoint used in the update
equation is the average of two consecutive full algebra points. A similar
interpretation can be given for the isospectral integrators given by Eq. 17.
Defining the variables

$\displaystyle g_{n,c_{i}}$
$\displaystyle=g_{n,r_{i}}+\frac{h_{i}}{2}k_{n,i}^{g}=\frac{g_{n,r_{i}}+g_{n,r_{i+1}}}{2}$
$\displaystyle p_{n,c_{i}}$
$\displaystyle=g_{n,r_{i}}+\frac{h_{i}}{2}k_{n,i}^{p}=\frac{p_{n,r_{i}}+p_{n,r_{i+1}}}{2}$

we have that the intermediate half points
$\mu_{n,c_{i}}=g_{n,c_{i}}^{\dagger}p_{n,c_{i}}$ result from the averaging of
the intermediate points $(g_{n,r_{i}},p_{n,r_{i}})$ and
$(g_{n,r_{i+1}},p_{n,r_{i+1}})$ at the cotangent bundle level.

##  4 Variational symplectic diagonally implicit Runge-Kutta methods

From the general theory of reduction one expects that the isospectral
symplectic Runge-Kutta methods of Modin and Vivani that were derived through
Hamiltonian reduction and the momentum map can also be obtained through a
reduced discrete variational principle. In the following theorem we show that
this is indeed the case for isospectral SDIRKs. The Cayley transform, which
serendipitously appeared in the previous section, will be an integral part for
this. Our result also has a surprising connection to an earlier variational
Lie group integrator in the work by Gawlik et al. gawlik2011geometric (7), see
Remark 3.

###### Theorem 4.1

Let $G$ be a matrix Lie Group with Lie algebra $\mathfrak{g}$ and
$L:TG\to\mathbb{R}$ be a left invariant Lagrangian with left trivialization
$\ell:\mathfrak{g}\to\mathbb{R}$ given by

$\displaystyle\ell(\xi)=\frac{1}{2}\left\langle
B^{-1}(\xi^{\dagger}),\xi\right\rangle$

where $B:\mathfrak{g}^{*}\to\mathfrak{g}$ is smooth and invertible. Let
$\\{g_{n,r_{1}},g_{n,r_{2}},\ldots,g_{n,r_{s}}\\}_{n=0,\ldots,N-1}\subset G$
for a positive integer $s$ with $g_{n}=g_{n,r_{0}}=g_{n-1,r_{s}}$ and let
$\\{\xi_{n,c_{1}},\ldots,\xi_{n,c_{s}}\\}_{n=0,\ldots,N-1}\subset\mathfrak{g}$
be the associated sequence in the Lie algebra such that

$\displaystyle
g_{n,r_{i}}=g_{n,r_{i-1}}\tau\big{(}h_{i}\,\xi_{n,c_{i}}\big{)}$ (19)

with $\\{b_{i},c_{i}\\}$ being the coefficients of the Butcher tableau of an
$s$-stage SDIRK scheme. If the sequence $\\{g_{n,r_{i}}\\}$ is stationary for
the action $s_{d}:G^{N}\to\mathbb{R}$

$\displaystyle
s_{d}(\\{g_{n,r_{i}}\\})=\sum_{n=0}^{N-1}\sum_{i=1}^{s}h_{i}\,l(\xi_{n,c_{i}})$

under discrete variations of $\\{g_{n,r_{i}}\\}$ with fixed endpoints, then
the sequece $\\{\xi_{n,c_{i}}\\}$ and associated
$\mu_{n,c_{i}}=B^{-1}(\xi_{n,c_{i}}^{\dagger})\in\mathfrak{g}^{*}$ have to
satisfy the discrete Euler-Poincaré equations

$\displaystyle\begin{aligned}
\left(\mathrm{d}\tau^{-1}_{h_{i+1}\xi_{n,c_{i+1}}}\right)^{*}\mu_{n,c_{i+1}}&=\left(\mathrm{d}\tau^{-1}_{-h_{i}\xi_{n,c_{i}}}\right)^{*}\mu_{n,c_{i}}\qquad
i=1,\ldots,s-1\\\\[3.0pt]
\left(\mathrm{d}\tau^{-1}_{h_{1}\,\xi_{n+1,c_{1}}}\right)^{*}\mu_{n+1,c_{1}}&=\left(\mathrm{d}\tau^{-1}_{-h_{s}\,\xi_{n,c_{s}}}\right)^{*}\mu_{n,c_{s}}\qquad
i=s.\end{aligned}$ (20)

Furthermore, in the case of a quadratic Lie group and given a initial
condition $\mu_{0}$, the sequence $\\{\mu_{n}\\}$ defined by

$\displaystyle\mu_{n+1}=\left(\mathrm{Ad}_{\prod_{i=1}^{s}\tau(h_{j}\,\xi_{n,c_{j}})}\right)^{*}\mu_{n}=\mathrm{Ad}^{*}_{\left(\prod_{i=1}^{s}\tau(h_{j}\,\xi_{n,c_{j}})\right)^{-1}}\,\mu_{n}.$
(21)

agrees with the isospectral SDIRK method of the previous section.

Computationally, Eq. 20 allows to implicitly update the intermediate reduced
momenta $\mu_{n,c_{i}}$ and with Eq. 21 the next full time step can be
obtained. Eq. 20 is of Euler-Poincaré type by considering
$\mathrm{d}\tau_{-h_{i}\xi_{n,c_{i}}}^{-1}$ as the finite time analogue of the
infinitesimal coadjoint action $\mathrm{ad}_{\xi_{n,c_{i}}}^{*}$ in the
continuous Euler-Poincaré equations.

###### Proof

The stationarity of the discrete action under variations of the group sequence
$\\{g_{n,r_{i}}\\}$ with fixed endpoints means

$\displaystyle
0=\frac{d}{d\epsilon}\Big{|}_{\epsilon=0}s_{d}\left(\\{g_{n,r_{i}}^{\epsilon}\\}\right)=\sum_{n=0}^{N-1}\sum_{i=1}^{s}\frac{d}{d\epsilon}h_{i}\,l(\xi^{\epsilon}_{n,c_{i}})\,\Big{|}_{\epsilon=0}=\sum_{n=0}^{N-1}\sum_{i=1}^{s}h_{i}\,\left\langle\frac{\delta
l}{\delta\xi_{n,c_{i}}},\delta\xi_{n,c_{i}}\right\rangle$

with $g^{\epsilon}_{0}=g_{0}$, $g^{\epsilon}_{N}=g_{N}$ and
$\\{g^{0}_{n,r_{i}}\\}=\\{g_{n,r_{i}}\\}$ for all $n,i,\epsilon$. Through Eq.
19, a discrete variation in the group $\\{\delta
g_{n,r_{i}}\\}=\\{\frac{d}{d\epsilon}g^{\epsilon}_{n,r_{i}}\big{|}_{\epsilon=0}\\}$
induces a variation in the corresponding algebra sequence. Writing
$\eta=g^{-1}\delta g$ and using the properties of retraction maps, the induced
algebra variations are given by

$\displaystyle\delta\xi_{n,c_{1}}=\frac{1}{h_{1}}\left(-\mathrm{d}\tau^{-1}_{h_{1}\xi_{n,c_{1}}}(\eta_{n,r_{0}})+\mathrm{d}\tau^{-1}_{-h_{1}\xi_{n,c_{1}}}(\eta_{n,r_{1}})\right).$

Stationarity of the discrete action thus amounts to

$\displaystyle 0$
$\displaystyle=\sum_{n=0}^{N-1}\sum_{i=1}^{s}\,\left\langle\frac{\delta
l}{\delta\xi_{n,c_{i}}},\left(-\mathrm{d}\tau^{-1}_{h_{i}\xi_{n,c_{i}}}(\eta_{n,r_{i-1}})+\mathrm{d}\tau^{-1}_{-h_{i}\xi_{n,c_{i}}}(\eta_{n,r_{i}})\right)\right\rangle.$

Since $g_{0}^{\epsilon}=g_{0}$ and $g_{N}^{\epsilon}=g_{N}$ for every
$\epsilon$, we have $\delta g_{0}=\delta g_{N}=0$ and $\eta_{0}=\eta_{N}=0$.
Eq. 20 then follows by duality and using a standard summation by parts
argument. From $\mu_{n,c_{i}}=B^{-1}(\xi_{n,c_{i}}^{\dagger})$ it follows that
Eq. 21 and Eq. 15 are equivalent.

Figure 1: Comparison of isospectral midpoint rule (left) and the implicit
midpoint rule from Gawlik et al. gawlik2011geometric (7) (right) for change in
energy (top) and change in the eigenvalues (bottom).

###### Remark 2

In the continous case, identifying $\mathfrak{g}$ and $\mathfrak{g}^{*}$, the
Euler-Poincaré equation for an isospectral flow is

$\displaystyle\dot{\mu}=[\xi,\mu]$

The solution $\mu(t)$ and the group curve satisfying
$g^{-1}(t)\dot{g}(t)=\xi(t)$ are related by (hairer2006geometric, 8, IV.3.2)

$\displaystyle g(t)\,\mu(t)g(t)^{-1}=g(0)\,\mu(0)g(0)^{-1}\quad\text{ for all
$t$}.$

This ensures that the $\mu(t)$ are related by a similarity transformation for
all $t$ and therefore all possess the same spectrum. Eq. 21 are equivalent to
the discrete equation

$\displaystyle g_{n+1}\,\mu_{n+1}g_{n+1}^{-1}=g_{n}\,\mu_{n}g_{n}^{-1}$

which provides a direct verification that our integrator preserves the
isospectral property.

###### Remark 3

The simplest integrator that results from the above theorem is the implicit
midpoint rule studied by Vivani viviani2019minimal (19). Interestingly, it
also appeared in disguise in the work by Gawlik et al. gawlik2011geometric (7)
on variational discretizations of (possibly infinite dimensional) systems on
Lie groups. Combining a forward half-step of the implicit midpoint rule from
the intermediate point $\tilde{\mu}_{k}$ with a backward half-step from
$\tilde{\mu}_{k+1}$ one obtains the following implicit update rule for the
half-steps

$\displaystyle\mathrm{d}\mathrm{cay}^{-1}_{hB(\tilde{\mu}_{k+1})}\tilde{\mu}_{k+1}$
$\displaystyle=\mathrm{d}\mathrm{cay}^{-1}_{-hB(\tilde{\mu}_{k})}\tilde{\mu}_{k}$
(22)

or, equivalently,

$\displaystyle\left(\mathrm{Id}-\frac{hB(\tilde{\mu}_{k+1})}{2}\right)\tilde{\mu}_{k+1}\left(\mathrm{Id}+\frac{hB(\tilde{\mu}_{k+1})}{2}\right)$
$\displaystyle=\left(\mathrm{Id}+\frac{hB(\tilde{\mu}_{k})}{2}\right)\tilde{\mu}_{k}\left(\mathrm{Id}-\frac{hB(\tilde{\mu}_{k})}{2}\right).$

This are Eq. 4.9 and Eq. 4.12 in the work by Gawlik et al. gawlik2011geometric
(7). However, since only the half time steps are considered in this work their
integrator is not isospectral and also yields significantly stronger
oscillations in the energy, see Fig. 1.

##  5 Numerical Experiments

In this section we present numerical results for the isospectral SDIRKs from
Theorem 4.1 for the rigid body, the Toda lattice, and the matrix model
discretization of the barotropic vorticity equation on $S^{2}$. We begin with
some details on the computations. Our implementation is available online at
[will be made available upon publication of the manuscript].

###  5.1 Implementation of isospectral SDIRKs for quadratic Lie groups

For the case of isospectral systems defined on quadratic Lie algebras, Eq. 20
and Eq. 21 are equivalent to Eq. 17. Given an SDIRK and using the same
notation as in Sec. 2.3, Theorem 4.1 amounts to

$\displaystyle\begin{aligned}
\mu_{n,r_{i-1}}&=\left(\text{Id}-\frac{h_{i}}{2}B\left(\mu_{n,c_{i}}\right)\right)\mu_{n,c_{i}}\left(\text{Id}+\frac{h_{i}}{2}B\left(\mu_{n,c_{i}}\right)\right)\\\\[3.0pt]
\mu_{n,r_{i}}&=\left(\text{Id}+\frac{h_{i}}{2}B\left(\mu_{n,c_{i}}\right)\right)\mu_{n,c_{i}}\left(\text{Id}-\frac{h_{i}}{2}B\left(\mu_{n,c_{i}}\right)\right)\end{aligned}$
(23)

which allows to update the $\mu_{n,r_{i}}$ and $\mu_{n,c_{i}}$ by leap
frogging between them. We summarize the overall computations in Algo. 1.

In the experiments, we will consider the $2^{\textrm{nd}}$-order implicit
midpoint rule and the isospectral integrators obtained from a $2$-stage
$2^{\textrm{nd}}$-order SDIRK kalogiratou2012diagonally (10) and a $7$-stage
$4^{\textrm{th}}$-order SDIRK kalogiratou2012diagonally (10).

Figure 2: Change in energy and change in eigenvalues for the rigid body for
implicit midpoint rule (top), a $2$-stage $2^{\textrm{nd}}$-order scheme
(middle) and a $7$-stage $4^{\textrm{nd}}$-order scheme (bottom).

###  5.2 Rigid body

The ideal rigid body is the classical example of an isospectral system. Its
dynamics are described by Euler’s equations, which, in matrix form, are
(Marsden1999, 11, Ch. 15)

$\displaystyle\dot{W}=[\nabla H(W)^{\dagger},W]$ (24)

where $W\in\mathfrak{so}^{*}(3)$ and the Hamiltonian is
$H(W)=\frac{1}{2}\left\langle\mathcal{I}^{-1}W,W\right\rangle$ with
$\mathcal{I}$ a symmetric $3\times 3$ matrix.

Fig. 2 shows results for the three isospectral integrators we consider for a
random skew-symmetric matrix with entries in $[-1,1]$ as initial condition and
with a time step of $h=0.01$. The superiority of the $4^{\textrm{th}}$-order
integrator is evident in the conservation of the Hamiltonian, which is almost
on the order of machine precision. The conservation of the angular momentum
(not shown) is of the same order as those for the eigenvalues.

###  5.3 Toda lattice

The Toda lattice describes the pairwise interactions of particles on a line
exerting exponential forces on each other Toda1967 (18). When the line is
periodic, the Lax pair formulation of the Toda lattice (hairer2006geometric,
8, Ch. IV.3.2) results in the isospectral flow $\dot{L}=\left[B(L),L\right]$
with Hamiltonian $H(L)=2\text{Tr}(L^{2})$ for

Figure 3: Change in energy (left) and change in eigenvalues (right) for the
isospectral implicit midpoint rule and the $4^{\textrm{th}}$-order SDIRK for
the Toda lattice.

$\displaystyle L=\begin{pmatrix}a_{1}&b_{1}&0&\ldots&b_{n}\\\
b_{1}&a_{2}&b_{2}&\ldots&0\\\ 0&b_{2}&a_{3}&\ldots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
b_{n}&0&0&\ldots&a_{n}\end{pmatrix},\quad
B(L)=\begin{pmatrix}0&b_{1}&0&\ldots&-b_{n}\\\ -b_{1}&0&b_{2}&\ldots&0\\\
0&-b_{2}&0&\ldots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
b_{n}&0&0&\ldots&0\end{pmatrix}$

While the Toda lattice is isospectral, it is not Lie-Poisson. It can, however,
be generalized to an isospectral flow in $\mathfrak{gl}(n)$ by introducing

$\displaystyle B(W)=\begin{pmatrix}0&W_{12}&0&\ldots&-W_{1n}\\\
-W_{21}&0&W_{23}&\ldots&0\\\ 0&0-W_{32}&0&\ldots&0\\\
\vdots&\vdots&\vdots&\ldots&\vdots\\\ W_{n1}&0&0&\ldots&0\\\ \end{pmatrix}$

which corresponds to the extended Hamiltonian
$\tilde{H}(W)=-\text{Tr}\left(W^{\dagger}B(W)\right)+H(W)$. The dynamics of
the extended Lie-Poisson system are
$\dot{W}=[\mathrm{d}\tilde{H}(W)^{\dagger},W]$ which has the form of Eq. 4.

In Fig. 3 we show experimental results for the periodic Toda lattice for the
implicit midpoint rule and the $7$-stage $4^{\textrm{th}}$-order Runge-Kutta
scheme. In both cases we obtain excellent conservation of the eigenvalues and
the extended Hamiltonian. The initial condition $W_{0}$ is the same as used in
modin2019lie (13) where $a_{i}=b_{i}=(-1)^{i},i=1,\ldots,4$ and we used
$h=0.1$ for the time step.

###  5.4 Matrix model of barotropic vorticity equation on $S^{2}$

The barotropic vorticity equation on the sphere is given by

$\displaystyle\dot{\zeta}=\\{\zeta,\Delta^{-1}\zeta\\},$ (25)

with $\zeta$ being vorticity and $\\{\,,\\}$ the canonical Poisson bracket on
$S^{2}$. Eq. 25 is a simple model for weather and climate dynamics and plays
an important role in various applications.

Figure 4: Change in energy (top), eigenvalues (middle) and enstrophy (bottom)
for the isospectral implicit midpoint rule and the $4^{\textrm{th}}$-order
SDIRK for the matrix model of the barotropic vorticity equation. We used a
truncation at $L=16$ ($N=33$).

A structure preserving discretization for it is the matrix model by Zeitlin
Zeitlin2004 (21) that goes back to early work for the torus Zeitlin1991 (20)
and results by Hoppe on matrix harmonics for $S^{2}$ Hoppe1998 (9). Recently,
Modin and Vivani modin2019lie (13) used it to study turbulence on the sphere.
The governing equation for the matrix model is

$\displaystyle\dot{W}=N^{3/2}\,[\hat{\Delta}W,W]$ (26)

where $\hat{\Delta}$ is a discrete Laplacian Hoppe1998 (9) and
$W\in\mathfrak{sl}(N)$ with $N=L+1$ is the representation of vorticity. It is
related to the spatial vorticity field by

$\displaystyle W=\sum_{l=1}^{L}\sum_{m=-l}^{l}\zeta_{lm}\,T_{lm}^{(L)}$

where the $\zeta_{lm}$ are the spectral coefficients of vorticity with respect
to spherical harmonics and the $T_{lm}$ are matrix harmonics Hoppe1998 (9).

Eq. 26 is an isospectral flow on $\mathfrak{sl}(N)$ with conserved quantities

$\displaystyle C_{n}=\mathrm{tr}(W^{n}),\quad n=1,\ldots,N,$

which are the discrete analoges to the integrated powers of vorticity that are
conserved for the continuous barotropic vorticity equation Khesin1998 (1).

In Fig. 4 we report our experimental results for the $2^{\textrm{nd}}$-order
implicit midpoint rules and the $4^{\textrm{th}}$-order SDIRK. The initial
condition $W_{0}$ was a randomly generated matrix in $\mathfrak{sl}(N)$ and we
use $h=0.0025$ (note that the matrix model results in a time rescaling, cf.
(modin2020casimir, 14, sec. 2.4)). Both isospectral time integration schemes
preserve the eigenvalues of $W$ up to machine precision, as expected, and this
carries over to the Casimirs, as exemplified for the discrete enstrophy
$C_{2}$.

Algorithm 1 Isospectral SDIRK

1:Input: $s$, $(b_{1},\ldots,b_{s}),N_{\text{steps}},h,\mu_{0}$,

2: $\mu_{0,r_{0}}=\mu_{0}$

3:for $n\leftarrow 1,\,N_{\text{steps}}$ do $\triangleright$ Time-stepping

4: for $i\leftarrow 1,\,s$ do

5: $h_{i}\leftarrow b_{i}\,h$

6: Solve
$\mu_{n,r_{i-1}}=\left(\text{Id}-\frac{h_{i}}{2}B\left(\mu_{n,c_{i}}\right)\right)\mu_{n,c_{i}}\left(\text{Id}+\frac{h_{i}}{2}B\left(\mu_{n,c_{i}}\right)\right)$
for $\mu_{n,c_{i}}$ $\triangleright$ nonlinear solve

7:
$\mu_{n,r_{i}}\leftarrow\left(\text{Id}+\frac{h_{i}}{2}B\left(\mu_{n,c_{i}}\right)\right)\mu_{n,c_{i}}\left(\text{Id}-\frac{h_{i}}{2}B\left(\mu_{n,c_{i}}\right)\right)$
$\triangleright$ explicit update

8: end for

9:end for

##  6 Conclusion

In this work, we showed that isospectral symplectic diagonally implicit Runge-
Kutta methods (SDIRKs) can be derived from the discrete Hamilton’s principle
of BouRabee2007 (4, 7). The resulting numerical schemes are the simplest and
most easily implemented isospectral higher order Runge-Kutta methods. We
demonstrate this with our numerical results for the rigid body, the Toda
lattice and the discrete 2D Euler fluid. Our results for a
$4^{\textrm{th}}$-order SDIRK yield a time integration that preserve energy
and isospectrality almost up to machine precision. In the available precision,
it is hence a faithful discretization of the system.

Some interesting directions for future work are as follow. Theorem 4.1 only
applies to SDIRKs and a generalization to arbitrary symplectic Runge-Kutta
methods is hence open. It is currently also not clear to us why the derivation
of Gawlik et al. gawlik2011geometric (7) yields the midpoints of the
isospectral implicit midpoint rules by Vivani viviani2019minimal (19). We
believe that this might shed more light onto the overall structure of discrete
variational principles for isospectral systems. An interesting objective for
future work is also to extend our results to systems whose configuration space
is a semi-direct product. In some cases, such as the heavy top, then also an
infinite number of conserved quantities exist. Bobenko and Suris Bobenko1999
(2) studied this already for the special case of the Lagrange top although
integrability is no longer given in the general case.

###### Acknowledgements.

We would like to acknowledge helpful discussions with Mathieu Desbrun.

##  Appendix A Right invariance

The canonical hamiltonian system in Eq. 5 resulting from the reconstruction of
the isospectral system in Eq. 4 is left invariant. A right invariant version
of the canonical hamiltonian system and the algorithm in Eq. 17 can be
achieved as follows.

Consider the left action of the Lie group G on $T^{*}G$ given by
$h\cdot(g,p)=(hg,(h^{\dagger})^{-1})$. It induces the momentum map
$\mathcal{J}(g,p)=pg^{\dagger}$. Defining the right invariant Hamiltonian on
$T^{*}G$ as $\tilde{H}(g,p)=H\circ\mathcal{J}(g,p)$, where $H$ is the
Hamiltonian for the Lie-Poisson system in Eq. 4, we have the canonical right
invariant Hamiltonian system on $T^{*}G$

$\displaystyle\dot{g}$ $\displaystyle=\nabla
H(g^{\dagger}p)g=B(g^{\dagger}p)^{\dagger}g$ $\displaystyle\dot{p}$
$\displaystyle=-\nabla H(g^{\dagger}p)^{\dagger}p=-B(g^{\dagger}p)p$

where again $B(\cdot)=\nabla H(\cdot)^{\dagger}$. Applying a symplectic
diagonally implicit Runge-Kutta method to the above equation and proceeding as
in section 3, we have a right invariant version of the algorithm in Eq. 17,

$\displaystyle\mu_{n,r_{i-1}}$
$\displaystyle=\mathrm{d}\mathrm{cay}^{-1}_{-h_{i}\,B({\mu}_{n,c_{i}}))}\,\mu_{n,c_{i}}$
$\displaystyle\mu_{n,r_{i}}$
$\displaystyle=\mathrm{d}\mathrm{cay}^{-1}_{h_{i}\,B({\mu}_{n,c_{i}}))}\,\mu_{n,c_{i}}.$

Notice that the correct time direction is only obtain when one distinguishes
the left and right invariant case. Theorem 4.1 also has a right invariant
version where the update equation in Eq. 19 becomes
$g_{n,r_{i}}=\tau(h_{i}\xi_{n,c_{i}})g_{n,r_{i-1}}$.

## References

  * (1) V. I. Arnold and B. A. Khesin  “Topological Methods in Hydrodynamics”, Applied Mathematical Sciences  New York: Springer, 1998
  * (2) A. I. Bobenko and Yu. B. Suris  “Discrete Time Lagrangian Mechanics on Lie Groups,¶with an Application to the Lagrange Top”  In _Communications in Mathematical Physics_ 204.1, 1999, pp. 147–188  DOI: 10.1007/s002200050642
  * (3) N. Bou-Rabee and J. E. Marsden  “Hamilton–Pontryagin Integrators on Lie Groups Part I: Introduction and Structure-Preserving Properties”  In _Found. Comput. Math._ 9.2 Berlin, Heidelberg: Springer-Verlag, 2009, pp. 197–219
  * (4) N. M. Bou-Rabee  “Hamilton-Pontryagin Integrators on Lie Groups”, 2007  URL: https://resolver.caltech.edu/CaltechETD:etd-06052007-153115
  * (5) M. Calvo, A. Iserles and A. Zanna  “Numerical solution of isospectral flows”  In _Mathematics of computation_ 66.220, 1997, pp. 1461–1486
  * (6) F. Diele, L. Lopez and T. Politi  “One Step Semi-Explicit Methods Based on the Cayley Transform for Solving Isospectral Flows”  In _J. Comput. Appl. Math._ 89.2 NLD: Elsevier Science Publishers B. V., 1998, pp. 219–223  DOI: 10.1016/S0377-0427(97)00236-7
  * (7) E. S. Gawlik et al.  “Geometric, variational discretization of continuum theories”  In _Physica D: Nonlinear Phenomena_ 240.21 Elsevier, 2011, pp. 1724–1760
  * (8) E. Hairer, C. Lubich and G. Wanner  “Geometric Numerical Integration”, Springer Series in Computational Mathematics  Springer-Verlag, 2006
  * (9) J. Hoppe and S.-T. Yau  “Some Properties of Matrix Harmonics on S2”  In _Communications in Mathematical Physics_ 195.1 Springer-Verlag, 1998, pp. 67–77  DOI: 10.1007/s002200050379
  * (10) Z Kalogiratou, Th Monovasilis and TE Simos  “Diagonally implicit symplectic Runge-Kutta methods with special properties”  In _AIP Conference Proceedings_ 1479.1, 2012, pp. 1387–1390  American Institute of Physics
  * (11) J. E. Marsden and T. S. Ratiu  “Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems”, Texts in Applied Mathematics  New York: Springer-Verlag, 1999
  * (12) J. E. Marsden and M. West  “Discrete mechanics and variational integrators”  In _Acta Numerica_ 10 Cambridge University Press, 2001, pp. 357–514
  * (13) K. Modin and M. Viviani  “Lie–Poisson methods for isospectral flows”  In _Foundations of Computational Mathematics_ Springer, 2019, pp. 1–33
  * (14) Klas Modin and Milo Viviani  “A Casimir preserving scheme for long-time simulation of spherical ideal hydrodynamics”  In _Journal of Fluid Mechanics_ 884 Cambridge University Press, 2020
  * (15) J. B. Moore, R. E. Mahony and U. Helmke  “Numerical Gradient Algorithms for Eigenvalue and Singular Value Calculations”  In _SIAM Journal on Matrix Analysis and Applications_ 15.3, 1994, pp. 881–902  DOI: 10.1137/S0036141092229732
  * (16) J. Moser and A. P. Veselov  “Discrete versions of some classical integrable systems and factorization of matrix polynomials”  In _Communications in Mathematical Physics_ 139.2 Springer, 1991, pp. 217 –243  DOI: cmp/1104203302
  * (17) B. Owren  “Lie group integrators”  In _Discrete Mechanics, Geometric Integration and Lie–Butcher Series_ , Springer Proceedings in Mathematics & Statistics 267  Springer, 2018  DOI: https://doi.org/10.1007/978-3-030-01397-4˙2
  * (18) M. Toda  “Vibration of a Chain with Nonlinear Interaction”  In _Journal of the Physical Society of Japan_ 22.2, 1967, pp. 431–436  DOI: 10.1143/JPSJ.22.431
  * (19) M. Viviani  “A minimal-variable symplectic method for isospectral flows”  In _BIT Numerical Mathematics_ Springer, 2019, pp. 1–18
  * (20) V. Yu. Zeitlin  “Finite-mode analogs of 2D ideal hydrodynamics: Coadjoint orbits and local canonical structure”  In _Physica D: Nonlinear Phenomena_ 49.3 North-Holland, 1991, pp. 353–362  DOI: 10.1016/0167-2789(91)90152-Y
  * (21) V. Yu. Zeitlin  “Self-Consistent Finite-Mode Approximations for the Hydrodynamics of an Incompressible Fluid on Nonrotating and Rotating Spheres”  In _Physical Review Letters_ 93.26 American Physical Society, 2004, pp. 264501  DOI: 10.1103/PhysRevLett.93.264501


# A simple but effective bushfire model:
analysis and real-time simulations

Serena Dipierro  S. D., Department of Mathematics and Statistics, University
of Western Australia, 35 Stirling Highway, WA 6009 Crawley, Australia.
<EMAIL_ADDRESS>,  Enrico Valdinoci  E. V., Department of
Mathematics and Statistics, University of Western Australia, 35 Stirling
Highway, WA 6009 Crawley, Australia<EMAIL_ADDRESS>,  Glen
Wheeler  G. W., School of Mathematics and Applied Statistics, University of
Wollongong, Northfields Avenue, NSW 2500 Wollongong, Australia.
<EMAIL_ADDRESS>and  Valentina-Mira Wheeler  V.-M. W., School of Mathematics
and Applied Statistics, University of Wollongong, Northfields Avenue, NSW 2500
Wollongong, Australia<EMAIL_ADDRESS>
###### Abstract.

We introduce a simple mathematical model for bushfires accounting for
temperature diffusion in the presence of a combustion term which is activated
above a given ignition state. The model also takes into consideration the
effect of the environmental wind and of the pyrogenic flow.

The simplicity of the model is highlighted from the fact that it is described
by a single scalar equation, containing only four terms, making it very handy
for rapid and effective numerical simulations which run in real-time.

In spite of its simplicity, the model is in agreement with data collected from
bushfire experiments in the lab, as well as with spreading of bushfires that
have been observed in the real world.

Moreover, the equation describing the temperature evolution can be easily
linked to a geometric evolution problem describing the level sets of the
ignition state.

Supported by Australian Research Council FL190100081 and DE190100379.

##  1\. The bushfire model

The study of bushfire propagation requires, in general, the understanding of
very complex and intertwined phenomena, involving fluid mechanics, physics,
chemistry, meteorology, and thermodynamics. In this spirit, a thorough
bushfire analysis should combine weather and climate models with fuel and
terrain models, describing accurately flaming combustion phenomena and eddy
formation in different spatial and temporal scales (see e.g. [BAK]). The
“ultimate” universal bushfire model is not available yet, and it is customary
to approach the problem under different perspectives, so that often the
different models are classified into the categories of “physical” (or
“theoretical”), “semi-physical” (or “semi-empirical”), and “empirical”, see
e.g. [PERRY, PASTOR]; sometimes “mathematical” models have been considered as
a separate category closely related to “simulation”, see e.g. [SULL1, SULL2,
SULL3].

Of course, on the one hand, it is desirable to describe bushfires by a model
which is as complete and accurate as possible, but, on the other hand, it
would also be important to have “simple” (though not simplistic) models
allowing for rigorous mathematical treatments, a qualitative understanding of
the main processes involved, as well as fast simulations making the model
reliable but also concretely usable in real-time.

With these considerations in mind, the model that we propose consists of a
single scalar equation, in which the effects coming from fluid dynamics and
chemical combustion are encoded by state parameters and interaction kernels.
More specifically, the equation on which we focus our attention has the form

$\begin{split}\partial_{t}u&=c\Delta
u+\int_{\mathds{R}^{n}}\big{(}u(y,t)-\Theta(y,t)\big{)}_{+}\,K(x,y)\,dy+\left(\left(\omega+\frac{\beta(u)\nabla
u}{|\nabla u|^{\gamma}}\right)\cdot\nabla u\right)_{-}\,,\end{split}$ (1.1)

where we used the notation for the “positive and negative parts” of a function

$f_{+}(x):=\max\\{0,f(x)\\}$ and $f_{-}(x):=\max\\{0,-f(x)\\}$. (1.2)

In this model, $u=u(x,t)$ represents the temperature of the environment at the
space point $x\in\mathds{R}^{n}$, at time $t$ (in many applications, one can
focus on the case $n=2$, and deal with a given region with suitable, e.g.
Dirichlet, conditions).

As customary, one assumes that temperature spreads through a heat equation, in
which $c=c(x,t)$ denotes the environmental diffusion coefficient (the case of
constant $c$ modeling a homogeneous environment).

The specific model for bushfires accounts for three additional phenomena,
namely combustion, wind, and convection.

The combustion is described by an effective ignition temperature
$\Theta=\Theta(x,t)$, above which the fuel burns. The case of a constant
$\Theta(x,t)=\Theta_{0}$ represents the ideal situation of fuel homogeneously
spread across the environment, with infinite abundance of fuel at every point.
In this case, the fuel burns at temperature $\Theta_{0}$ and affects nearby
points through the interaction potential $K$ (notice that when the temperature
$u$ of the fuel is below the effective ignition temperature $\Theta$, the
combustion effect disappears, in view of the positive part of the function
considered in the integral).

Situations in which the availability of fuel varies in space and time are of
course modeled by nonconstant functions $\Theta$. We stress that the model can
also include a memory effect, in which the burned fuel turns into ash and is
not any longer subject to combustion. This can be done, for instance, by
considering the case in which

$\Theta(x,t):=\begin{cases}\displaystyle\tan\left(\frac{\pi}{2F}\int_{0}^{t}\big{(}u(x,\tau)-\bar{\Theta}(x)\big{)}_{+}\,d\tau\right)&{\mbox{if
}}\displaystyle\int_{0}^{t}\big{(}u(x,\tau)-\bar{\Theta}(x)\big{)}_{+}\,d\tau<F,\\\
+\infty&{\mbox{otherwise.}}\end{cases}$ (1.3)

In this situation, the constitutional ignition temperature at a given point
$x$ is $\bar{\Theta}(x)$ (the case of a constant ignition temperature
$\bar{\Theta}(x)=\Theta_{0}$ accounting for a uniformly distributed fuel).
Above this ignition temperature, the fuel burns, and the amount of fuel burned
is described by the quantity

$B(x,t):=\int_{0}^{t}\big{(}u(x,\tau)-\bar{\Theta}(x)\big{)}_{+}\,d\tau.$
(1.4)

In this situation, $F=F(x,t)$ accounts for the fuel available (for instance,
the case $F=F(x)$ corresponds to the situation of an initial amount of fuel
$F(x)$ available at the point $x$, and the case of $F$ also depending on time
allows one to consider situations in which the fuel is recharged over time).

We highlight the difference between the constitutional ignition temperature
$\bar{\Theta}$ and the effective ignition temperature $\Theta$. Namely, while
$\bar{\Theta}$ can be thought as a specific property of the combustible (which
is modeled to ignite above this temperature, in the ideal situation of
infinite amount of fuel available), the quantity $\Theta$ takes into account
the past state of the fuel and the propellant $B$ already burned. Indeed,
comparing (1.3) and (1.4), the effective ignition temperature in this case can
be written in the form

$\Theta=\tan\left(\frac{\pi B}{2F}\right),$

hence when the quantity of burned combustible $B$ reaches the threshold $F$,
the fuel gets extinguished, and the corresponding effective ignition
temperature diverges to $+\infty$, making the interaction ignition term in
(1.1) vanish, due to the presence of the positive part in the equation.

See also [MR2542721] for models taking into account the environmental
temperature with a fuel ignition term.

The effect of the environmental wind is modeled by a vector
$\omega=\omega(x,t)$ inducing a transport term in equation (1.1).

The convection is encoded by an additional wind term of the form
$\frac{\beta(u)\nabla u}{|\nabla u|^{\gamma}}$. The function $\beta$ can be
thought as supported above the ignition temperature, so that this term is
active only in the burning region, directed along the gradient temperature.
The exponent $\gamma$ can be taken in the interval $[0,2]$ to further modulate
the intensity of this term with respect to the temperature gradient.

The negative part in the term regarding the environmental wind and the
convective term accounts for the fact that wind can spread the fire faster,
but cannot stop it, hence these terms provide a contribution only when their
direction is oriented along the gradient temperature, which is in turn a proxy
for the direction of the front. In this spirit, the role of the negative part
in the last term of (1.1) is to spread the fire in the direction of the wind:
as the simulations confirm, without the negative part, the burning region will
tend to mainly translate in the direction of the wind.

Terms related to a fire-induced wind (also known as “fire wind” or “pyrogenic
flow”) have been introduced in the literature to account for the significant
buoyant upflow over the fire region, created by the strong temperature
gradient, producing local low pressure, which acts as a sink with horizontal
pressure gradient and draws in the surrounding air (see [SMITH]).

The possibility of modeling the wind advection by a transport term, in various
forms, has also been previously considered in the literature, see e.g.
[MR4439510, equation (2.2)]. In terms of concrete applications, the
combination of pyrogenic and environmental wind can be a delicate issue:
indeed, on the one hand, in case of wind-driven fires it is expected the
effect of pyrogenic wind to be less prominent (see [BEER]), on the other hand
pyrogenic winds can be produced even by relatively small wildfires (as
observed in [LAREAU]). See [HILTON1, HILTON2, HILTON201812, MR4439510] and the
references therein for different models describing bushfires accounting for
wind and convection effects.

The model that we propose is also prone to modifications and generalizations
to capture additional details. For instance, the interaction kernel $K$ is
supposed here to be a nonnegative function (the case in which
$K(x,y)=\bar{K}(|x-y|)$ describing an ignition interaction depending only on
the distance between sites), but one can replace the term $K(x,y)\,dy$ by a
more general measure integral: for instance, integrating against a Delta
Function would return a source term of the type
$\big{(}u(x,t)-\Theta(x,t)\big{)}_{+}$.

In the same spirit, the terms linear in $(u-\Theta)_{+}$ can be replaced by
nonlinear modifications, to capture nonlinear scalings in fire propagation and
radiation effects, the interaction kernels may also depend on the temperature,
and the Laplacian can be substitute with other diffusive operators, e.g. to
model various types of anisotropy and inhomogeneity in the terrain, the
negative part in the wind term can be dropped, e.g. to model possible
containment effects of the wind.

The Laplacian can also be replaced, or complemented, by nonlinear or nonlocal
operators, with the aim of modeling anomalous types of diffusion, and also by
operators of mixed orders, for instance to describe the combination of two
types of diffusion (such as the spread of heat along the ground surface and
the boundary effect on the ground coming from the spread of heat in the air).
For the sake of simplicity, we do not present all these generalizations in
full detail here, but we plan to address them specifically in future projects.

As a sanity check, let us remark that, since both the positive and negative
parts of a function are, in our notation (1.2), nonnegative, the solutions of
(1.1) are always supersolutions of the heat equation (which is consistent,
from the physical point of view, with the fact that bushfires can only
increase the environmental temperature).

Albeit non-standard, the partial differential equation in (1.1) is of
parabolic type and is endowed with a solid mathematical theory, developed in
[PAPER2].

##  2\. Numerical validation and real-time simulations

The model presented in (1.1) is simple enough to allow for an agile and
effective numerical implementation, producing sufficiently precise results in
real-time. As concrete examples, we confronted our numerical outputs with
supervised experiments of fires performed in the lab as well as with real data
collected on the occasion of natural bushfires.

Numerical simulations were written in C++ with a standard finite-element
method, using the library FreeFEM++ [HECHT]. Our simulations used the weak
formulation of the PDE (1.1) for $u$ on a unit square with Dirichlet zero
boundary conditions. For the convolution kernel $K(x,y)$ we used an
approximation to the Dirac mass, which allows close non-local interactions and
facilitates fast computation. For parameters, we used $c=10^{-3}$,
$\Theta(y,t)=1$, and $\beta(u)=0$. Figure 1 uses wind with $\omega=(-1,0.4)$,
otherwise $\omega=0$. For the finite element method we used piecewise
continuous elements and divided the domain into a mesh with 10,000 individual
elements. Our simulations were run on a standard laptop without special
processing power.

We remark that we did not make special effort to change parameters in (1.1)
drastically to fit pictures, but instead tried to make only natural, broad
choices. In a real-world scenario, we expect that parameters are at best
approximated, so fine-tuning them would be counterproductive in our view.

Figure 1.  Observed isochrones for spot fire E14 investigated in [CRUZ],
compared to isochrones from our simulation taken at evenly spaced iterations
of the scheme. The red line is the fire front with temperature above ignition,
and temperature decreasing below ignition as the colours progress from blue to
green, then yellow.

Specifically, in Figure 1 we compare our numerical simulations with one of the
outdoor spot fire experiments conducted in Ridgewood by the CSIRO in NSW
Australia [CRUZ]. Named E14, the isochrones of the fire within 90 seconds from
its ignition are reproduced on the left, as recorded in [CRUZ]. For fire E14,
the authors in [CRUZ] note a change of wind direction by approximately 40
degrees during its burning. Our simulation has a constant wind direction.

On the right, we present the output of our numerical simulation using the
model in (1.1). The specifics of this numerical validation considered the fire
emerging from a point, with wind WNW and constant (24 degrees above the West
direction). There is a Dirichlet zero boundary condition imposed on the
boundary of the fire domain. There is a slight difference in shape from our
simulation compared to the experiment, which is likely due to the volatile
wind in the experiment.

Figure 2. Experiment CF-14 from [VIEGAS], with a 30 degree cone angle. Frames
taken from an infrared camera at thirty second intervals. Simulation frames
taken at evenly spaced iterations of the numerical simulation.

In Figures 2, 3, and 4 we compare our numerical simulations with experiments
of fires ran in labs. Namely, Figure 2 reproduces (on the top) the thermal
images of the experiment in [VIEGAS] of a V-shaped fire, to be confronted (on
the bottom) with our numerical simulations of the model in (1.1). This shows
that our model works well even in the challenging case of fire-line merging;
the movement of the vertex of the fire is quite close to the observed vertex,
cf. [W15].

Figure 3. Experimental vs. Simulation Frames. Experiment conducted at the
Pyrotron CSIRO Canberra February 2024. Experimental frames taken at one-minute
intervals. A magenta outline is traced around the approximate fire line.
Simulation frames taken at evenly spaced iterations of the numerical
simulation.

In Figures 3 and 4, we report about an experiment recently taken place at
Pyrotron CSIRO Canberra, tracing the fire lines emanating from an initial
burning square. The photos taken in the lab are confronted in Figures 3 with
our simulations.

This comparison is retaken also in Figure 4, taking the approximated burned
area as a measure of validation. The matching is almost perfect after up to 3
minutes, and a small deviation takes place at minute 4 and minute 5. This
small deviation has also to be expected, since our simulations are based on
ideal Dirichlet conditions, while the experimental fire may be affected by
slightly approximated boundary conditions allowing for a faster diffusion near
the boundary. Note also the coarseness of the fuel used for the experimental
fire.

Figure 4. Comparison of area within fire line for the experiment vs the
simulation displayed in Figure 3.

##  3\. The hybrid geometric problem

For constant ignition temperature $\Theta_{0}$, one can consider the level set
$\\{u=\Theta_{0}\\}$ (or, more precisely, the boundary of
$\\{u>\Theta_{0}\\}$) as the moving front of the fire. The velocity of this
front in the external normal direction, in light of (1.1), is

$-cH+\frac{c\partial^{2}_{\nu}u}{|\nabla u|}+\frac{1}{|\nabla
u|}\int_{\mathds{R}^{n}}\big{(}u(y,t)-\Theta_{0}\big{)}_{+}\,K(x,y)\,dy+\left(\left(\omega+\frac{\beta(u)\nabla
u}{|\nabla u|^{\gamma}}\right)\cdot\nabla u\right)_{-}\,,$ (3.1)

where $H$ is the mean curvature of $\partial\\{u>\Theta_{0}\\}$ and $\nu$ the
external unit normal.

The evolution equation in (3.1) has an interesting structure since it mixes a
geometric term of curvature type with a forcing term arising from the global
pattern of the solution of (1.1). This “hybrid” structure reflects the
physical property that the moving front of the bushfire is not modeled solely
by its geometry, since it is influenced by all the structural and
environmental parameters, such as temperature, wind, and convection.

Similar, but structurally different, equations have been obtained, or imposed
a-priori, in other bushfire models, see e.g. [MR4439510, equation (2.6)], in
which the outward normal velocity of a fire line in the plane was set to be
equal (up to normalizing constants) to a constant (inherited by free boundary
problems of Stefan type and accounting for oxygen gradient), minus the
curvature (which was proposed in [MARK] as a corrective term to match the
experimental data and acts as a smoothing term on perturbations of the fire
line, see e.g. [MR4439510, Figure 5]), plus a dynamical term requiring the
solution of an advection-diffusion for the oxygen concentration.

Interestingly, in our model, the normal velocity of the front is not imposed
a-priori, but follows from the constitutive equation (1.1). To see how
formally deduce (3.1) from (1.1) when $\Theta=\Theta_{0}$, let us consider the
level set $\\{u=\Theta_{0}\\}$ and suppose, up to a rigid motion, that it
passes through the origin of $\mathds{R}^{n}$ at time $t=0$, and that, in a
neighborhood of the origin, the burning region $\\{u>\Theta_{0}\\}$
corresponds to $\\{x_{n}>\gamma(x^{\prime},t)\\}$ for some
$\gamma:\mathds{R}^{n-1}\times\mathds{R}\to\mathds{R}$, with $\gamma(0,0)=0$
and

$\nabla_{x^{\prime}}\gamma(0,0)=0.$ (3.2)

Then, differentiating the relation
$u(x^{\prime},\gamma(x^{\prime},t),t)=\Theta_{0}$ with respect to the space
variables, we obtain that, for each $j$, $k\in\\{1,\dots,n-1\\}$,

$\displaystyle\partial_{x_{j}}u(x^{\prime},\gamma(x^{\prime},t),t)+\partial_{x_{n}}u(x^{\prime},\gamma(x^{\prime},t),t)\,\partial_{x_{j}}\gamma(x^{\prime},t)=0$
and
$\displaystyle\partial_{x_{j}x_{k}}^{2}u(x^{\prime},\gamma(x^{\prime},t),t)+\partial_{x_{j}x_{n}}^{2}u(x^{\prime},\gamma(x^{\prime},t),t)\,\partial_{x_{k}}\gamma(x^{\prime},t)+\partial_{x_{k}x_{n}}^{2}u(x^{\prime},\gamma(x^{\prime},t),t)\,\partial_{x_{j}}\gamma(x^{\prime},t)$
$\displaystyle\qquad+\partial_{x_{n}}^{2}u(x^{\prime},\gamma(x^{\prime},t),t)\,\partial_{x_{j}}\gamma(x^{\prime},t)\,\partial_{x_{k}}\gamma(x^{\prime},t)+\partial_{x_{n}}u(x^{\prime},\gamma(x^{\prime},t),t)\,\partial^{2}_{x_{j}x_{k}}\gamma(x^{\prime},t)=0.$

We now evaluate the last equation at $x^{\prime}=0$ and $t=0$, using (3.2),
and we conclude that, for each $j$, $k\in\\{1,\dots,n-1\\}$,

$\partial_{x_{j}x_{k}}^{2}u(0,0)+\partial_{x_{n}}u(0,0)\,\partial^{2}_{x_{j}x_{k}}\gamma(0,0)=0$

and therefore, choosing $j=k$ and summing up,

$\begin{split}0&=\sum_{k=1}^{n-1}\partial_{x_{k}}^{2}u(0,0)+\partial_{x_{n}}u(0,0)\,\Delta\gamma(0,0)\\\
&=\Delta u(0,0)-\partial_{x_{n}}^{2}u(0,0)+|\nabla u(0,0)|\,H(0,0)\end{split}$
(3.3)

Besides, differentiating the relation
$u(x^{\prime},\gamma(x^{\prime},t),t)=\Theta_{0}$ with respect to the time
variables, we obtain that

$\displaystyle\partial_{x_{n}}u(x^{\prime},\gamma(x^{\prime},t),t)\,\partial_{t}\gamma(x^{\prime},t)+\partial_{t}u(x^{\prime},\gamma(x^{\prime},t),t)=0$

and therefore

$|\nabla u(0,0)|\,\partial_{t}\gamma(0,0)+\partial_{t}u(0,0)=0.$

From this, (1.1) and (3.3), omitting the evaluation at $(0,0)$ for short, we
arrive at

$\displaystyle-\partial_{t}\gamma$ $\displaystyle=$
$\displaystyle\frac{\partial_{t}u}{|\nabla u|}$ $\displaystyle=$
$\displaystyle\frac{c\Delta u}{|\nabla u|}+\frac{1}{|\nabla
u|}\int_{\mathds{R}^{n}}\big{(}u(y,t)-\Theta(y,t)\big{)}_{+}\,K(x,y)\,dy+\left(\left(\omega+\frac{\beta(u)\nabla
u}{|\nabla u|^{\gamma}}\right)\cdot\frac{\nabla u}{|\nabla u|}\right)_{-}$
$\displaystyle=$ $\displaystyle-cH+\frac{c\partial_{x_{n}}^{2}u}{|\nabla
u|}+\frac{1}{|\nabla
u|}\int_{\mathds{R}^{n}}\big{(}u(y,t)-\Theta(y,t)\big{)}_{+}\,K(x,y)\,dyx$
$\displaystyle\qquad\qquad+\left(\left(\omega+\frac{\beta(u)\nabla u}{|\nabla
u|^{\gamma}}\right)\cdot\frac{\nabla u}{|\nabla u|}\right)_{-}\,,$

which is (3.1).

##  4\. Conclusions

We have put forth a model for bushfires relying on partial differential
equations. In its simplest form, the model takes the form

$\begin{split}\partial_{t}u&=c\Delta
u+\int_{\mathds{R}^{n}}\big{(}u(y,t)-\Theta(y,t)\big{)}_{+}\,K(x,y)\,dy+\left(\left(\omega+\frac{\beta(u)\nabla
u}{|\nabla u|^{\gamma}}\right)\cdot\nabla u\right)_{-}\,,\end{split}$

where $u$ represents the temperature of the environment, $c$ the environmental
diffusion coefficient, $\Theta$ the effective ignition temperature, $\omega$ a
parameter accounting for wind, and $\beta$ a function modulating a convective
effect.

The front of the fire can be modeled by the level set of the temperature
function corresponding to the ignition temperature. The velocity of the front
in the normal direction takes the form

$-cH+\frac{c\partial^{2}_{\nu}u}{|\nabla u|}+\frac{1}{|\nabla
u|}\int_{\mathds{R}^{n}}\big{(}u(y,t)-\Theta_{0}\big{)}_{+}\,K(x,y)\,dy+\left(\left(\omega+\frac{\beta(u)\nabla
u}{|\nabla u|^{\gamma}}\right)\cdot\nabla u\right)_{-}\,,$

where $H$ is the mean curvature of the burning region.

Specific features of this setting include:

  * •

The basic model can be easily generalized to comprise additional phenomena,
including a memory effect to identify regions which have been already burned,
variations of the environment in space and time, different types of heat
diffusions along the ground and in the air, etc.

  * •

Albeit non-standard, the partial differential equation is endowed with a solid
mathematical theory,

  * •

The description of the geometric front combines a smoothing effect, encoded by
the mean curvature term (with a positive sign, allowing for the descriptions
of moving fronts starting, for instance, by domains with corners),

  * •

The model can be efficiently implemented through numerical simulations,

  * •

The simulations run in real-time and are in agreement both with the data
collected in supervised lab experiments and with those coming from real
bushfire events.

## References