# A new variational discretization technique for initial value problems
bypassing governing equations

ARAlexander Rothkopf  JNJan Nordström  Faculty of Science and Technology,
University of Stavanger, 4021, Stavanger, Norway  Department of Mathematics,
Applied Mathematics, Linköping University, SE-581 83, Linköping, Sweden
Department of Mathematics and Applied Mathematics, University of Johannesburg,
P.O. Box 524, Auckland Park 2006, Johannesburg, South Africa

###### Abstract

Motivated by the fact that both the classical and quantum description of
nature rest on causality and a variational principle, we develop a novel and
highly versatile discretization prescription for classical initial value
problems (IVPs). It is based on an optimization (action) functional with
doubled degrees of freedom, which is discretized using a single regularized
summation-by-parts (SBP) operator. Formulated as optimization task it allows
us to obtain classical trajectories without the need to derive an equation of
motion. The novel regularization we develop in this context is inspired by the
weak imposition of initial data, often deployed in the modern treatment of
IVPs and is implemented using affine coordinates. We demonstrate numerically
the stability, accuracy and convergence properties of our approach in systems
with classical equations of motion featuring both first and second order
derivatives in time. onvergence properties of our approach in systems with
classical equations of motion featuring both first and second order
derivatives in time.

Initial Value Problem, Summation By Parts, Variational Principle,

###### keywords:

Research Article - PreprintFP

##  1 Introduction

The numerical treatment of dynamical phenomena in classical and quantum
systems is at the core of progress in natural sciences and engineering. In
computational fluid dynamics blazek2015computational  or electrodynamics
taflove2005computational , a set of coupled partial differential equations is
solved on a predefined geometric domain with boundary conditions, starting
from an initial condition, in order to predict trajectories of point particles
or configurations of fields. In the study of atomic properties, linear and
non-linear variants of the Schrödinger equation or more generally Lindblad
equations breuer2002theory  of multiple entangled particles are solved as
initial value problems. For an understanding of the nuclei of atoms on the
other hand, an ensemble of fluctuating quantum fields of a non-linear variant
of Maxwell’s equations (Yang-Mills theory) needs to be simulated on a
hypercubic grid (lattice QCD) gattringer2009quantum .

Much progress has been made in developing accurate and cost effective
discretization schemes for partial differential equations over the past two
decades. Due to their ease of implementation, finite difference schemes have
long enjoyed popularity, but historically were challenged when confronted with
intricate simulation geometries. It took the development of summation-by-parts
(SBP) finite difference operators (for reviews see e.g. svard2014review ;
fernandez2014review ; lundquist2014sbp ), to elevate finite difference schemes
to a similar level of versatility as traditional functional basis approaches,
such as Galerkin schemes nordstrom2017roadmap . The SBP approach both in
spatial dimensions, as well as in time lundquist2014sbp ;
nordstrom2013summation ; nordstrom2016summation  provides proofs of stability
for finite difference based discretization schemes via the so-called energy
method and is easily extended to higher order approximations.

Implementing the integration-by-parts property of the underlying continuum
IVP, summation-by-parts operators are an example of so-called mimetic
discretizations. It has been shown that SBP operators form a versatile
framework, which encompasses various other numerical approximation techniques
besides finite differences svard_stable_2007 , such as finite volume schemes
nordstrom_finite_2003 , spectral element carpenter_entropy_2014 , flux
reconstruction ranocha_summation-by-parts_2016  and both continuous
abgrall_analysis_2020  and discontinuous Galerkin (dG) gassner_skew-
symmetric_2013 ; hesthaven_stable_1996  schemes.

A crucial development in the numerical treatment of differential equations is
the concept of weak boundary or initial conditions. It acknowledges that the
solution of a discretized PDE not only in the interior of the domain but also
on the boundary ( or initial time slice ) need only be as accurate as the
order of the discretization. By allowing the solution to deviate from the
initial or boundary conditions within the tolerance of the discretization, one
obtains a new lever, which one can exploit in the construction of
discretization schemes. The simultaneous approximation term (SAT) approach
carpenter1994time  e.g. implements weak boundary or initial conditions by the
addition of appropriately designed penalty terms to the differential equation
of interest. In recent studies it has been shown how to absorb part of these
penalty terms into a redefinition of the SBP operators, in order to reduce
their null-space to the corresponding physical dimension, leading to so called
null-space consistent SBP operators svard2019convergence ;
linders_properties_2020 ; svard_convergence_2021 ; ranocha2021new .

In spite of this substantial progress in the numerical treatment of IVPs,
challenges of both conceptual and technical nature remain. The treatment of
intrinsic constraints, e.g. the divergence constraint for Maxwell’s equations
and the discretization of second order systems are two examples. The third one
is the derivation of the equations of motion of linear electrodynamics or non-
linear Yang-Mills theory in terms of the gauge potentials (see e.g.
Ipp:2018hai ). It is relevant for quantum theory and depends on a choice of
gauge. The choice of setting the zeroth component of the four-potential to
zero $A_{0}=0$ renders the role of Gauss’ law opaque, since $A_{0}$ is
actually the Lagrange multiplier that preserves this intrinsic constraint. On
the other hand, the discretized Lagrangian of these theories remains
manifestly gauge invariant and $A_{0}$ retains its central role. If one could
solve the associated initial value problem, i.e. determine future field
configurations directly on the level of the action, without the need to derive
the equation of motion, no choice of gauge is necessary and the manifested
gauge invariance would render Gauss’ law automatically fulfilled.

In the treatment of initial value problems for second order ODEs with the SBP-
SAT technique it was found that using the same regularized SBP operator for
the first and second derivative does not lead to stable procedures. Instead,
different SBP operators for the position and velocity degree of freedom had to
be defined nordstrom2016summation . On the other hand in the action
formulation of second order systems, at most first order derivatives act on
the degrees of freedom arnold_mathematical_1989 . This reduction of the order
of the derivatives compared to the corresponding governing equations by at
least one power is a general feature of the action formulation. A lower order
derivative operator is more robust against rough data than a higher order one,
which requires more regularity gustafsson2007high . By expressing the theory
directly on the level of the action one can reduce the highest order of
derivatives and in particular for second order systems, a single modified SBP
operator suffices to obtain a regularized and null-space consistent
discretization.

In this paper we develop a discretization scheme for one-dimensional initial
value problems, based on a generalized variational principle applied directly
to the action of the system of interest. I.e. we will determine the classical
trajectory of the system, given an initial condition, without the need to
derive an equation of motion. To this end we deploy the SBP technique to
approximate derivatives occurring in the continuum formulation of the problem
and will take inspiration by the SAT approach to regularize the resulting
difference operators. This is achieved by the use of affine coordinates to
absorb the whole penalty term including data into a redefined null-space
consistent SBP operator.

Of course all examples treated in this study have well known ODEs as equation
of motion, which can be solved with established numerical methods. As proof of
principle, the intention of this study is merely to establish the validity of
the direct action based approach and prepare the ground for its application to
higher dimensional theories in future work, where its utility is more obvious.
Our direct action approach allows us to formulate the system using any
geometry for which SBP operators have been developed and therefore allows us
to avoid the intricate construction of dual meshes elcott_building_2005 ,
which underlie implementations of other variational approaches to IVPs, such
as discrete exterior calculus stern_geometric_2015 .

It is important to note that the main point of this paper is to present a new
technique for solving IVPs without utilizing governing equations.

The paper is structured as follows: In section 2 we review the continuum
formulation of the conventional variational principle of classical physics for
initial value problems, featuring equations of motion with second order
derivatives in time. In the subsequent section 3 we introduce our
discretization prescription, discuss the need for regularization and construct
a regularized SBP operator based on initial value data. To extend the
applicability of our discretization scheme to a more general class of systems,
we consider a generalized variational principle in section 4 and show that our
approach successfully captures differential equations of motion containing
also first order derivatives. We close with a brief summary in section 5.

##  2 Continuum formalism for second order equations of motion

The classical physics of closed systems (i.e. systems that are not in contact
with their environment) is conveniently captured via their Lagrangian. In
point mechanics the Lagrangian is a functional, which depends on the
trajectory of the point mass $x(t)$ and its velocity $\dot{x}(t)$. In a field
theory, such as in electromagnetism, it is formulated in terms of the vector
fields $A_{\mu}(x)$ and their derivatives $\partial_{\nu}A_{\mu}(x)$. In the
simple systems under consideration here, the Lagrangian can be written as the
difference between the kinetic energy of the system and its potential energy.
Taking a point mass in a constant gravitational field as an explicit example
we have

$\displaystyle{\cal L}=T-V=\frac{1}{2}m\dot{x}^{2}(t)-mgx(t).$ (1)

In the 19th century, physicists discovered that the trajectory of a classical
particle can be obtained from a variational principle (see e.g.
ref.goldstein_classical_2001 ). If a point mass starts out at position
$x(t_{1})$ at $t_{1}$ and ends up at position $x(t_{2})$ at $t_{2}$ then the
classical trajectory that is realized in nature between those points is given
by the critical point of the classical action111More specifically, for short
times, where the classical trajectory has not yet reached any turning point,
the action exhibits an actual minimum. In the presence of turning points in
the classical trajectory it in general represents a saddle point of the
action. (see e.g. ref. arnold_mathematical_1989 )

$\displaystyle S[x(t),\dot{x}(t)]=\int_{t_{1}}^{t_{2}}dt\,{\cal
L}[x(t),\dot{x}(t)].$ (2)

This observation is stated as a boundary value problem, where the start and
end point of the trajectory are specified. While it is of conceptual interest
that such a formulation exists, we cannot use it to determine the classical
trajectory itself in a causal fashion, since in order to formulate the
variational principle, we already need to know where the point mass will end
up at $t_{2}$. When setting up an experiment, we are of course only in control
of the initial position and velocity.

To proceed, physicists conventionally convert the above boundary value problem
into an initial value problem using the following strategy: one derives a set
of differential equations that are equivalent to the variational principle and
which can be solved as initial value problem. These are the celebrated Euler-
Lagrange equations. We wish to inspect the variation of the Lagrangian using a
slightly deformed path $x(t)+\delta x(t)$. Here $\delta x(t)$ is an arbitrary
function with the only condition that it goes to zero at $t_{1}$ and $t_{2}$,
as the points $x(t_{1})$ and $x(t_{2})$ are fixed. Varying the action we
obtain

$\displaystyle\delta{\cal S}$
$\displaystyle=\int_{t_{1}}^{t_{2}}dt\Big{\\{}\frac{\partial{\cal L}}{\partial
x}\delta x+\frac{\partial{\cal
L}}{\partial\dot{x}}\delta\dot{x}\Big{\\}}=\int_{t_{1}}^{t_{2}}dt\Big{\\{}\frac{\partial{\cal
L}}{\partial x}\delta x+\frac{\partial{\cal
L}}{\partial\dot{x}}\frac{d}{dt}\delta x\Big{\\}}$ (3)
$\displaystyle=\int_{t_{1}}^{t_{2}}dt\Big{\\{}\frac{\partial{\cal L}}{\partial
x}-\frac{d}{dt}\frac{\partial{\cal L}}{\partial\dot{x}}\Big{\\}}\delta
x+\left.\Big{[}\frac{\partial{\cal L}}{\partial\dot{x}}\delta
x\Big{]}\right|_{t_{1}}^{t_{2}},$ (4)

where the last line results from integration by parts (IBP). Since the
variation $\delta x$ by construction vanishes on the boundary, the term in the
square brackets also vanishes. If we inspect the critical point of the
functional ${\cal S}$, defined by $\delta{\cal S}=0$, we find that it is
equivalent to the term in the curly brackets equalling zero, since $\delta x$
can be any (well behaved) function between $t_{1}$ and $t_{2}$.

Figure 1: Differences between the variational principle as boundary value
problem (left) and as initial value problem (right). In the BVP setting there
may exist multiple paths that fulfill the boundary conditions, but only one of
them, $x_{\rm cl}(t)$, represents an extremum of the action. In the IVP
setting, a doubling of the degrees of freedom is required in order to allow
the path at time $t_{2}$ to vary freely. Figure adapted from ref.
galley_classical_2013 .

In other words, if we assume the validity of the variational principle, i.e.
that the classical trajectory follows from the critical point of the action,
then this trajectory must fulfill the Euler-Lagrange equations, which are just
the terms inside the curly brackets set to zero

$\displaystyle\left.\frac{\delta{\cal S}[x,\dot{x}]}{\delta
x}\right|_{x=x_{\rm classical}}=0\quad\underset{\rm IVP}{\overset{\rm
BVP}{\iff}}\quad\frac{\partial{\cal L}}{\partial
x}-\frac{d}{dt}\frac{\partial{\cal L}}{\partial\dot{x}}=0.$ (5)

Since for more complex systems (with internal constraints etc.) it is often
easier to formulate the action than to derive the Euler-Lagrange equations as
initial value problem, our goal here is to formulate and solve the initial
value problem as a variational problem directly on the level of the action. To
this end we follow the reasoning of ref. galley_classical_2013 , which
establishes the continuum formalism for the variational IVP.

Retracing the train of thought of ref.galley_classical_2013 , we first note
that the equivalence between the Euler-Lagrange equations (which provide the
correct classical equations of motion) and the stationarity of the action
requires that the variation of the path vanishes at initial $t_{1}$ and final
time $t_{2}$ (see eq. 4). Since we do not know $x(t_{2})$ apriori and we only
know $x(t_{1})$ and $\dot{x}(t_{1})$, we must instead find a way how to
formulate the variational principle in such a way that the value of $x$ at
$t_{2}$ does not need to be fixed. This can be accomplished by doubling the
degrees of freedom with one trajectory $x_{1}(t)$ describing a forward path
and one trajectory $x_{2}(t)$ describing a backward path. Intuitively we will
use the doubled degrees of freedom in a way reminiscent of the shooting
method, extended such that the aim is to find the trajectory which returns to
the starting point given by the initial conditions (see fig. 1 for a sketch of
the difference in the approaches).

In order to achieve the necessary cancellation of the boundary terms between
the forward and backward path at $t_{2}$, ref.galley_classical_2013
constructs a new joint action for the two degrees of freedom $x_{1}(t)$ and
$x_{2}(t)$ as

$\displaystyle S_{\rm IVP}[x_{1}(t),\dot{x}_{1}(t),x_{2}(t),\dot{x}_{2}(t)]$
$\displaystyle=\int_{t_{1}}^{t_{2}}dt\Big{(}{\cal
L}[x_{1}(t),\dot{x}_{1}(t)]-{\cal L}[x_{2}(t),\dot{x}_{2}(t)]\Big{)},$ (6)
$\displaystyle=\int_{t_{1}}^{t_{2}}dt\Big{(}{\rm
L}[x_{1}(t),\dot{x}_{1}(t),x_{2}(t),\dot{x}_{2}(t)]\Big{)}.$ (7)

The Lagrangian housing the backward path $x_{2}$ is introduced with a relative
minus sign, which, as we will show, allows the boundary terms arising in the
variation of ${\cal L}[x_{1}(t),\dot{x}_{1}(t)]$ and ${\cal
L}[x_{2}(t),\dot{x}_{2}(t)]$ to cancel. Let’s carry out the variation of this
new $S_{\rm IVP}$ explicitly, which yields twice as many terms

$\displaystyle\delta{\cal S}$ $\displaystyle=\int
dt\Big{(}\Big{\\{}\frac{\partial{\rm L}}{\partial
x_{1}}-\frac{d}{dt}\frac{\partial{\rm L}}{\partial\dot{x}_{1}}\Big{\\}}\delta
x_{1}-\Big{\\{}\frac{\partial{\rm L}}{\partial
x_{2}}-\frac{d}{dt}\frac{\partial{\rm L}}{\partial\dot{x}_{2}}\Big{\\}}\delta
x_{2}\Big{)}$ (8) $\displaystyle+\left.\Big{[}\frac{\partial{\rm
L}}{\partial\dot{x}_{1}}\delta
x_{1}\Big{]}\right|_{t_{1}}^{t_{2}}-\left.\Big{[}\frac{\partial{\rm
L}}{\partial\dot{x}_{2}}\delta x_{2}\Big{]}\right|_{t_{1}}^{t_{2}}.$ (9)

In order to see how the cancellations come about, it is advantageous to change
coordinates, going over to relative $x_{-}=x_{1}-x_{2}$ and centered
coordinates $x_{+}=(x_{1}+x_{2})/2$. This change is not necessary, but
expressed in $x_{-}$ and $x_{+}$ the new variational principle can be
formulated in a very concise form and the relation between the functional and
the resulting differential equation of motion becomes much more lucid.

We vary the action using $x_{\pm}(t)+\delta x_{\pm}(t)$. The new path
deformations $\delta x_{\pm}(t)$ vanish at the initial time $t_{1}$, as the
original deformations are set to zero there $\delta x_{1}(t_{1})=\delta
x_{2}(t_{1})=0$. As the action is now a functional of the newly introduced
paths ${\cal S}_{\rm IVP}[x_{+}(t),\dot{x}_{+}(t),x_{-}(t),\dot{x}_{-}(t)]$,
its variation produces the following expression

$\displaystyle\delta{\cal S}_{\rm IVP}$ $\displaystyle=\int
dt\Big{(}\Big{\\{}\frac{\partial{\rm L}}{\partial
x_{+}}-\frac{d}{dt}\frac{\partial{\rm L}}{\partial\dot{x}_{+}}\Big{\\}}\delta
x_{+}+\Big{\\{}\frac{\partial{\rm L}}{\partial
x_{-}}-\frac{d}{dt}\frac{\partial{\rm L}}{\partial\dot{x}_{-}}\Big{\\}}\delta
x_{-}\Big{)}$ $\displaystyle+\left.\Big{[}\frac{\delta{\rm
L}}{\delta\dot{x}_{+}}\delta x_{+}+\frac{\delta{\rm
L}}{\delta\dot{x}_{-}}\delta x_{-}\Big{]}\right|_{t_{1}}^{t_{2}}.$ (10)

In order to correctly cancel the boundary contribution $\frac{\delta{\rm
L}}{\delta\dot{x}_{-}}\delta x_{-}$ at $t_{2}$, we see that the values of
$x_{1}(t_{2})$ and $x_{2}(t_{2})$ have to agree, i.e. $x_{-}(t_{2})=0$. It is
important to note that the paths $x_{1}$ and $x_{2}$ themselves are not fixed
to a certain value at $t_{2}$, since we do not know that value apriori. I.e.
the forward and backward paths need to be connected, corresponding to the
condition $x_{1}(t_{2})=x_{2}(t_{2})$.

What happens to the other boundary term $\frac{\delta{\rm
L}}{\delta\dot{x}_{+}}\delta x_{+}$? Since
$\dot{x}_{1}(t)=\dot{x}_{+}(t)+\frac{1}{2}\dot{x}_{-}(t)$ and
$\dot{x}_{2}(t)=\dot{x}_{+}(t)-\frac{1}{2}\dot{x}_{-}(t)$, we find the
following expression for the derivative of the joint Lagrangian ${\rm L}$

$\displaystyle\frac{\delta{\rm L}}{\delta\dot{x}_{+}}=\frac{\delta{\rm
L}}{\delta\dot{x}_{1}}\frac{\partial\dot{x}_{1}}{\partial\dot{x}_{+}}+\frac{\delta{\rm
L}}{\delta\dot{x}_{2}}\frac{\partial\dot{x}_{2}}{\partial\dot{x}_{+}}=\frac{\delta{\cal
L}}{\delta\dot{x}_{1}}-\frac{\delta{\cal
L}}{\delta\dot{x}_{2}}=\pi_{1}-\pi_{2}.$ (11)

In the second equality we have explicitly written ${\rm L}$ as the difference
between the individual Lagrangians ${\cal L}$ for the forward path $x_{1}$ and
backward path $x_{2}$. In the last step we furthermore introduced the
conjugate momenta of the paths, defined as $\pi_{1,2}=\delta{\cal
L}/\delta\dot{x}_{1,2}$. This relation between the functional derivative with
respect to $x_{+}$ and the difference between the momenta on the forward and
backward path tells us that we can make the remaining boundary term in eq. 10
vanish at $t_{2}$ if we construct our paths such that the difference between
the momenta $\pi_{1}(t_{2})-\pi_{2}(t_{2})=0$ vanishes at time $t_{2}$.

For the systems considered here, which exhibit second order derivatives in
their equation of motion, we have the kinetic term $T=\frac{1}{2}\dot{x}^{2}$
in the Lagrangian ${\cal L}$. This term leads to the identification
$\pi_{1,2}=\dot{x}_{1,2}$. In turn we find that if we require that in addition
to the values of the paths at $t_{2}$ also the derivatives are identified
$\dot{x}_{1}(t_{2})=\dot{x}_{2}(t_{2})$, both boundary terms in eq. 10 will
vanish. This establishes the necessary conditions for joining the forward and
backward path

$\displaystyle
x_{1}(t_{2})=x_{2}(t_{2}),\qquad\dot{x}_{1}(t_{2})=\dot{x}_{2}(t_{2}),$ (12)

in order to relate the extremum of the joint functional ${\rm L}$ to the
Euler-Lagrange equation expressions in the curly brackets in eq. 10.

We had to introduce doubled degrees of freedom to correctly cancel the
boundary terms that arise from the fact that for an IVP the value of the
classical path is unknown at time $t_{2}$. In the end there only exists a
single classical trajectory and we hence must undo the proliferation of
degrees of freedom. To this end ref. galley_classical_2013  introduces what
they call the physical limit, which enforces $x_{1}(t)-x_{2}(t)=x_{-}(t)=0$ at
all times. When applied to the equations of motion resulting from eq. 10, i.e.

$\displaystyle\frac{\partial L}{\partial x_{\pm}}=\frac{d}{dt}\frac{\partial
L}{\partial\dot{x}_{\pm}}$ (13)

only those equations independent of $x_{-}$ survive. Since the functional
$L={\cal L}[x_{1},\dot{x}_{1}]-{\cal L}[x_{2},\dot{x}_{2}]$ is constructed
from a difference of the Lagrangians on the forward and backward paths it will
always contain at least a linear dependence on $x_{-}$ and $\dot{x}_{-}$. Thus
only the equation in eq. 13, in which the derivative with respect to $x_{-}$
is taken can survive.

Combining the variation of the joint action of the forward and backward path
with the physical limit, we thus arrive at the following concise formulation
of the variational principle for a classical initial value problem

$\displaystyle\left.\frac{\delta S_{\rm IVP}[x_{\pm}]}{\delta
x_{-}}\right|_{x_{-}=0,x_{+}=x_{\rm class}}=0.$ (14)

Note that in deriving eq. 14, integration by parts (IBP) took center stage.
This fact motivates the use of summation-by-parts (SBP) operators in the
discretization of the variational principle in the next section 3. Indeed, if
the discretization is able to exactly mimic IBP, all steps up to this point
follow through also in the discrete setting (see eq. 38 in section 4.1).

Formulating classical mechanics as variational problem offers further insight
derived from Noether’s theorem. Following ref. sieberer2016keldysh  one can
show that Noether’s theorem for an action with doubled degrees of freedom can
be established and it provides two important results. Using as starting point
the action ${\cal S}_{\rm IVP}$ and using only integration by parts and the
swapping of differentiation and variation, it follows that the sum of the
energy of the forward and the backward path is preserved in time, as is the
difference between the two. This establishes that even though the additional
backward path has been added to the system the energy associated with it
remains bounded and the system is in fact stable. For more details and
discussion see A.

Take as explicit example the point mass in a constant gravitational field. Its
Lagrangian is ${\cal L}=\frac{1}{2}m\dot{x}^{2}(t)-mgx(t)$ and the Euler-
Lagrange equation reads

$\displaystyle\ddot{x}_{\rm class}(t)=-g,\qquad x_{\rm
class}(t)=-\frac{1}{2}gt^{2}+\dot{x}(0)t+x(0),$ (15)

which is nothing but Newtons law in terms of acceleration and can be solved in
a straight forward manner. We will take $g$ to be positive to indicate that
gravity is acting downwards.

Using the formalism based on the doubled degrees of freedom we have instead

$\displaystyle S_{\rm IVP}$ $\displaystyle=\int
dt\,\Big{(}\frac{1}{2}m(\dot{x}^{2}_{1}(t)-\dot{x}^{2}_{2}(t))-mg(x_{1}(t)-x_{2}(t))\Big{)},$
(16) $\displaystyle=\int
dt\Big{(}m\dot{x}_{+}(t)\dot{x}_{-}(t)-mgx_{-}(t)\Big{)}.$ (17)

In computing the variation of the action, we carried out one integration by
parts, which, in effect, allows us to re-express $\delta{\cal S}$ as depending
solely on the variation of the paths and not on their derivatives (see eq. 3).
Similarly we can integrate by parts here to move the time derivative on
$x_{-}$ in the kinetic term to $x_{+}$

$\displaystyle S_{\rm IVP}=\int
dt\Big{(}-m\ddot{x}_{+}(t)x_{-}(t)-mgx_{-}(t)\Big{)}.$ (18)

Since we identify both the values and derivatives of the paths at $t_{2}$ no
boundary terms contribute. Taking the functional derivative of eq. 18 with
respect to $x_{-}$, setting the result to zero and identifying $x_{+}=x_{\rm
class}$ in the physical limit yields exactly the conventional Euler-Lagrange
equation

$\displaystyle\ddot{x}_{+}(t)=\ddot{x}_{\rm class}(t)=-g.$ (19)

We have by now seen how the continuum variational principle for IVPs is
derived and have acquired intuition in a simple system what form the joint
action ${\cal S}_{\rm IVP}$ takes on in terms of $x_{+}$ and $x_{-}$. In the
remainder of the paper we will only work on the level of the joint action and
not need to refer to the equation of motion anymore. Let us briefly mention
that the corresponding functional $L[x_{+},\dot{x}_{+},x_{-},\dot{x}_{-}]$ for
a large variety of systems with second order equations of motion of the form
$\ddot{x}+f(x)=0$ can be written as

$\displaystyle{\cal S}_{\rm IVP}=\int
dt\Big{(}m\dot{x}_{+}(t)\dot{x}_{-}(t)-f(x_{+})x_{-}(t)\Big{)}.$ (20)

In section 4, after having established the discrete formalism for systems with
a second order equation of motion in time, we will consider a generalized
variational principle also derived in ref. galley_classical_2013 , which will
allow us to extend the discrete treatment to systems with differential
equations of motion containing also single derivatives in time.

##  3 Variational IVP based on SBP operators in time

###  3.1 A naive SBP discretization of the model boundary value problem

As a first step, let us formulate the discretized variational principle in its
conventional form as boundary value problem. The point mass in a constant
gravitational field will again serve as an explicit example.

Our goal here is to discretize the action for the single forward path $x(t)$

$\displaystyle{\cal S}=\int
dt\,\Big{(}\frac{1}{2}m\dot{x}^{2}(t)-mgx(t)\Big{)},$ (21)

with Dirichlet boundary conditions $x(0)=x_{i}$ and $x(1)=x_{f}$, in order to
compute the classical trajectory at its critical point. To this end we
introduce the path ${\bf x}=(x(0),x(\Delta t),x(2\Delta t),\ldots)^{\rm T}$
resolved at $N_{t}$ points with time step $\Delta t=1/(N_{t}-1)$. The integral
can be approximated by a quadrature rule, whose particular form is captured in
a (diagonal) positive definite matrix $\mathds{H}$ and which defines an inner
product on discretized paths $({\bf x},{\bf x}^{\prime})={\bf x}^{\rm
T}\mathds{H}{\bf x}^{\prime}$.

Remember that the derivation of the continuum variational principle in eq. 4
and eq. 9 required us to carry out integration by parts. In order to guarantee
the equivalence between the stationarity of the action and the solution of the
Euler-Lagrange equation, we must therefore use a discretization that mimics
IBP in the discrete setting. Hence we discretize the time derivative with a
summation by parts operator $\mathds{D}=\mathds{H}^{-1}\mathds{Q}$, compatible
with our choice of $\mathds{H}$, where $\mathds{Q}^{\rm
T}+\mathds{Q}=\mathds{E}_{N}-\mathds{E}_{0}={\rm diag}[-1,0,\ldots,0,1]$.

The lowest order SBP discretization scheme SBP21 of order two in the interior
and order one on the boundary ensues when choosing the trapezoid rule for
integration

$\mathds{H}^{[2,1]}=\Delta t\left[\begin{array}[]{ccccc}1/2&&&&\\\ &1&&&\\\
&&\ddots&&\\\ &&&1&\\\
&&&&1/2\end{array}\right],\quad\mathds{D}^{[2,1]}=\frac{1}{2\Delta
t}\left[\begin{array}[]{ccccc}-2&2&&&\\\ -1&0&1&&\\\ &&\ddots&&\\\ &&-1&0&1\\\
&&&-2&2\end{array}\right].$

The next higher order SBP scheme SBP42 is fourth order in the interior and
second order on the boundary

$\displaystyle\mathds{H}^{[4,2]}=\Delta
t\left[\begin{array}[]{ccccccc}\frac{17}{48}&&&&&&\\\ &\frac{59}{48}&&&&&\\\
&&\frac{43}{48}&&&&\\\ &&&\frac{49}{48}&&&\\\ &&&&&1&\\\ &&&&&&\ddots\\\
\end{array}\right]$ (28)

$\displaystyle\mathds{D}^{[4,2]}=\frac{1}{\Delta
t}\left[\begin{array}[]{ccccccccc}-\frac{24}{17}&\frac{59}{34}&-\frac{4}{17}&-\frac{3}{34}&&&&&\\\
-\frac{1}{2}&0&\frac{1}{2}&0&&&&&\\\
\frac{4}{43}&-\frac{59}{86}&0&\frac{59}{86}&-\frac{4}{43}&&&&\\\
\frac{3}{98}&0&-\frac{59}{86}&0&\frac{32}{49}&-\frac{4}{49}&&&\\\
&&\frac{1}{12}&-\frac{2}{3}&0&\frac{2}{3}&-\frac{1}{12}&&\\\
&&&&&&&\ddots&\end{array}\right].$ (35)

In this section we will show explicit results based on the SBP21 operator and
include the outcomes from the SBP42 operator in our scaling tests.

The discretized action, on which the variational boundary value problem rests,
reads

$\displaystyle\mathds{S}_{\rm BVP}=\frac{1}{2}m\big{(}\mathds{D}{\bf
x}\big{)}^{\rm T}\mathds{H}\big{(}\mathds{D}{\bf x}\big{)}-mg\mathds{1}^{\rm
T}\mathds{H}{\bf x}+\lambda_{1}(x(0)-x_{i})+\lambda_{2}(x(1)-x_{f}).$ (36)

We have added two Lagrange multipliers that are treated as additional
dynamical degrees of freedom of our system and in turn enforce the boundary
conditions of the numerical solution. This procedure may appear to introduce
the boundary conditions strongly, however it does not amount to an apriori
replacement of $x(t_{1})$ and $x(t_{2})$ by $x_{i}$ and $x_{f}$. During the
procedure to locate the critical point of $\mathds{S}_{\rm BVP}$, we find that
the minimization algorithms approach the extremum of the functional globally.
I.e. the boundary conditions are fulfilled to machine precision for the actual
solution, while deviations are possible at intermediate steps.

Note that when adding Lagrange multipliers to introduce constraints to an
optimization functional, the relevant extremum may become a saddle point. If
only algorithms are available that locate the minima of a functional, one can
circumvent this issue by optimizing the norm of the gradient of the functional
instead, for which the saddle point becomes a local minimum. In this study we
use as a preconditioning step the gradient-free Nelder-Mead minimizer with a
subsequent application of the Newton method and the Interior Point method
222Standard implementations of the aforementioned algorithms in the NMinimize
and FindMinimum functions provided by the software Mathematica 12.0
Mathematica  have been used..

Having derived the continuum Euler-Lagrange equations for the point mass in a
constant gravitational field before in eq. 15, we compute the explicit
solution for the classical trajectory in the time interval $t\in[0,1]$, based
on a value of $m=1$, $g=1$ and initial conditions $x(0)=1$, $\dot{x}(0)=0.3$.
In that case the point mass reaches the position $x(1)=0.8$ at time $t_{2}=1$.
Supplying these values to eq. 36, we can solve for the extremum and, as shown
the left panel of fig. 2, obtain a solution (red dots) that recovers the
correct solution of the Euler-Lagrange equations (gray solid line).

Figure 2:  Using a SBP21 operator in time for the boundary value formulation
we find (left) the discretized path ${\bf x}$ (red circles) ($N_{t}=32$) that
optimize the functional eq. 36. The discretized path ${\bf x}_{1}$ (red
circles) and ${\bf x}_{2}$ (blue crosses) ($N_{t}=32$) that optimize the
functional eq. 37 for the initial value formulation, corresponding to the
discretized $S_{\rm IVP}$ is shown in the right panel. Note that only half of
the path elements reproduce the correct solution. Continuum solution of the
Euler-Lagrange equation $\ddot{x}(t)=-g$ is shown as solid gray line.

While we succeed in recovering the correct solution, this approach, as
mentioned before, is conceptually not satisfactory, since the formulation of
the BVP relied on information about $x(1)$ obtained from the prior solution of
the Euler-Lagrange equations as initial value problem.

###  3.2 A naive SBP discretization of the model initial value problem

Let us continue by turning our attention to discretizing the continuum
formulation of the variational principle for initial value problems, which is
based on two paths. Introducing discretized paths ${\bf
x}_{1}=(x_{1}(0),x_{1}(\Delta t),x_{1}(2\Delta t),\ldots)^{\rm T}$ and
correspondingly ${\bf x}_{2}$ and using the same symbols as before for the
integration $\mathds{H}$ and summation-by-parts difference operators
$\mathds{D}$, we arrive at the following action

$\displaystyle\mathds{S}_{\rm IVP}$
$\displaystyle=\Big{\\{}\frac{1}{2}(\mathds{D}{\bf x}_{1})^{\rm
T}\mathds{H}(\mathds{D}{\bf x}_{1})-g\mathds{1}^{\rm T}\mathds{H}{\bf
x}_{1}\Big{\\}}-\Big{\\{}\frac{1}{2}(\mathds{D}{\bf x}_{2})^{\rm
T}\mathds{H}(\mathds{D}{\bf x}_{2})-g\mathds{1}^{\rm T}\mathds{H}{\bf
x}_{2}\Big{\\}}$
$\displaystyle+\lambda_{1}(x_{1}(0)-x_{i})+\lambda_{2}((\mathds{D}{\bf
x}_{1})(0)-\dot{x}_{i})$
$\displaystyle+\lambda_{3}(x_{1}(N_{t})-x_{2}(N_{t}))+\lambda_{4}((\mathds{D}{\bf
x}_{1})(N_{t})-(\mathds{D}{\bf x}_{2})(N_{t})).$ (37)

Here we have introduced four Lagrange multipliers to both enforce the initial
conditions for position $x_{i}$ and derivative $\dot{x}_{i}$ of the forward
path $(\lambda_{1},\lambda_{2})$, as well as to enforce the correct
identification of the position and derivatives at the last point of the
forward and backward path $(\lambda_{3},\lambda_{4})$. We consider all
$\lambda_{i}$’s as dynamical degrees of freedom, such that the constraints are
enforced exactly on the final solution of the optimization problem, while
permitting deviations from the constraints at intermediate steps.

Before we continue to determine the optimal paths according to
$\mathds{S}_{\rm IVP}$, we show that this discrete functional yields the
correct equations of motion according to the stationarity condition eq. 14, if
SBP operators are used. Focusing on the terms in curly brackets in eq. 37, we
introduce the discretized ${\bf x}_{1}={\bf x}_{+}+(1/2){\bf x}_{-}$ and ${\bf
x}_{2}={\bf x}_{+}-(1/2){\bf x}_{-}$ so that

$\displaystyle\mathds{S}$ $\displaystyle=\frac{1}{2}(\mathds{D}{\bf
x}_{+})^{\rm T}\mathds{H}(\mathds{D}{\bf x}_{-})+\frac{1}{2}(\mathds{D}{\bf
x}_{-})^{\rm T}\mathds{H}(\mathds{D}{\bf x}_{+})-g\mathds{1}^{\rm
T}\mathds{H}{\bf x}_{-}$ $\displaystyle=(\mathds{D}{\bf x}_{+})^{\rm
T}\mathds{H}(\mathds{D}{\bf x}_{-})-g\mathds{1}^{\rm T}\mathds{H}{\bf x}_{-}$
$\displaystyle=-(\mathds{D}\mathds{D}{\bf x}_{+})^{\rm T}\mathds{H}{\bf
x}_{-}+(\mathds{D}{\bf x}_{+})^{\rm T}(\mathds{E}_{N}-\mathds{E}_{0}){\bf
x}_{-}-g\mathds{1}^{\rm T}\mathds{H}{\bf x}_{-}$ (38)

Here we have used the symmetry of $\mathds{H}$ to arrive at the second line
and explicitly exploited the SBP property of $\mathds{D}$ in the third line.
As we enforce the initial conditions and identify the forward and backward
path at the final time, both boundary terms involving ${\bf x}_{-}$ vanish.
Mimicking the continuous derivation, let us take the derivative of
$\mathds{S}$ with respect to the i-th component of the vector ${\bf x}_{-}$,
which yields the following expression

$\displaystyle\Big{(}-(\mathds{D}\mathds{D}{\bf x}_{+})^{\rm
T}-g\mathds{1}^{\rm T}\Big{)}\mathds{H}{\bf e}_{i}=0.$ (39)

Since $\mathds{H}$ is diagonal, eq. 39 establishes the discrete equation of
motion $(\mathds{D}\mathds{D}{\bf x}_{+})=-g\mathds{1}$, a faithful
representation of the continuum result $\ddot{x}(t)=-g$.

Let us continue to determining the optimal paths ${\bf x}_{1}$ and ${\bf
x}_{2}$ that correspond to the critical point of $\mathds{S}_{\rm IVP}$ using
$\mathds{H}^{[2,1]}$ and $\mathds{D}^{[2,1]}$ on $N_{t}=32$ discrete points,
we find the solution shown in the bottom panel of fig. 2. We plot the values
of the forward path as red circles, while those of the backward path are given
as blue crosses. Note that they lie on top of each other, which tells us that
the optimal solution fulfills the physical limit condition ${\bf x}_{1}={\bf
x}_{2}$.

On the other hand we also immediately see that only around half of the points
on each path agree with the correct solution from the Euler-Lagrange equations
(gray solid line). The other half lies significantly below the correct
solution, forming a highly oscillatory structure. For an even number of grid
points the last point of ${\bf x}_{1}$ and ${\bf x}_{2}$ lies on the correct
trajectory, while for an odd number of points, the path ends on the
oscillatory structure below.

We have identified the origin of these oscillatory structures to arise from
the particular structure of the null-space of the finite difference operator.
In the kinetic terms of $\mathds{S}_{\rm IVP}$ both $\mathds{D}$ and
$\mathds{D}^{\rm T}$ appear. The study of null-space consistency of the lowest
order SBP21 operator considered here, reveals that it contains exactly two
zero eigenvalues. The space of right eigenvectors of $\mathds{D}$, associated
with this doubly degenerate eigenvalue, is only one-dimensional. Both
eigenvectors are proportional to the constant function.

However when we study the form of the left eigenvectors of $\mathds{D}$, or
equivalently the right eigenvectors of $\mathds{D}^{\rm T}$, we find that
those projecting into the null space are not at all constant but highly
oscillatory, reminiscent of the so-called $\pi$-mode. An example of these
eigenvectors is shown in fig. 3.

Figure 3: Right (red circles) and left (blue crosses) eigenvector associated
with a zero eigenvalue in $\mathds{D}^{[2,1]}$, based on $N_{t}=32$ points.
Note the highly oscillatory character of the latter.

Such unphysical oscillatory solutions have recently been identified to also
interfere in determining the solutions of differential equations with non-
trivial boundary conditions in one- and multiple dimensions in
ranocha_discrete_2020 . In the context of the variational problem considered
here, the oscillatory solutions did not affect the solution when formulated as
a boundary value problem in section 3.1. The fixing of the boundary at $t_{2}$
in the BVP formulation apparently prevents the oscillatory solution. On the
other hand the IVP action eq. 37 clearly accommodates these oscillatory paths.

An accurate discretization scheme for the IVP system action must therefore be
able to avoid the appearance of unphysical oscillatory modes and several
strategies to do so have been explored in the literature. One class of
strategies consists of modifying the first order derivative operator
$\mathds{D}$ by adding higher order derivative operators to it. A conventional
SBP finite difference operator of order $p$ needs to fulfill the derivative
property $\mathds{D}^{(p)}{\bf x}^{n}=(n-1){\bf x}^{n-1}$ exactly only on
monomials up to order $n=p-1$. Therefore, adding a higher derivative operator
$\Delta t\,\mathds{D}^{(p+2)}$, scaled by the grid spacing does not affect
this property, as it annihilates all lower order monomials. In addition, this
correction term vanishes in the limit of taking $\Delta t\to 0$. In the
context of upwind schemes one e.g. adds the symmetric second derivative
operator to the SBP first derivative, turning it into an upwind derivative. If
one deals with complex functions one may instead add the symmetric second
derivative multiplied with the imaginary unit. This modification is known as
adding a Wilson term wilson_confinement_1974  in the physics
literature333Wilson derived that regularization after investigating the
Green’s function of the differential operator that defines the equation of
motion of the system. What he found is that in Fourier space the Green’s
function exhibits not only a pole corresponding to the physical trajectory,
but due to the finite grid spacing a second pole appeared at the end of the
Brillouin zone, which introduces exactly the oscillatory mode we observed
above..

Both of these approaches present challenges, which we wish to avoid here. By
turning the central stencil into an upwind stencil, we lose the symmetry of
the system, which adversely affects the accuracy of the solutions. Introducing
a purely imaginary modification on the other hand requires the difference
operator to act on complex functions to be meaningful. One may contemplate the
possibility to complexify the functions involved in the variational problem,
which while only cumbersome in the classical case will lead to conceptual
problems when trying to use the discretization in the context of quantum path
integrals (c.f. sign problem).

We therefore wish to explore a different route to remove the unphysical zero
modes of the operator $\mathds{D}$, taking inspiration from more recent works
on null-space consistent SBP operators, such as in refs. svard2019convergence
; linders_properties_2020 ; svard_convergence_2021 ; ranocha2021new . The
central ingredient in these approaches is to exploit the weak formulation of
boundary and initial conditions. Concretely, when boundary conditions are
enforced weakly via a penalty term, this penalty term can be partially
absorbed into the derivative operator to remove the zero modes of that
operator. On the level of differential equations, the strategy works as
follows. Consider the following IVP, the differential equation for the
exponential function

$\displaystyle\frac{d}{dt}u(x)=\lambda u(x),\quad u(0)=u_{0},$ (40)

which in its discretized form reads

$\displaystyle\mathds{D}{\bf u}=\lambda{\bf
u}+\sigma_{0}\mathds{H}^{-1}\mathds{E}_{0}\big{(}{\bf u}-{\bf g}\big{)}.$ (41)

Here we have added a so-called SAT penalty term on the RHS, which includes the
matrix $\mathds{E}_{0}={\rm diag}[1,0,\ldots,0]$ that singles out the first
entry in the discretized functions ${\bf u}$ and ${\bf g}$. The former ${\bf
u}$ refers to the solution of the differential equation and the latter ${\bf
g}=(u_{0},0,\ldots,0)$ contains the initial value as its first entry. Note
that $\mathds{H}^{-1}$ contains ${\Delta t}^{-1}$, which contributes with
increasing weight as ${\Delta t}\to 0$. The parameter $\sigma_{0}$ in the SBP-
SAT approach is tuned to satisfy stability properties and its optimal value is
found to be $\sigma_{0}=-1$, a choice we adopt in the following. The standard
approach developed in the conventional SBP-SAT treatment of IVPs consists of
absorbing the penalty term proportional to ${\bf u}$ into a redefined
$\tilde{\mathds{D}}=\mathds{D}-\sigma_{0}\mathds{H}^{-1}E_{0}$, which does not
feature any zero modes anymore. That operator is now non-singular
ruggiu_eigenvalue_2020  and may be inverted to obtain the solution ${\bf u}$.
In the next section we will develop a similar strategy applicable to the
variational problem.

###  3.3 Regularized SBP discretization of the model initial value problem

Taking inspiration from the work on regularizing SBP operators in differential
equations, we set out to absorb information about the initial conditions into
the SBP operator as means of regularization. In the functional of eq. 37 we do
not have an equality sign, such as in our example (41), to rearrange terms.
Instead we must find a way to incorporate the whole penalty term in
$\mathds{D}$. Note that the penalty term contains one expression that is
proportional to the function that the SBP operator acts on and one expression
proportional to a constant. I.e. we have to modify the difference operator to
include a shift. In other words, we are dealing with an affine transformation.

Figure 4: (Left) Eigenvalue spectrum of the unregularized SBP21 operator
$\mathds{D}^{[2,1]}$ (blue crosses) for $N_{t}=32$ grid points and
corresponding spectrum of the regularized operator $\bar{\mathds{D}}^{[2,1]}$
(red circles), which does not feature any zero modes. The zoomed inset on the
right reveals the presence of the zero modes in the unregularized SBP21
operator $\mathds{D}^{[2,1]}$.

There exists an elegant way to express affine transformations using so-called
affine coordinates. One defines $\bar{A}[{\bf b}]\bar{\bf x}=A{\bf x}+{\bf
b}$, where $\bar{A}[{\bf b}]$ refers to the matrix $A$ amended by one more row
and column with $1$ placed in the lower right corner. The additional column
available in $\bar{A}[{\bf b}]$ is filled with the values of ${\bf b}$. The
vector $\bar{\bf x}$ is just ${\bf x}$ amended by one more entry with value
one. For our application to the variational formulation of the IVP we
therefore define a new $\bar{\mathds{D}}$ using as shift the vector containing
the initial values ${\bf b}=\sigma_{0}\mathds{H}^{-1}E_{0}{\bf g}$ where ${\bf
g}={\rm diag}[x_{i},x_{i}+\Delta t\,\dot{x}_{i},0,\cdots,0]$. For the SBP21
operator444For a higher order SBP operator, the values of ${\bf g}$ need to be
chosen, so that $[\mathds{D}{\bf g}](0)=\dot{x}_{i}$ and $[{\bf g}](0)=x_{i}$.
the explicit expression we obtain reads

$\displaystyle\bar{\mathds{D}}^{[2,1]}=\left[\begin{array}[]{cccccc}-\frac{1}{\Delta
t}-\sigma_{0}\frac{2}{\Delta t}&\frac{1}{\Delta
t}&&&&\sigma_{0}\frac{2}{\Delta t}x_{i}\\\ -\frac{1}{2\Delta
t}&0&\frac{1}{2\Delta t}&&&0\\\ &&\ddots&&&\vdots\\\ &&-\frac{1}{2\Delta
t}&0&\frac{1}{2\Delta t}&0\\\ &&&-\frac{1}{\Delta t}&\frac{1}{\Delta t}&0\\\
0&&\ldots&&0&1\\\ \end{array}\right].$ (48)

Note that in this paper we choose the parameter $\sigma_{0}=-1$, whenever a
penalty term is incorporated in ${\bar{\mathds{D}}}$. This choice is motivated
by the fact that in the conventional treatment of IVPs using the SBP-SAT
approach, this value leads to a minimal discretization error (see e.g. ref.
lundquist2014sbp ). We find that $\sigma_{0}=-1$ is optimal for our approach
too, as only in this case the correct classical solution is recovered.

All zero modes of the original operator $\mathds{D}^{[2,1]}$ have been lifted
in $\bar{\mathds{D}}^{[2,1]}$ and the resulting spectrum of eigenvalues $\nu$
is shown in fig. 4. Note that $\bar{\mathds{D}}^{[2,1]}$ still correctly
annihilates the constant function, as long as it is compatible with the
initial conditions $x(0)=x_{i}$. In affine coordinates this annihilation does
not lead to a resulting zero vector, but a vector that contains vanishing
entries, except for the final one associated with the single real eigenvalue
of value one, shown in fig. 4.

Figure 5: Numerical solution of the discretized path ${\bf x}_{1}$ (red
circles) and ${\bf x}_{2}$ (blue crosses) ($N_{t}=32$) that optimize the
functional eq. 49 corresponding to the discretized $S_{\rm IVP}$. We use the
regularized SBP21 operator in time. Continuum solution of the Euler-Lagrange
equation $\ddot{x}(t)=-g$ is shown as solid gray line. Note that the solution
successfully avoids oscillatory contamination.

When formulating the action with the modified SBP operator, we obtain

$\displaystyle\mathds{S}_{\rm IVP}$
$\displaystyle=\Big{\\{}\frac{1}{2}(\bar{\mathds{D}}{\bar{\bf x}}_{1})^{\rm
T}\bar{\mathds{H}}(\bar{\mathds{D}}\bar{{\bf x}}_{1})-g\mathds{1}^{\rm
T}{\mathds{H}}{\bf
x}_{1}\Big{\\}}-\Big{\\{}\frac{1}{2}(\bar{\mathds{D}}\bar{{\bf x}}_{2})^{\rm
T}\bar{\mathds{H}}(\bar{\mathds{D}}\bar{{\bf x}}_{2})-g\mathds{1}^{\rm
T}\mathds{H}{\bf x}_{2}\Big{\\}}$
$\displaystyle+\lambda_{1}(x_{1}(0)-x_{i})+\lambda_{2}((\mathds{D}{\bf
x}_{1})(0)-\dot{x}_{i})$
$\displaystyle+\lambda_{3}(x_{1}(N_{t})-x_{2}(N_{t}))+\lambda_{4}((\mathds{D}{\bf
x}_{1})(N_{t})-(\mathds{D}{\bf x}_{2})(N_{t})).$ (49)

In order to implement the inner product in affine coordinates, we define
$\bar{\mathds{H}}$, which denotes the matrix $\mathds{H}$, amended by one
extra row and column of values zero. The last entry of the vector
$\bar{\mathds{D}}{\bar{\bf x}}_{1,2}$ serves only to implement the shift in
affine coordinates, hence it can be discarded via $\bar{\mathds{H}}$ since the
regularized SBP operator has already acted on the path.

Note that here we again add Lagrange multipliers as dynamical degrees of
freedom, to fulfill the initial conditions. One may ask whether enforcing the
initial conditions in this way neutralizes the effect of the regularization.
We emphasize that this is not the case. Minimization algorithms approach the
extremum of the functional globally, allowing the regulator to remain
effective and to avoid the oscillatory solutions.

Another question of both conceptual and practical relevance is whether the
functional eq. 49 houses one or multiple different local extrema. For the case
of the point particle in a constant gravitational field, we find that the
answer is that the functional is convex and a thus any local extremum is also
a global extremum. Let us determine the curvature of eq. 49 with respect to
the individual entries of the paths ${\bf x}_{1,2}$. The matrix
$A_{ij}=\partial^{2}\mathds{S}_{\rm IVP}/\partial x_{i}\partial
x_{j}=[\bar{\mathds{D}}^{\rm T}\bar{\mathds{H}}\bar{\mathds{D}}]_{ij}$ is
indeed positive semi-definite, as can be checked explicitly using a computer
algebra system. However, in general convexity is not automatic and needs to be
checked on a case-by-case basis555None of the functionals treated in this
study suffered from multiple extrema, allowing the Newton and Quasi-Newton
methods implemented in Mathematica Mathematica  to arrive at a single solution
independent of starting point..

For the point particle in a constant gravitational field, the solutions ${\bf
x}_{1}$ and ${\bf x}_{2}$ obtained with eq. 49 and the regularized SBP21
operator are shown in fig. 5 as red circles and blue crosses
respectively.666An explicit implementation of all examples discussed in this
manuscript can be found as open-access Mathematica script at the Zenodo
repository rothkopf_mathematica_2022 .. The regularization has successfully
removed the contamination by an unphysical oscillatory mode and we are able to
reproduce the correct classical solution.

The initial conditions, implemented using Lagrange multipliers and expressed
in a variational formulation, have provided us with a novel discretization
prescription for a wide range of classical systems whose differential
equations of motion contain second order derivatives. We emphasize that we did
not have to derive the equation of motion to compute the classical trajectory
here.

Let us take a look at the accuracy and convergence properties of the
discretization scheme constructed in this section. The simple model of a point
mass in a constant gravitational field again serves as explicit example. To
this end we compute the optimal path ${\bf x}_{1}$ according to the
appropriately regularized eq. 49 using both the regularized SBP21 and
regularized SBP42 operator on different grids using $N_{t}\in[16,\ldots,512]$
points. We compare the values of the path at the final step $t_{2}=1$ to the
analytically known solution and compute the absolute error between them. These
errors are shown in fig. 6.

Figure 6: Deviation of the value (left plot) and the derivative (right plot)
between the optimal numerical solution of eq. 49 and the true solution at the
final time $t_{2}=1$. The deviation for the regularized SBP21 operator is
given as blue crosses, the one for the SBP42 operator as red circles. Note
that while in the SBP21 case the solution improves steadily with diminishing
$\Delta t\sim 1/N_{t}$, the results for the SBP42 operator are already exact
to machine precision. For a more detailed view of the SBP21 behavior see fig.
7

The lowest order SBP21 approximation (blue crosses) exhibits steady
improvement in the residual deviation from the true solution, as the grid
spacing is reduced. In fig. 7 we zoom in on the SBP21 errors and fit them with
a power-law ansatz, which on the log-log plot appears as a straight line. The
best fit exponent $\Delta t^{2.03}$ we obtain, tells us that our
discretization scheme achieves second order accuracy in the solution values.
Interestingly, when considering the SBP42 operator, we find that irrespective
of the grid spacing we are able to reproduce the true solution down to machine
precision (which in our case, using the Mathematica software package, was set
to $10^{-30}$). This result is reassuring, as by construction the SBP operator
and the corresponding quadrature rule are able to differentiate and integrate
polynomials up to second order exactly. Since the solution of the point mass
in the constant gravitational field is a parabola, we do not find any residual
dependence on the grid spacing.

Figure 7: Deviation of the value between the optimal numerical solution of eq.
49 and the true solution at the final time $t_{2}=1$ (blue crosses) for the
SBP21 operator. The best power law fit to the behavior is shown as the gray
solid line which corresponds to $\Delta t^{2.03}$.

###  3.4 Discretization of non-linear second-order IVPs

While instructive, our model example described a very simple linear system.
Let us use the formalism established in eq. 20 to apply our discretization
prescription to a system, which features a genuinely non-linear differential
equation of motion instead. The action functional behind the ODE

$\displaystyle\ddot{x}(t)+\kappa x^{3}(t)=0,$ (50)

is given by the following expression

$\displaystyle{\cal S}_{\rm IVP}=\int
dt\Big{(}\dot{x}_{+}(t)\dot{x}_{-}(t)-\kappa x_{+}^{3}(t)x_{-}(t)\Big{)},$
(51)

which we must subsequently discretize. The resulting optimization functional
reads

$\displaystyle\mathds{S}_{\rm IVP}$
$\displaystyle=\Big{\\{}\frac{1}{2}(\bar{\mathds{D}}{\bar{\bf x}}_{1})^{\rm
T}\bar{\mathds{H}}(\bar{\mathds{D}}\bar{{\bf
x}}_{1})\Big{\\}}-\Big{\\{}\frac{1}{2}(\bar{\mathds{D}}\bar{{\bf x}}_{2})^{\rm
T}\bar{\mathds{H}}(\bar{\mathds{D}}\bar{{\bf x}}_{2})\Big{\\}}$
$\displaystyle-\kappa\Big{(}({{\bf x}}_{1}+{{\bf
x}}_{2})/2\Big{)}^{3}\mathds{H}({{\bf x}}_{1}-{{\bf x}}_{2})$
$\displaystyle+\lambda_{1}(x_{1}(0)-x_{i})+\lambda_{2}((\mathds{D}{\bf
x}_{1})(0)-\dot{x}_{i})$
$\displaystyle+\lambda_{3}(x_{1}(N_{t})-x_{2}(N_{t}))+\lambda_{4}((\mathds{D}{\bf
x}_{1})(N_{t})-(\mathds{D}{\bf x}_{2})(N_{t})).$ (52)

The third power in the second line is understood as acting element-wise on the
entries of $({{\bf x}}_{1}+{{\bf x}}_{2})/2$. The continuum trajectory for the
choice $\kappa=20$ is given as the gray solid line in fig. 8. We plot it
together with the numerical solutions ${\bf x}_{1}$ (red circles) and ${\bf
x}_{2}$ (blue crosses) based on the SBP21 operator along $N_{t}=32$ grid
points.

Figure 8: Numerical solution of the discretized path ${\bf x}_{1}$ (red
cricles) and ${\bf x}_{2}$ (blue crosses) ($N_{t}=32$) that optimize the
functional eq. 52 corresponding to the discretized $S_{\rm IVP}$. We use the
regularized SBP21 summation-by-parts operator in time. Continuum solution of
the Euler-Lagrange equation $\ddot{x}(t)+\kappa x^{3}(t)=0$ for $\kappa=20$ is
shown as solid gray line. Note that the solution successfully avoids
oscillatory contamination.

Let us consider how the new discretization prescription performs on this
second order non-linear problem. The continuum solution is given in terms of
the Jacobi elliptic function and its inverse, i.e. it is not polynomial. Thus
neither the SBP21 nor the SBP42 operator are able to reproduce it exactly. We
find minute oscillations around the true solution exhibited by the data points
in fig. 8, which diminish monotonously with grid refinement, as expected from
a stable procedure. The deviation of the numerical solution ${\bf x}_{1}$ from
the true solution is shown in the left panel of fig. 9 as blue crosses for the
SBP21 operator and as red circles for the SBP42 operator.

Figure 9: Deviation of the value (left plot) and the derivative (right plot)
between the optimal numerical solution of eq. 52 and the true solution at the
final time $t_{2}=1$. The deviation for the regularized SBP21 operator is
given as blue crosses with best fit $\Delta t^{2.12}$, the one for the
regularized SBP42 operator as red circles with $\Delta t^{3.35}$. The best
power law fit to the data is shown as gray lines. Both orders of the solution
improve steadily with diminishing $\Delta t\sim 1/N_{t}$. We find that our
approach exhibits convergence in the values with one higher power in the grid
spacing than for the derivative, where for SBP21 we find as best fit $\Delta
t^{1.06}$ and for SBP42 we have $\Delta t^{1.87}$.

One finds that compared to the linear system, the error made in the non-linear
system is around one order of magnitude larger at the same lattice spacing for
the SBP21 operator. However the order of convergence remains close to
quadratic with a best fit $\Delta t^{2.12}$ shown by the gray lines in fig. 9.
For the SBP42 operator, we find that the error at $N_{t}=16$ improves by one
order of magnitude and convergence to the continuum limit proceeds with
$\Delta t^{3.35}$. The expected behavior from solutions of ODEs
svard2019convergence ; svard_convergence_2021  in this case is $\Delta t^{3}$,
and while our method seems to perform better in this scenario, we believe that
to be a coincidence.

In line with established results for the conventional numerical treatment of
differential equations, we see in fig. 9 that the error in the derivative
$\dot{x}=[\mathds{D}{\bf x}_{1}]({N_{t}})$ of the numerical solution shows
convergence with one full order less than the values of the solution itself.
As plotted in the right panel of fig. 9 we obtain for the regularized SBP21
operator convergence for the derivative according to $\Delta t^{1.06}$, while
the SBP42 operator exhibits $\Delta t^{1.87}$.

Having established the applicability and convergence properties of our novel
discretization approach for both a linear and non-linear second order
differential equation of motion, let us continue to treat systems with
equations of motion that feature different orders of time derivatives.

##  4 Discretization of first order derivative terms

So far we have considered the simplest case of physical systems with classical
equations of motion that contain a second derivative in time. These follow
naturally from the conventional formulation of the continuum variational
principle, based on an action that is written in terms of a Lagrangian. As has
been shown in galley_classical_2013 , the variational principle is able to
accommodate a much larger variety of systems, including those with dissipative
forces, which are not time-reversal invariant. Such systems exhibit equations
of motion, which contain also first order derivatives in time. The crucial
step is to realize that, one may generalize the classical variational
principle by adding to the Lagrangian $L$ another functional $\Lambda$ that
may depend on both the forward and backward path777It is interesting to note
that the term $\Lambda$ also arises naturally in the classical limit of the
Schwinger-Keldysh contour formalism of the quantum path integral for
dissipative systems, after integrating out the environment degrees of freedom.
In that context it is known as the so called Feynman-Vernon influence
functional. and their derivatives as follows

$\displaystyle S_{\rm GIVP}[x_{1}(t),\dot{x}_{1}(t),x_{2}(t),\dot{x}_{2}(t)]=$
$\displaystyle\int_{t_{1}}^{t_{2}}dt\Big{(}{\cal
L}[x_{1}(t),\dot{x}_{1}(t)]-{\cal
L}[x_{2}(t),\dot{x}_{2}(t)]+\Lambda[x_{1}(t),\dot{x}_{1}(t),x_{2}(t),\dot{x}_{2}(t)]\Big{)}.$
(53)

Ref. galley_classical_2013 , aided by ref.berges_quantum_2007 , shows in
detail that the classical equations of motion also for this generalized
variational principle are obtained by going over to relative
$x_{-}=x_{1}-x_{2}$ and centered coordinates $x_{+}=(x_{1}+x_{2})/2$ with the
defining equation

$\displaystyle\left.\frac{\delta S_{\rm GIVP}[x_{\pm}]}{\delta
x_{-}}\right|_{x_{-}=0,x_{+}=x_{\rm class}}=0.$ (54)

The stability properties of these systems, as shown explicitly in ref.
sieberer2016keldysh , can also be formulated in terms of a generalized
Noether’s theorem in which e.g. the time-dependence of the total system energy
is correctly captured. This immediately invites us to apply the discretization
prescription developed in the previous section to two systems often considered
in the literature, the one which features the defining equation of the
exponential function as equation of motion, as well as the damped harmonic
oscillator. By considering these two examples, we acquire intuition in how to
construct the appropriate continuum functional $\Lambda$ of eq. 53, in order
to describe systems, which feature a differential equation of motion also
containing first order derivatives.

###  4.1 A purely first order system

Figure 10: Numerical solution of the discretized path ${\bf x}_{1}$ (red
circles) and ${\bf x}_{1}$ (blue crosses) ($N_{t}=32$) that optimize the
functional eq. 58 corresponding to the discretized $S_{\rm GIVP}$. We use the
regularized SBP21 operator in time. Continuum solution of the Euler-Lagrange
equation $\dot{x}(t)-\frac{5}{2}x(t)=0$ is shown as solid gray line.

Our goal here is to determine the classical trajectory of a system, which
features as its equation of motion the defining equation of the exponential
function

$\displaystyle\dot{x}(t)-\kappa x(t)=0.$ (55)

Since there are no second order derivatives present in eq. 55 we do not need
to supply the standard kinetic term to L in (53). The term linear in $x(t)$
can be thought of as arising from a potential contribution in ${\cal L}$,
which must contain one power of $x_{+}$ and one power of $x_{-}$ similar to
our argument in eq. 20. The new ingredient is the term that features a single
time derivative. It has to emerge from $S_{\rm GIVP}$ after functional
differentiation with respect to $x_{-}$. This behavior is achieved by choosing
the following Lagrangian and $\Lambda$ functional

$\displaystyle{\cal L}=-\frac{1}{2}\kappa
x^{2}(t),\quad\Lambda=\dot{x}_{+}x_{-},$ (56)

which amounts to the joint Lagrangian

$\displaystyle{\rm L}=-\kappa x_{+}(t)x_{-}(t)+\dot{x}_{+}x_{-}.$ (57)

Note that if one rewrites eq. 57 explicitly in terms of $x_{1}$ and $x_{2}$,
the contribution from $\Lambda$ indeed does not factorize into terms that
depend on $x_{1}$ or $x_{2}$ separately.

Figure 11: Deviation of the value (left plot) and the derivative (right plot)
between the optimal numerical solution of eq. 58 and the true solution at the
final time $t_{2}=1$. The deviation for the regularized SBP21 operator is
given as blue crosses with best fit $\Delta t^{2.03}$, the one for the
regularized SBP42 operator as red circles and $\Delta t^{2.95}$. The best
power law fit to the data is shown as gray lines. Both orders of the solution
improve steadily with diminishing $\Delta t\sim 1/N_{t}$. We find that our
approach exhibits convergence in the values with one higher power in the grid
spacing than for the derivative, where for SBP21 we find $\Delta t^{1.06}$ and
for SBP42 we have $\Delta t^{1.99}$.

Using the strategy developed in the previous section, let us write down the
discretized action functional, keeping in mind that for a first order equation
only the initial position needs to be supplied at the beginning of ${\bf
x}_{1}$. Correspondingly only the position information needs to be matched at
the end of the paths

$\displaystyle\mathds{S}_{\rm GIVP}$
$\displaystyle=\Big{\\{}-\frac{1}{2}\kappa({\bf x}_{1})^{\rm T}\mathds{H}{\bf
x}_{1}\Big{\\}}-\Big{\\{}-\frac{1}{2}\kappa({\bf x}_{2})^{\rm T}\mathds{H}{\bf
x}_{2}\Big{\\}}$
$\displaystyle+\frac{1}{2}\Big{(}\bar{\mathds{D}}\big{(}{\bar{\bf
x}}_{1}+{\bar{\bf x}}_{2}\big{)}\Big{)}^{\rm T}\bar{\mathds{H}}({\bf
x}_{1}-{\bf x}_{2})$
$\displaystyle+\lambda_{1}(x_{1}(0)-x_{i})+\lambda_{3}(x_{1}(N_{t})-x_{2}(N_{t})).$
(58)

The solutions ${\bf x}_{1}$ (red circles) and ${\bf x}_{2}$ (blue crosses) to
this equation on the interval $t\in[0,1]$ discretized with $N_{t}=32$
equidistant steps $\Delta t$ and a $\kappa=\frac{5}{2}$ with initial condition
$x_{i}=1$ produces the data shown in fig. 10.

In fig. 11 we plot the difference between the classical trajectory obtained
from eq. 58 and the true solution at the final time $t_{2}=1$ focussing on the
value itself in the left plot and the derivative in the right plot. The
results for the regularized SBP21 operator are given as blue crosses, those
for the SBP42 operator as red circles. Since the solution is not a simple
polynomial, the SBP42 operator cannot exactly integrate it. Both the SBP21 and
SBP42 cases show the same convergence rates, as observed in the conventional
formulation of IVPs (c.f. ref. svard_convergence_2021 ). We find again that
the convergence is one order higher in the values of the solution than in the
derivative of the solution. The SBP21 operator yields a $\Delta t^{2.03}$
improvement for the values of $x_{1}(1)$, while the SBP42 operator exhibits
$\Delta t^{2.95}$.

###  4.2 The damped harmonic oscillator

Figure 12: Numerical solution (red dots) of the discretized path ${\bf x}_{1}$
(red circles) and ${\bf x}_{2}$ (blue crosses) ($N_{t}=32$) that optimize the
functional eq. 62 corresponding to the discretized $S_{\rm GIVP}$. We use the
regularized SBP21 operator in time. Continuum solution of the Euler-Lagrange
equation $\mu\ddot{x}(t)+\xi\dot{x}(t)+\kappa x(t)=0$ is shown as solid gray
line.

Figure 13: Deviation of the value (left plot) and the derivative (right plot)
between the optimal numerical solution of eq. 62 and the true solution at the
final time $t_{2}=1$. The deviation for the regularized SBP21 operator is
given as blue crosses with best fit $\Delta t^{2.03}$, the one for the
regularized SBP42 operator as red circles with $\Delta t^{3.04}$. The best
power law fit to the data is shown as gray lines. Both orders of the solution
improve steadily with diminishing $\Delta t\sim 1/N_{t}$. We find again that
our approach exhibits convergence in the values with one higher power in the
grid spacing than for the derivative, where for SBP21 we find $\Delta
t^{1.01}$ and for SBP42 we have $\Delta t^{2.06}$.

As final challenge let us now turn to a physics system, which exhibits both
first and second order derivatives in its equation of motion: the damped
harmonic oscillator. The damped harmonic oscillator is characterized by an
Euler-Lagrange equation that reads

$\displaystyle\mu\ddot{x}(t)+\xi\dot{x}(t)+\kappa x(t)=0.$ (59)

Even though this system underlies a wealth of experimentally relevant
phenomena, the conventional formulation of classical mechanics is unable to
accommodate it in terms of a classical Lagrangian. In the generalized approach
of ref. galley_classical_2013  the kinetic and conservative force term are
captured by the Lagrangian and the dissipative term is included via the
$\Lambda$ functional. Remembering that only terms linear in $x_{-}$ survive
the stationarity condition and using the intuition we built in the preceding
sections, we can now write down the corresponding expressions for the
functionals of the damped harmonic oscillator

$\displaystyle{\cal L}=\frac{1}{2}\mu\dot{x}^{2}(t)-\kappa
x^{2}(t),\quad\Lambda=-\xi\dot{x}_{+}x_{-},$ (60)

which correspond to the joint Lagrangian

$\displaystyle{\rm L}=\mu\dot{x}_{+}(t)\dot{x}_{-}-2\kappa
x_{+}(t)x_{-}(t)-\xi\dot{x}_{+}x_{-}.$ (61)

Inserting the above into eq. 54 immediately yields eq. 59. The discretized
joint action functional based on eq. 61 reads

$\displaystyle\mathds{S}_{\rm GIVP}$
$\displaystyle=\Big{\\{}\frac{1}{2}\mu(\bar{\mathds{D}}{\bar{\bf x}}_{1})^{\rm
T}\bar{\mathds{H}}(\bar{\mathds{D}}\bar{{\bf x}}_{1})-\frac{1}{2}\kappa({\bf
x}_{1})^{\rm T}\mathds{H}{\bf
x}_{1}\Big{\\}}-\Big{\\{}\frac{1}{2}\mu(\bar{\mathds{D}}{\bar{\bf
x}}_{2})^{\rm T}\bar{\mathds{H}}(\bar{\mathds{D}}\bar{{\bf
x}}_{2})-\frac{1}{2}\kappa({\bf x}_{2})^{\rm T}\mathds{H}{\bf x}_{2}\Big{\\}}$
$\displaystyle-\xi\frac{1}{2}\Big{(}\bar{\mathds{D}}\big{(}{\bar{\bf
x}}_{1}+{\bar{\bf x}}_{2}\big{)}\Big{)}^{\rm T}\bar{\mathds{H}}({\bf
x}_{1}-{\bf x}_{2})$
$\displaystyle+\lambda_{1}(x_{1}(0)-x_{i})+\lambda_{2}((\mathds{D}{\bf
x}_{1})(0)-\dot{x}_{i})$
$\displaystyle+\lambda_{3}(x_{1}(N_{t})-x_{2}(N_{t}))+\lambda_{4}((\mathds{D}{\bf
x}_{1})(N_{t})-(\mathds{D}{\bf x}_{2})(N_{t})).$ (62)

Searching for the extremum of this functional numerically using the parameters
$\mu=0.5$, $\kappa=1$, $\xi=0.00071$ (c.f. ref. tsang_slimplectic_2015 ) with
initial conditions $x_{i}=1$ $\dot{x}_{i}=0$ on discretized paths with
$N_{t}=32$ steps and regularized SBP21 operator leads to the results for ${\bf
x}_{1}$ (red circles) and ${\bf x}_{2}$ (blue crosses) shown in fig. 12.

Figure 14: (Left panel) Numerical solution (red dots) of the discretized path
${\bf x}_{1}$ ($N_{t}=2048$, $\Delta t=0.1$) that optimizes the functional eq.
62 corresponding to the discretized $S_{\rm GIVP}$ using the SBP42 operator in
time. The path ${\bf x}_{2}$ takes on the same values. Continuum solution of
the Euler-Lagrange equation $\mu\ddot{x}(t)+\xi\dot{x}(t)+\kappa x(t)=0$ is
given as solid gray line. (Right panel) Comparison of the optimal solutions at
a different number of grids points $N_{t}=512$ (green crosses), $N_{t}=1024$
(blue crosses) and $N_{t}=2048$ (red circles) at late times.

In the left panel of fig. 13 we present the deviation between the numerically
determined critical point of eq. 58 and the true solution at the final time
$t_{2}=1$. The deviation of the derivative is given in the right plot. Blue
crosses denote the SBP21 operator case, while red circles refer to the SBP42
operator. The power law fits show that also for the damped harmonic oscillator
the convergence order agrees with the expectations from the conventional
formulation of IVPs svard_convergence_2021  and is one order higher in the
values of the solution than in the derivative of the solution. The SBP21
operator yields a $\Delta t^{2.03}$ improvement for the values of $x_{1}(1)$,
while the SBP42 operator exhibits $\Delta t^{3.04}$.

So far we have only investigated the short time behavior. However it is late-
time stability that plays the most important role for the utility of a
discretization scheme to the description of physical processes in practice.
For initial boundary value problems this calls for so-called error-bounded
schemes nordstrom_long_2018 ; kopriva_error_2017 ; nordstrom_error_2008 . In
the literature the damped harmonic oscillator is often used as a non-trivial
test-bed to evaluate the late-time stability and accuracy of numerical
solvers. Let us therefore determine the numerical solution to eq. 62 up to
$t_{2}=204.8$, which is shown in the left panel of fig. 14, based on the
regularized SBP42 operator at $\Delta t=0.1$. This choice of $t_{2}$ allows
the system to pass through multiple oscillations and to show a visible
reduction of the oscillation amplitude. As our SBP in time discretization is
inherently implicit, we find numerically that while the solution degrades in
accuracy as we increase the grid spacing $\Delta t\in\\{0.1,0.133,0.2,0.4\\}$
it remains bounded for all times. The behavior of the discrete solution for
different grid spacings $\Delta t$ at late times, is shown in the right panel
of fig. 14. One finds that the solution converges as the grid spacing is
decreased. The most pertinent error introduced by the discretization procedure
appears to be an artificial phase shift, the dispersion error, which however
vanishes as the continuum limit $\Delta t\to 0$ is approached.

A common quality criterion for numerical solvers in the physical sciences is
the reproduction of the system energy. We consider here as energy the
following Hamiltonian ${\cal H}=T+V=\frac{1}{2}\mu\dot{x}^{2}(t)+\kappa
x^{2}(t)$, which is plotted for different grid spacings in fig. 15. One can
clearly see that the discretization procedure leads to an overall shift in the
value of the energy and the appearance of oscillations around the mean value,
known as dispersion and diffusion errors (for a detailed exposition of the
dispersion errors of SBP operators see e.g. ref. linders2015uniformly  and
references therein). However, both the shift, as well as the oscillations
vanish with grid refinement and no artificial energy deviation or instability
is observed.

Figure 15: Numerical solution of the energy of the solution path ${\bf x}_{1}$
from functional eq. 62 using the regularized SBP42 operator. The energy is
plotted at a different number of grids points $N_{t}=512$ (green crosses),
$N_{t}=1024$ (blue crosses) and $N_{t}=2048$ (red circles) and compared to the
continuum solution from the Euler-Lagrange equation, which is shown as solid
black line. Note that refining the grid consistently leads to smaller
oscillations and the numerical solution approaches the continuum result.

##  5 Summary

We have presented a new and unified discretization strategy for a wide range
of initial value problems, formulated in terms of a generalized variational
problem. The classical trajectory is obtained as the critical point of an
action functional with doubled degrees of freedom, without the need to derive
equations of motion. Our approach is based on summation-by-parts finite
difference operators, which are regularized through a penalty term, associated
with the initial conditions using affine coordinates. We introduce the initial
conditions, as well as the identification of the forward and backward path in
the IVP in section 3.3 using Lagrange multipliers. The approach has been
successfully applied to systems with purely second order, first order and
mixed time derivatives in their equations of motion. Explicit scaling tests
for each of these systems confirmed that our prescription converges to the
true solution under grid refinement, with the same convergence rates, as
expected from the conventional formulation of IVPs. Convergence in the values
of the solution of the systems investigated here occurs with one higher order
in the grid spacing $\Delta t$ compared to the derivative of the solution. Our
approach is based on a single realization of the first order SBP finite
difference operator and does not require a separate treatment of systems with
equations of motion containing first or second order derivatives in time.

Our study describes a genuinely novel and versatile approach to the
computation of trajectories of classical systems without the need to resort to
their equation of motion. The extension to higher dimensions for the treatment
of partial differential equations is the natural next step and is work in
progress. We are looking forward to applying the formalism to the
discretization of theories with intrinsic constraints, in particular gauge
theories, such as classical electromagnetism and Yang-Mills theory.

## Acknowledgements

A. R. thanks Will Horowitz for stimulating discussions and gladly acknowledges
support by the Research Council of Norway under the FRIPRO Young Research
Talent grant 286883. J. N. was supported by the Swedish Research Council grant
nr. 2018-05084 and 2021-05484. The study has benefitted from computing
resources provided by UNINETT Sigma2 - the National Infrastructure for High
Performance Computing and Data Storage in Norway under project NN9578K-QCDrtX
”Real-time dynamics of nuclear matter under extreme conditions”

## Competing interests

The authors declare that they have no competing interests.

## Author’s contributions

  * •

A. Rothkopf: project inception, development of the discretization and
regularization prescription, scaling tests, writing and editing

  * •

J. Nordström: guidance on the SBP-SAT formalism, development of the
discretization formalism, literature review, writing and editing

##  Appendix A Noether Theorem for generalized variational problems

Following ref. sieberer2016keldysh  we retrace here how Noether’s theorem
connects global symmetries of the system with conserved quantities, the so
called Noether charges in the case of doubled degrees of freedom. Note that
transformations that act the same on the forward and backward path, due to the
construction of $L$, will lead to vanishing Noether currents. To identify a
finite conserved quantity, we need to consider transformations that act
differently on the paths $x_{1,2}$. The relevant symmetry in our case is the
invariance under opposite time translations. To formalize this statement, one
converts an infinitesimal time shift $\tau$ in the argument of the paths via a
Taylor expansion to ${\cal
T}_{\tau}x_{1,2}(t)=x_{1,2}(t)\pm\tau\dot{x}_{1,2}(t)+{\cal O}(\tau^{2})$. The
Lagrangian ${\cal L}$ is a scalar, just as the paths and therefore transforms
the same way as $\delta{\cal L}[x_{1,2}]/\delta\tau=\pm d{\cal
L}[x_{1,2}]/dt$. One thus obtains

$\displaystyle\delta L$ $\displaystyle=\Big{(}\frac{\partial{\cal L}}{\partial
x_{1}}\delta x_{1}+\frac{\partial{\cal
L}}{\partial\dot{x}_{1}}\delta\dot{x}_{1}-\frac{\partial{\cal L}}{\partial
x_{2}}\delta x_{2}-\frac{\partial{\cal
L}}{\partial\dot{x}_{2}}\delta\dot{x}_{2}\Big{)}$ (63)
$\displaystyle=\Big{(}\frac{\partial{\cal L}}{\partial x_{1}}\delta
x_{1}-\frac{d}{dt}\frac{\partial{\cal L}}{\partial\dot{x}_{1}}\delta
x_{1}-\frac{\partial{\cal L}}{\partial x_{2}}\delta
x_{2}+\frac{d}{dt}\frac{\partial{\cal L}}{\partial\dot{x}_{2}}\delta
x_{2}+\frac{d}{dt}\Big{\\{}\frac{\partial{\cal L}}{\partial\dot{x}_{1}}\delta
x_{1}-\frac{\partial{\cal L}}{\partial\dot{x}_{2}}\delta
x_{2}\Big{\\}}\Big{)}$

where we have used integration by parts to arrive at the second line. The
first four terms are nothing but the Euler Lagrange equations we obtained in
eq. 9, which vanish identically, so that we are left with

$\displaystyle\frac{\delta
L}{\delta\tau}=\frac{d}{dt}\Big{(}\frac{\partial{\cal
L}}{\partial\dot{x}_{1}}\dot{x}_{1}+\frac{\partial{\cal
L}}{\partial\dot{x}_{2}}\dot{x}_{2}\Big{)}=\frac{d}{dt}\Big{(}{\cal
L}[x_{1},\dot{x}_{1}]+{\cal L}[x_{2},\dot{x}_{2}]\Big{)}$ (64)

The last equal sign arises from the fact that the Lagrangian itself transforms
as a scalar. We thus arrive at the final expression for the conserved quantity
as

$\displaystyle\frac{d}{dt}\Big{(}\pi_{1}\dot{x}_{1}+\pi_{2}\dot{x}_{2}-{\cal
L}[x_{1},\dot{x}_{1}]-{\cal
L}[x_{2},\dot{x}_{2}]\Big{)}=\frac{d}{dt}\Big{(}H_{+}+H_{-}\Big{)}=0$ (65)

This equation states that the total Hamiltonian, i.e. the total energy of the
forward and backward path degrees of freedom is conserved in time. Had we
instead started with time translations that acted the same on the forward and
backward contour ${\cal
T^{\prime}}_{\tau}x_{1,2}(t)=x_{1,2}(t)+\tau\dot{x}_{1,2}(t)+{\cal
O}(\tau^{2})$, the corresponding Noether charge would be the difference
between the energies on the forward and backward path, i.e.

$\displaystyle\frac{d}{dt}\Big{(}\pi_{1}\dot{x}_{1}-\pi_{2}\dot{x}_{2}-{\cal
L}[x_{1},\dot{x}_{1}]+{\cal
L}[x_{2},\dot{x}_{2}]\Big{)}=\frac{d}{dt}\Big{(}H_{+}-H_{-}\Big{)}=0$ (66)

telling us that not only is the total energy preserved but the energy
difference between the forward and backward path must remain the same over
time.

It is interesting to realize that the variational principle of eq. 14, which
has been derived here in a fully classical context, identically arises as the
classical limit of IVPs in quantum field theory. As discussed in detail in
ref. berges_quantum_2007 , formulating initial value problems in the language
of Feynman’s path integral necessitates the introduction of doubled degrees of
freedom. The combination of the forward and backward path are referred to as
the Schwinger-Keldysh time contour. It turns out that the relative path
$x_{-}$ is related to the quantum contributions and taking $x_{-}\to 0$ is
therefore intimately related to the classical limit. The classical limit of
taking $\hbar\to 0$ actually enforces $x_{-}=0$. Ref. berges_quantum_2007
shows that the variation of the joint action with respect to $x_{-}$ is the
relevant expression that describes how the fluctuating quantum paths collapse
onto the deterministic classical trajectory, which indeed emerges after taking
the limit $x_{-}\to 0$. In the context of Noether’s theorem, as discussed in
ref. sieberer2016keldysh , transformations that act equally on forward and
backward path are associated with quantum Noether currents, which do not have
a finite expectation value. On the other hand transformations that couple the
forward and backward path can be considered as quantum transformations, which
lead to classical Noether currents that in turn can have a finite expectation
value even in the classical limit, as we saw in eq. 65.

## References

  * (1) Blazek, J.: Computational Fluid Dynamics: Principles and Applications. Butterworth-Heinemann, ??? (2015)
  * (2) Taflove, A., Hagness, S.C., Piket-May, M.: Computational Electromagnetics: the Finite-difference Time-domain Method vol. 3. Elsevier Amsterdam, The Netherlands, ??? (2005)
  * (3) Breuer, H.-P., Petruccione, F., et al.: The Theory of Open Quantum Systems. Oxford University Press on Demand, ??? (2002)
  * (4) Gattringer, C., Lang, C.: Quantum Chromodynamics on the Lattice: an Introductory Presentation vol. 788. Springer, ??? (2009)
  * (5) Svärd, M., Nordström, J.: Review of summation-by-parts schemes for initial–boundary-value problems. Journal of Computational Physics 268, 17–38 (2014)
  * (6) Fernández, D.C.D.R., Hicken, J.E., Zingg, D.W.: Review of summation-by-parts operators with simultaneous approximation terms for the numerical solution of partial differential equations. Computers & Fluids 95, 171–196 (2014)
  * (7) Lundquist, T., Nordström, J.: The SBP-SAT technique for initial value problems. Journal of Computational Physics 270, 86–104 (2014)
  * (8) Nordström, J.: A roadmap to well posed and stable problems in computational physics. Journal of Scientific Computing 71(1), 365–385 (2017)
  * (9) Nordström, J., Lundquist, T.: Summation-by-parts in time. Journal of Computational Physics 251, 487–499 (2013)
  * (10) Nordström, J., Lundquist, T.: Summation-by-parts in time: the second derivative. SIAM Journal on Scientific Computing 38(3), 1561–1586 (2016)
  * (11) Svärd, M., Carpenter, M.H., Nordström, J.: A stable high-order finite difference scheme for the compressible Navier–Stokes equations, far-field boundary conditions. Journal of Computational Physics 225(1), 1020–1038 (2007). doi:10.1016/j.jcp.2007.01.023. Accessed 2022-05-19
  * (12) Nordström, J., Forsberg, K., Adamsson, C., Eliasson, P.: Finite volume methods, unstructured meshes and strict stability for hyperbolic problems. Applied Numerical Mathematics 45(4), 453–473 (2003). doi:10.1016/S0168-9274(02)00239-8. Accessed 2022-05-19
  * (13) Carpenter, M.H., Fisher, T.C., Nielsen, E.J., Frankel, S.H.: Entropy Stable Spectral Collocation Schemes for the Navier–Stokes Equations: Discontinuous Interfaces. SIAM Journal on Scientific Computing 36(5), 835–867 (2014). doi:10.1137/130932193. Publisher: Society for Industrial and Applied Mathematics. Accessed 2022-05-19
  * (14) Ranocha, H., Öffner, P., Sonar, T.: Summation-by-parts operators for correction procedure via reconstruction. Journal of Computational Physics 311, 299–328 (2016). doi:10.1016/j.jcp.2016.02.009. Accessed 2022-09-15
  * (15) Abgrall, R., Nordström, J., Öffner, P., Tokareva, S.: Analysis of the SBP-SAT Stabilization for Finite Element Methods Part I: Linear Problems. Journal of Scientific Computing 85(2), 43 (2020). doi:10.1007/s10915-020-01349-z. Accessed 2022-05-19
  * (16) Gassner, G.J.: A Skew-Symmetric Discontinuous Galerkin Spectral Element Discretization and Its Relation to SBP-SAT Finite Difference Methods. SIAM Journal on Scientific Computing 35(3), 1233–1253 (2013). doi:10.1137/120890144. Publisher: Society for Industrial and Applied Mathematics. Accessed 2022-05-19
  * (17) Hesthaven, J.S., Gottlieb, D.: A Stable Penalty Method for the Compressible Navier–Stokes Equations: I. Open Boundary Conditions. SIAM Journal on Scientific Computing 17(3), 579–612 (1996). doi:10.1137/S1064827594268488. Publisher: Society for Industrial and Applied Mathematics. Accessed 2022-05-19
  * (18) Carpenter, M.H., Gottlieb, D., Abarbanel, S.: Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodology and application to high-order compact schemes. Journal of Computational Physics 111(2), 220–236 (1994)
  * (19) Svärd, M., Nordström, J.: On the convergence rates of energy-stable finite-difference schemes. Journal of Computational Physics 397, 108819 (2019)
  * (20) Linders, V., Nordström, J., Frankel, S.H.: Properties of Runge-Kutta-Summation-By-Parts methods. Journal of Computational Physics 419, 109684 (2020). doi:10.1016/j.jcp.2020.109684. Accessed 2022-05-19
  * (21) Svärd, M., Nordström, J.: Convergence of energy stable finite-difference schemes with interfaces. Journal of Computational Physics 429, 110020 (2021). doi:10.1016/j.jcp.2020.110020. Accessed 2022-05-23
  * (22) Ranocha, H., Nordström, J.: A new class of a stable summation by parts time integration schemes with strong initial conditions. Journal of Scientific Computing 87(1), 1–25 (2021)
  * (23) Ipp, A., Müller, D.: Implicit schemes for real-time lattice gauge theory. Eur. Phys. J. C 78(11), 884 (2018). doi:10.1140/epjc/s10052-018-6323-x. 1804.01995
  * (24) Arnold, V.I.: Mathematical Methods of Classical Mechanics, 2nd edition edn. Springer, New York (1989)
  * (25) Gustafsson, B.: High Order Difference Methods for Time Dependent PDE. Springer Series in Computational Mathematics. Springer, ??? (2007). https://books.google.no/books?id=guaT5yPiLAgC
  * (26) Elcott, S., Schröder, P.: Building your own DEC at home. In: ACM SIGGRAPH 2005 Courses. SIGGRAPH ’05, p. 8. Association for Computing Machinery, New York, NY, USA (2005). doi:10.1145/1198555.1198667. https://doi.org/10.1145/1198555.1198667 Accessed 2022-08-28
  * (27) Stern, A., Tong, Y., Desbrun, M., Marsden, J.E.: Geometric Computational Electrodynamics with Variational Integrators and Discrete Differential Forms. vol. 73, pp. 437–475 (2015). arXiv:0707.4470 [physics]. http://arxiv.org/abs/0707.4470 Accessed 2022-08-28
  * (28) Goldstein, H., Poole, C., Safko, J.: Classical Mechanics, 3rd edition edn. Pearson, San Francisco Munich (2001)
  * (29) Galley, C.R.: Classical Mechanics of Nonconservative Systems. Physical Review Letters 110(17), 174301 (2013). doi:10.1103/PhysRevLett.110.174301. Publisher: American Physical Society. Accessed 2022-01-04
  * (30) Sieberer, L.M., Buchhold, M., Diehl, S.: Keldysh field theory for driven open quantum systems. Reports on Progress in Physics 79(9), 096001 (2016)
  * (31) Inc., W.R.: Mathematica, Version 13.0.0. Champaign, IL, 2021. https://www.wolfram.com/mathematica
  * (32) Ranocha, H., Ostaszewski, K., Heinisch, P.: Discrete Vector Calculus and Helmholtz Hodge Decomposition for Classical Finite Difference Summation by Parts Operators. Communications on Applied Mathematics and Computation 2(4), 581–611 (2020). doi:10.1007/s42967-019-00057-2. Accessed 2022-05-01
  * (33) Wilson, K.G.: Confinement of quarks. Physical Review D 10(8), 2445–2459 (1974). doi:10.1103/PhysRevD.10.2445. Publisher: American Physical Society. Accessed 2022-05-19
  * (34) Ruggiu, A.A., Nordström, J.: Eigenvalue Analysis for Summation-by-Parts Finite Difference Time Discretizations. SIAM Journal on Numerical Analysis 58(2), 907–928 (2020). doi:10.1137/19M1256294. Publisher: Society for Industrial and Applied Mathematics. Accessed 2022-05-19
  * (35) Rothkopf, A.: Mathematica 12 implementation of a novel IVP discretization from a variational problem. Zenodo (2022). doi:10.5281/zenodo.6605375. https://zenodo.org/record/6605375 Accessed 2022-06-02
  * (36) Berges, J., Gasenzer, T.: Quantum versus classical statistical dynamics of an ultracold Bose gas. Physical Review A 76(3), 033604 (2007). doi:10.1103/PhysRevA.76.033604. Publisher: American Physical Society. Accessed 2022-01-30
  * (37) Tsang, D., Galley, C.R., Stein, L.C., Turner, A.: “SLIMPLECTIC” INTEGRATORS: VARIATIONAL INTEGRATORS FOR GENERAL NONCONSERVATIVE SYSTEMS. The Astrophysical Journal 809(1), 9 (2015). doi:10.1088/2041-8205/809/1/L9. Publisher: American Astronomical Society. Accessed 2022-01-04
  * (38) Nordström, J., Frenander, H.: On Long Time Error Bounds for the Wave Equation on Second Order Form. Journal of Scientific Computing 76(3), 1327–1336 (2018). doi:10.1007/s10915-018-0667-0. Accessed 2022-05-26
  * (39) Kopriva, D.A., Nordström, J., Gassner, G.J.: Error Boundedness of Discontinuous Galerkin Spectral Element Approximations of Hyperbolic Problems. Journal of Scientific Computing 72(1), 314–330 (2017). doi:10.1007/s10915-017-0358-2. Accessed 2022-05-26
  * (40) Nordström, J.: Error Bounded Schemes for Time-dependent Hyperbolic Problems. SIAM Journal on Scientific Computing 30(1), 46–59 (2008). doi:10.1137/060654943. Publisher: Society for Industrial and Applied Mathematics. Accessed 2022-05-26
  * (41) Linders, V., Nordström, J.: Uniformly best wavenumber approximations by spatial central difference operators. Journal of Computational Physics 300, 695–709 (2015)