Pseudo-Hamiltonian neural networks for learning partial differential
  equations
Pseudo-Hamiltonian neural networks (PHNN) were recently introduced for
learning dynamical systems that can be modelled by ordinary differential
equations. In this paper, we extend the method to partial differential
equations. The resulting model is comprised of up to three neural networks,
modelling terms representing conservation, dissipation and external forces, and
discrete convolution operators that can either be learned or be given as input.
We demonstrate numerically the superior performance of PHNN compared to a
baseline model that models the full dynamics by a single neural network.
Moreover, since the PHNN model consists of three parts with different physical
interpretations, these can be studied separately to gain insight into the
system, and the learned model is applicable also if external forces are removed
or changed.
Sølve Eidnes, Kjetil Olsen Lye
# Pseudo-Hamiltonian neural networks for learning partial differential equations

###### Abstract

Pseudo-Hamiltonian neural networks (PHNN) were recently introduced for learning dynamical systems that can be modelled by ordinary differential equations. In this paper, we extend the method to partial differential equations. The resulting model is comprised of up to three neural networks, modelling terms representing conservation, dissipation and external forces, and discrete convolution operators that can either be learned or be given as input. We demonstrate numerically the superior performance of PHNN compared to a baseline model that models the full dynamics by a single neural network. Moreover, since the PHNN model consists of three parts with different physical interpretations, these can be studied separately to gain insight into the system, and the learned model is applicable also if external forces are removed or changed.

## 1 Introduction

The field called physics-informed machine learning combines the strengths of physics-based models and data-driven techniques to achieve a deeper understanding and improved predictive capabilities for complex physical systems [35, 57]. The rapidly growing interest in this interdisciplinary approach is largely motivated by the increasing capabilities of computers to store and process large quantities of data, along with the decreasing costs of sensors and computers that capture and handle data from physical systems. Machine learning for differential equations can broadly be divided into two categories: the forward problem, which involves predicting future states from an initial state, and the inverse problem, which entails learning a system or parts of it from data. A wealth of recent literature exists on machine learning for the forward problem in the context of partial differential equations (PDEs). The proposed methods include neural-network-based substitutes for numerical solvers [20, 21, 55, 50], but also methods that can aid the solution process, e.g. by optimizing the discretization to be used in a solver [2]. The focus of this paper is however on the inverse problem, and much of the foundation for our proposed model can be found in recent advances in learning neural network models for ordinary differential equations (ODEs). Specifically, we build on recent works on models that incorporate Hamiltonian mechanics and related structures that underlie the physical systems we seek to model.
Greydanus et al. introduced Hamiltonian neural networks (HNN) in [26], for learning finite-dimensional Hamiltonian systems from data. They assume that the data \(q\in\mathbb{R}^{n}\), \(p\in\mathbb{R}^{n}\) is obtained from a canonical Hamiltonian system

\[\begin{pmatrix}\dot{q}\\ \dot{p}\end{pmatrix}=\begin{pmatrix}0&I_{n}\\ -I_{n}&0\end{pmatrix}\begin{pmatrix}\frac{\partial H}{\partial q}\\ \frac{\partial H}{\partial p}\end{pmatrix},\]

and aim to learn a neural network model \(H_{\theta}\) of the Hamiltonian \(H:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}\). This approach has since been further explored and expanded in a number of directions, which includes considering control input [61], dissipative systems [60, 25], constrained systems [24, 10], port-Hamiltonian systems [16, 18, 19], and metriplectic system formulations [36, 29]. A similar approach considering a Lagrangian formulation instead of a Hamiltonian is presented in [13]. Modifications of the models that focus on improved training from sparse and noisy data have been proposed in [12, 33, 15, 11].

In [23], we proposed pseudo-Hamiltonian neural networks (PHNN), which are models for learning a formulation that allows for damping of the Hamiltonian and external forces acting on the system. In addition, we allow for the Hamiltonian system to be non-canonical, and thus consider the formulation

\[\dot{x}=(S(x)-R(x))\nabla H(x)+f(x,t),\qquad x\in\mathbb{R}^{d}, \tag{1}\]

where \(S(x)=-S(x)^{T}\), and \(y^{T}R(x)y\geq 0\) for all \(y\). That is, \(S(x)\in\mathbb{R}^{d\times d}\) can be any skew-symmetric matrix and \(R(x)\in\mathbb{R}^{d\times d}\) can be any positive semi-definite matrix. Since we put no restrictions on the external forces, a pseudo-Hamiltonian formulation can in principle be obtained for any first-order ODE, which in turn can be obtained from any arbitrary-order ODE by a variable transformation. The formulation makes it possible to learn models that can be separated into internal dynamics and the external forces, i.e. \((\hat{S}_{\theta}(x)-\hat{R}_{\theta}(x))\nabla\hat{H}_{\theta}(x)\) and \(\hat{f}_{\theta}(x,t)\). This requires some sense of uniqueness in this separation, so certain restrictions need be put on the model to consider systems less general than (1). A major advantage that comes with this feature of the PHNN approach is that it makes it possible to learn a model of the system as if under ideal conditions even if data is sampled from a system affected by disturbances, if one assumes that an undisturbed system is closed and thus given only by the internal dynamics.

The motivating idea behind the present paper is to extend the framework of [23] to PDEs. In principe, one could always treat the spatially discretized PDE as a system of ODEs and apply the PHNN models of [23] to that. However, that would be disregarding certain structures we know to be present in the discretized PDE and would lead to inefficient models. Thus, compared to the ODE case, we consider different neural network architectures. Moreover, we will in the PDE setting impose some restrictions on the form of the external forces, in that we will not allow for them to depend on spatial derivatives of the solution. On the other hand, we will consider a more general form of the internal dynamics, where dissipation can result from a separate term and not just damping of the Hamiltonian.

Although HNN and extensions of this has attracted considerable attention in recent years, there has been very few studies on extending the methodology to PDEs. To our knowledge, the only prior works that consider HNN for PDEs are those of Matsubara et al. in [44] and Jin et al. in [32]. The latter has included a numerical example on the nonlinear Schrodinger equation. The former reference considers both Hamiltonian PDEs, exemplified by the Korteweg-de Vries equation, and an extension to dissipative PDEs, demonstrated by the Cahn-Hilliard equation. This paper has been a major inspiration for our work, especially on the neural network architecture we use to model the integrals in our PDE formulation. By generalizing to a wider class of PDEs that can
have conservative and dissipative terms at once, and also allowing for external forces, we largely expand the utility of this learning approach.

The extension of PHNN to PDEs we propose here should naturally also be put in context with other recent advances in learning of PDEs from data. Long et al. introduced PDE-Net in [40], and together with [44] this may be the work that is most comparable to what we present here. Their model is similar to the baseline model we will compare PHNN to in this paper, albeit less general. Their approach has two components: learning neural network models for the nonlinear terms in the PDE and identifying convolution operators that correspond to appropriate finite difference approximations of the spatial derivative operators present. They do however make considerable simplifying assumptions in their numerical experiments, e.g. only considering linear terms and a forward Euler discretization in time. Other works that have received significant attention are those that have focused on identifying coefficients of terms of the PDE, both in the setting where one assumes that the terms are known and approximate them by neural network models [50] and in the setting where one also identifies the terms present from a search space of candidates using sparse regression [51, 53, 34]. There has also been considerable recent research on learning operators associated with the underlying PDEs, where two prominent methods are Fourier neural operators (FNO) [1, 38] and deep operator network (DeepONet) [41]. These operators can e.g. map from a source term to solution states or from an initial state to future states; in the latter case, learning the operator equates to solving the forward problem of the PDE. The review paper [5] summarizes the literature on operator learning and system identification of PDEs, as well as recent developments on learning order reductions.

As will be demonstrated theoretically and experimentally in this paper, assuming a pseudo-Hamiltonian structure when solving the inverse problem for PDEs has both qualitative and quantitative advantages. The latter is shown by numerical comparisons to a baseline model on four test cases. The main qualitative feature of PHNN is that it is composed of up to six trainable submodels, which after training each can be studied for an increased understanding of the system we are modelling. Moreover, we could train a system affected by external forces and remove these from the model after training, so that we have a model unaffected by these disturbances. The code for this paper is built on the code for PHNN for ODEs developed for [23], and we have updated the GitHub repository [https://github.com/SINTEF/pseudo-hamiltonian-neural-networks](https://github.com/SINTEF/pseudo-hamiltonian-neural-networks) and the Python package phlearn with this extension to the PDE case.

The rest of this paper is organised as follows. In the next section, we explore the theoretical foundations upon which our method is based. Then the pseudo-Hamiltonian formulation and the class of PDEs we will learn are presented and discussed in Section 3. Section 4 is the centerpiece of the paper, as it is here we present the PHNN method for PDEs. We then dedicate a substantial portion of the paper to presenting and evaluating numerical results for various PDEs, in Section 5. The penultimate section is devoted to analysis of the results and our model, and a discussion of open questions to address in future research. We summarize the main results and draw some conclusions in the last section.

## 2 Background: Derivatives, discretizations and neural networks

Before delving into the pseudo-Hamiltonian formulation and the model we propose based on this, we will review and discuss some requisites for making efficient neural network models of systems
governed by PDEs.

### Learning dynamical systems

Consider the PDE

\[u_{t}=g(u^{P},x,t),\quad u\in H^{p}(\Omega),x\in\Omega\subseteq\mathbb{R}^{d},t\in \mathbb{R}, \tag{2}\]

first-order in time, where \(u^{P}\) means \(u\) itself and its partial derivatives up to the arbitrary order \(p\) with respect to the spatial variables \(x_{1},\ldots,x_{d}\). We seek to train a model \(\hat{g}_{\theta}\) of \(g\) so that solving \(u_{t}=\hat{g}_{\theta}(u^{P},x,t)\) leads to accurate predictions of the future states of the system. The universal approximation theorem [31, 14] states that \(g\) can be approximated with an arbitrarily small error by a neural network. In practice, we have to assume an abundance of observations of \(u_{t}\) and \(u^{P}\) at \(t\) and \(x\) to actually find a precise neural network approximation of \(g\). This brings us straight to one of the fundamental challenges of machine learning of differential equations: in a typical real-world setting, we can not expect to have data on the derivatives, neither temporal nor spatial. Thus we will have to depend on approximations obtained from discrete data. In this paper we will use sub- and superscript to denote discrete solution points in space resp. time. That is, \(u^{j}(x)=u(x,t^{j})\) and \(u_{i}(t)=u(x_{i},t)\), and we will suppress the arguments when they are not necessary. Let us consider the issue of time-discretization first, an issue shared by ODEs and PDEs alike, and defer the second issue to the next subsection.

In several of the papers introducing the most prominent recent methods for learning finite-dimensional dynamical systems, e.g. the original HNN paper [26] and the first paper by Brunton et al. on system identification [5], the derivatives of the solution are assumed to be known or approximated by finite differences. Approximating the time-derivative by the forward finite difference operator is equivalent to training using the forward Euler integrator, which is also what is done in the PDE-Net papers [40, 39]. There has however been several recent papers proposing more efficient training methods that incorporate other numerical integration schemes, see e.g. [12, 33, 15]. We follow [23, 45] and set up the training is such a way that we can use any mono-implicit integrator; that is, any integrator that relies explicitly on the solution in the times it integrates from and to. For the majority of the experiments in this paper, we use the implicit midpoint method, which is second-order, symplectic and symmetric. That is, we train the model \(\hat{g}_{\theta}\) by identifying the parameters \(\theta\) that minimize the loss function

\[\mathcal{L}_{g_{\theta}}=\bigg{\|}\frac{u^{j+1}-u^{j}}{\Delta t}-\hat{g}_{ \theta}\Big{(}\frac{(u^{P})^{j}+(u^{P})^{j+1}}{2},x,\frac{t^{j}+t^{j+1}}{2} \Big{)}\bigg{\|}_{2}^{2},\]

given for one training point \(u^{j}\) and barring regularization for now. This yields a considerable improvement over the forward Euler method at next to no additional computational cost, since the model in both cases is evaluated at only one point at each iteration of training. The option to use other integrators, including symmetric methods of order four and six, is readily implemented in the phlearn package, and we do demonstrate the need and utility of a fourth-order integrator in Section 5.3. For a thorough study of integrators especially suited for training neural network models of dynamical systems, we refer the reader to [45, 46].

### Spatial derivatives and convolution operators

Moving from finite-dimensional systems to infinite-dimensional systems introduces the issue of how to approximate spatial derivatives by the neural network models. Thankfully, a proposed solution
to this issue can be found in recent literature, as several works have noted the connection between finite difference schemes for differential equations and the convolutional neural network models originally developed for image analysis [8, 17, 52, 2, 9].

Given a function \(u\) and a kernel or filter \(w\), a discrete convolution is defined by

\[(u*w)(x_{i})=\sum_{j=-r}^{s}w_{j}u(x_{i-j}),\quad r,s\geq 0. \tag{3}\]

Here \(*\) is called the convolution operator, and the kernel \(w\) is a tensor containing trainable weights: \(w=\left[w_{-r},w_{-r+1},\ldots,w_{0},\ldots,w_{s-1},w_{s}\right]\). If the function \(u\) is periodic, so that \(u_{0}=u_{M}\), we obtain a circular convolution, which can be expressed by a circulant matrix applied on the vector \(u=\left[u_{0},\ldots,u_{M-1}\right]^{T}\).

A convolutional layer in a neural network can be represented as

\[y_{k}(u_{i})=\phi\big{(}(u*w_{k})(x_{i})+b_{k}\big{)}=\phi\bigg{(}\sum_{j=-r}^{ s}w_{kj}u_{i-j}+b_{k}\bigg{)},\]

where \(y_{k}(u_{i})\) is the output of the \(k\)-th feature map at point \(u_{i}\), \(w_{kj}\) are the weights of the kernel \(w_{k}\), \(b_{k}\) is the bias term, and \(\phi(\cdot)\) is an activation function. The width of the layer is \(r+s+1\), and this is usually referred to as the size of the convolution kernel, or filter. For our purpose it makes sense to have either \(r=0\) or \(r=s\), and the latter is the standard when convolutional neural networks are used in image analysis. Training the convolutional layer of a neural network constitutes of optimizing the weights and biases, which we collectively denoted by \(\theta\) in the previous subsection.

Similarly, a finite difference approximation of the \(n\)-th order derivative of \(u\) at a point \(x_{i}\) can also be expressed as applying a discrete convolution:

\[\frac{\mathrm{d}^{n}u(x_{i})}{\mathrm{d}x^{n}}\approx\sum_{j=-r}^{s}a_{j}u(x_ {i-j}), \tag{4}\]

where the finite difference weights \(a_{j}\) depend on the spatial grid. If we assume the spatial points to be equidistributed and let \(h:=x_{i+1}-x_{i}\), we have e.g.

\[\frac{\mathrm{d}u(x_{i})}{\mathrm{d}x} =\frac{u(x_{i+1})-u(x_{i})}{h}+\mathcal{O}(h),\] \[\frac{\mathrm{d}u(x_{i})}{\mathrm{d}x} =\frac{u(x_{i+1})-u(x_{i-1})}{2h}+\mathcal{O}(h^{2}),\] \[\frac{\mathrm{d}^{2}u(x_{i})}{\mathrm{d}x^{2}} =\frac{u(x_{i+1})-u(x_{i})+u(x_{i-1})}{h^{2}}+\mathcal{O}(h^{2}),\] \[\frac{\mathrm{d}^{3}u(x_{i})}{\mathrm{d}x^{3}} =\frac{u(x_{i+2})-2u(x_{i+1})+2u(x_{i-1})-u(x_{i-2})}{2h^{3}}+ \mathcal{O}(h^{2}).\]

Hence, a kernel size of two, with e.g. \(r=0\) and \(s=1\), is sufficient to obtain a first-order approximation of the first derivative, while a kernel size of three is sufficient and necessary to obtain second order approximations of first and second derivatives. Further, kernel size five is needed to approximate the third derivative. As noted by [52], higher order derivatives can be approximated either by increasing the kernel size or applying multiple convolution operations. In our models, we have designed neural networks where only the first layer is convolutional, and thus the kernel size restricts the order of the derivative we can expect to learn, while it also restricts the order of the approximations of these derivatives.
### Variational derivative

Given the function \(H\) depending on \(u\), \(x\) and the first derivative \(u_{x}\), let \(\mathcal{H}\) be the integral of \(H\) over the spatial domain:

\[\mathcal{H}[u]=\int_{\Omega}H(x,u,u_{x})\,\mathrm{d}x. \tag{5}\]

The variational derivative, or functional derivative, \(\frac{\delta\mathcal{H}}{\delta u}[u]\) of \(\mathcal{H}\) is defined by the property

\[\left\langle\frac{\delta\mathcal{H}}{\delta u}[u],v\right\rangle_{L_{2}}= \frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{\epsilon=0}\mathcal{H}[u+ \epsilon v]\quad\forall v\,\in H^{p}(\Omega). \tag{6}\]

When \(\mathcal{H}\) as here only depends on first derivatives, the variational derivative can be calculated explicitly by the relation

\[\frac{\delta\mathcal{H}}{\delta u}[u]=\frac{\partial H}{\partial u}-\frac{ \mathrm{d}}{\mathrm{d}x}\frac{\partial H}{\partial u_{x}},\]

assuming enough regularity in \(H\).

## 3 Pseudo-Hamiltonian formulation of PDEs

In this paper we consider the class of PDEs that can be written on the form

\[u_{t}=S(u^{P},x)\frac{\delta\mathcal{H}}{\delta u}[u]-R(u^{P},x)\frac{\delta \mathcal{V}}{\delta u}[u]+f(u^{P},x,t), \tag{7}\]

where \(S(u^{P},x)\) and \(R(u^{P},x)\) are operators that are skew-symmetric resp. positive semi-definite with respect to the \(L^{2}\) inner product, \(\mathcal{H}\) and \(\mathcal{V}\) are integrals of the form (5) and \(f:\mathbb{R}\times\mathbb{R}^{d}\times\mathbb{R}\rightarrow\mathbb{R}\). The naming of this class is a challenge, since similar but not identical classes have been called by a myriad of names in the literature. Ignoring the term \(f\), the class could be referred to as metriplectic PDEs, where the name is a portmanteau of metric and symplectic [28, 4]. The formulation is also similar to an infinite-dimensional variant of the General Equation for Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) formalism from thermodynamics [27, 47], except for \(f\) and the fact that \(R(u^{P},x)\) is positive instead of negative semi-definite. To be consistent with our previous work [23] and to make clear the connection to the vast recent literature on Hamiltonian neural networks, we say that (7) is the class of pseudo-Hamiltonian PDEs. This marks a generalization of the definition used in [23], in addition to the extension to infinite-dimensional systems, in that we here allow for \(\mathcal{H}\) and \(\mathcal{V}\) to be two different integrals. In the case \(\mathcal{V}=0\) and \(f(u^{P},x,t)=0\), we have the class of integral-preserving PDEs, which encompasses all (non-canonical) Hamiltonian PDEs [37]. That is, given the appropriate boundary conditions, e.g. periodic, the PDE will preserve the integral \(\mathcal{H}\), usually labeled as a Hamiltonian, or integral of motion, of the system. This follows from the skew-symmetry of \(S\):

\[\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}t}=\left\langle\frac{\delta\mathcal{H} }{\delta u}[u],\frac{\partial u}{\partial t}\right\rangle_{L^{2}}=\left\langle \frac{\delta\mathcal{H}}{\delta u}[u],S(u^{P},x)\frac{\delta\mathcal{H}}{ \delta u}[u]\right\rangle_{L^{2}}=0.\]

Similarly, if \(\mathcal{H}=0\) and \(f(u^{P},x,t)=0\) but \(\mathcal{V}\geq 0\), the PDE (7) will dissipate the integral \(\mathcal{V}\), and \(\mathcal{V}\) may be called a Lyapunov function.
**Example 1**.: Consider the KdV-Burgers (or viscous KdV) equation [56]

\[u_{t}+\eta uu_{x}-\nu u_{xx}-\gamma^{2}u_{xxx}=0. \tag{8}\]

This can be written on the form (16) with \(A\) and \(R\) both being the identity operator \(I\), \(S=\frac{\partial}{\partial x}\) and \(f(u,x,t)=0\), and

\[\mathcal{H}=-\int_{\Omega}\left(\frac{\eta}{6}u^{3}+\frac{\gamma^{2}}{2}u_{x} ^{2}\right)\,dx \tag{9}\]

and

\[\mathcal{V}=\frac{\nu}{2}\int_{\Omega}u_{x}^{2}\,dx. \tag{10}\]

We see this connection by deriving the variational derivatives

\[\frac{\delta\mathcal{H}}{\delta u}[u]=-(\frac{\eta}{2}u^{2}-\gamma^{2}u_{xx}) \tag{11}\]

and

\[\frac{\delta\mathcal{V}}{\delta u}[u]=-\nu u_{xx}. \tag{12}\]

We have that (8) reduces to the inviscid Burgers' equation for \(\eta=-1\) and \(\nu=\gamma=0\), the viscous Burgers' equation for \(\eta=-1\), \(\nu\neq 0\) and \(\gamma=0\), and the KdV equation for \(\nu=0\), \(\eta\neq 0\) and \(\gamma\neq 0\).

### Spatial discretization

In this section and this section only we will use boldface notation for vectors, to distinguish continuous functions and parameters from their spatial discretizations. Assume that values of \(u\) are obtained at grid points \(\mathbf{x}=[x_{0},\ldots,x_{M}]^{T}\). Following [22], we interpret these as quadrature points with non-zero quadrature weights \(\kappa=[\kappa_{0},\ldots,\kappa_{M}]^{T}\), and approximate the \(L_{2}\) inner product by a weighted discrete inner product:

\[\langle u,v\rangle=\int_{\Omega}u(x)v(x)\,\mathrm{d}x\approx\sum_{i=0}^{M} \kappa_{i}u(x_{i})v(x_{i})=u^{T}\mathrm{diag}(\kappa)v=:\langle u,v\rangle_{ \kappa}.\]

Let \(\mathbf{p}\) denote the discretization parameters that consist of \(\mathbf{x}\) and the associated \(\kappa\). Then, assuming that there exists a consistent approximation \(\mathcal{H}_{\mathbf{p}}(\mathbf{u})\) to \(\mathcal{H}[u]\) that depends on \(u\) evaluated at \(\mathbf{x}\), we define the discretized variational derivative by the analogue to (6)

\[\left\langle\frac{\delta\mathcal{H}_{\mathbf{p}}}{\delta\mathbf{u}}(\mathbf{u }),\mathbf{v}\right\rangle_{\kappa}=\frac{\mathrm{d}}{\mathrm{d}\epsilon} \bigg{|}_{\epsilon=0}\mathcal{H}_{\mathbf{p}}(\mathbf{u}+\epsilon\mathbf{v}) \quad\forall\mathbf{v}\,\in\mathbb{R}^{M+1}.\]

Thus, as shown in [22], we have a relationship between the discretized variational derivative and the gradient:

\[\frac{\delta\mathcal{H}_{\mathbf{p}}}{\delta\mathbf{u}}(\mathbf{u})=\mathrm{ diag}(\kappa)^{-1}\nabla_{\mathbf{u}}\mathcal{H}_{\mathbf{p}}(\mathbf{u}).\]

Furthermore, we approximate \(S(u^{P},x)\) and \(R(u^{P},x)\) by matrices \(S_{d}(\mathbf{u})\) and \(R_{d}(\mathbf{u})\) that are skew-symmetric resp. positive semi-definite with respect to \(\langle\cdot,\cdot\rangle_{\kappa}\). Then a spatial discretization of (7) is given by

\[\mathbf{u}_{t}=S_{d}(\mathbf{u})\frac{\delta\mathcal{H}_{\mathbf{p}}}{\delta \mathbf{u}}(\mathbf{u})-R_{d}(\mathbf{u})\frac{\delta\mathcal{V}_{\mathbf{p}} }{\delta\mathbf{u}}(\mathbf{u})+\mathbf{f}(\mathbf{u},\mathbf{x},t),\]
which may equivalently be written as

\[\mathbf{u}_{t}=S_{\mathbf{p}}(\mathbf{u})\nabla_{\mathbf{u}}\mathcal{H}_{\mathbf{ p}}(\mathbf{u})-R_{\mathbf{p}}(\mathbf{u})\nabla_{\mathbf{u}}\mathcal{V}_{\mathbf{p}}( \mathbf{u})+\mathbf{f}(\mathbf{u},\mathbf{x},t), \tag{13}\]

where \(S_{\mathbf{p}}(\mathbf{u}):=S_{d}(\mathbf{u})\operatorname{diag}(\kappa)^{-1}\) and \(R_{\mathbf{p}}(\mathbf{u}):=R_{d}(\mathbf{u})\operatorname{diag}(\kappa)^{-1}\) are skew-symmetric resp. positive semi-definite by the standard definitions for matrices.

Thus, upon discretizing in space, we obtain a system of ODEs (13) that is on a form quite similar to the generalized pseudo-Hamiltonian formulation considered in [23]. In fact, if \(\mathcal{V}=\mathcal{H}\), we obtain the system

\[\mathbf{u}_{t}=(S_{\mathbf{p}}(\mathbf{u})-R_{\mathbf{p}}(\mathbf{u}))\nabla_ {\mathbf{u}}\mathcal{H}_{\mathbf{p}}(\mathbf{u})+\mathbf{f}(\mathbf{u}, \mathbf{x},t).\]

Still, we do not recommend applying the PHNN method of [23] on this directly without taking into consideration what we know about \(\mathcal{H}_{\mathbf{p}}\). Specifically, we want to exploit that it is a discrete approximation of the integral (5), and can thus be expected to be given by a sum of \(M\) terms that each depend in the same way on \(u_{i}\) and the neighbouring points \(u_{i-1}\) and \(u_{i+1}\). Hence, as discussed in the next section, we will employ convolutional neural networks with weight sharing across the spatial discretization points.

**Example 2**.: Consider again the KdV-Burgers equation (8), on the domain \(\Omega=[0,P]\) with periodic boundary conditions \(u(0,t)=u(P,t)\). We assume that the \(M+1\) grid points are equidistributed and define \(h:=x_{i+1}-x_{i}=P/M\). We approximate the integrals (9) and (10) by

\[\mathcal{H}_{\mathbf{p}}=-h\sum_{i=0}^{M-1}\left(\frac{\eta}{6}u_{i}^{3}+ \frac{\gamma}{2}\left(\delta_{f}u_{i}\right)^{2}\right) \tag{14}\]

and

\[\mathcal{V}_{\mathbf{p}}=\frac{\nu}{2}h\sum_{i=0}^{M-1}\kappa_{i} \left(\delta_{f}u_{i}\right)^{2}, \tag{15}\]

where the operator \(\delta_{f}\) denotes forward difference, i.e. \(\delta_{f}u_{i}=\frac{u_{i+1}-u_{i}}{h}\). Furthermore, we approximate \(\partial_{x}\) by the matrix corresponding to the central difference approximation \(\delta_{c}\) defined by \(\delta_{c}u_{i}=\frac{u_{i+1}-u_{i-1}}{2h}\), i.e.

\[S_{d}=\frac{1}{2h}\begin{pmatrix}0&1&0&\cdots&0&-1\\ -1&0&1&0&\cdots\\ 0&-1&0&1&0&\cdots\\ &&\ddots&\ddots&\ddots&\\ &\cdots&0&-1&0&1\\ 1&0&\cdots&0&-1&0\end{pmatrix}\in\mathbb{R}^{M\times M},\]

where the first and last rows are adjusted according to the periodic boundary conditions.

To obtain (13), we have that \(S_{\mathbf{p}}=\frac{1}{h}S_{d}\) and \(R_{\mathbf{p}}=\frac{1}{h}I\), with \(I\) being the identity matrix, and take the gradients of the approximated integrals to find

\[\nabla_{u}\mathcal{H}_{\mathbf{p}} =-h\left(\frac{\eta}{2}\mathbf{u}^{2}-\gamma^{2}\delta_{c}^{2} \mathbf{u}\right),\] \[\nabla_{u}\mathcal{V}_{\mathbf{p}} =-h\,\nu\,\delta_{c}^{2}\mathbf{u},\]
where \(\mathbf{u}^{2}\) and \(\mathbf{u}^{3}\) denotes the element-wise square and cube of \(\mathbf{u}\), and \(\delta_{c}^{2}:=\delta_{f}\delta_{b}\) denotes the second order difference operator approximating the second derivative by \(\delta_{c}^{2}u_{i}=\frac{1}{2h}(u_{i+1}-2u_{i}+u_{i-1})\). Observe that \(\frac{\delta\mathcal{H}_{\mathbf{p}}}{\delta\mathbf{u}}(\mathbf{u})=\frac{1}{ h}\nabla_{\mathbf{u}}\mathcal{H}_{\mathbf{p}}(\mathbf{u})\) and \(\frac{\delta\mathcal{V}_{\mathbf{p}}}{\delta\mathbf{u}}(\mathbf{u})=\frac{1}{ h}\nabla_{\mathbf{u}}\mathcal{V}_{\mathbf{p}}(\mathbf{u})\) are consistent discrete approximations of (11) and (12). Moreover, they are second-order approximations of these variational derivatives, even though (14) and (15) are only first-order approximations of the integrals (9) and (10).

### Restricting the class by imposing assumptions

Without imposing any further restrictions, the formulation (7) can be applied to any PDE that is first-order in time and will not be unique for any system. Any contribution from the two first terms on the left-hand side could also be expressed in \(f\), and even if we restrict this term, the operators \(S\) and \(R\) are generally not uniquely defined for a corresponding integral. In the remainder of this paper, we will consider the cases where the operators \(S\) and \(R\) are linear and independent of \(x\) and \(u\), we assume \(R\) to be symmetric, and we will not let \(f\) depend on derivatives of \(u\). Furthermore, we apply the symmetric positive semi-definite operator \(A\) to the equation and require that this commutes with \(R\) and \(S\). We redefine \(S:=AS\), \(R:=AR\) and \(f(u,x,t):=Af(u,x,t)\), and thus get

\[Au_{t}=S\frac{\delta\mathcal{H}}{\delta u}[u]-R\frac{\delta\mathcal{V}}{\delta u }[u]+f(u,x,t), \tag{16}\]

where the new \(S\) is still skew-symmetric and the new \(R\) is still symmetric and positive semi-definite.

In the following we will denote the identity operator by \(I\) and the zero operator by \(0\), so that \(Iv=v\) and \(0v=0\) for any \(v\in L_{2}\). We note that the zero operator is positive semi-definite, symmetric and skew-symmetric, while the identity operator is symmetric and positive semi-definite, but not skew-symmetric.

## 4 The PHNN model for PDEs

Since we assume that the operators \(A\), \(S\) and \(R\) are independent of \(x\) and \(u\), the discretization of these operators will necessarily result in circulant matrices, given that \(u\) is periodic. That is, they can be viewed as discrete convolution operators. We thus set \(\hat{A}_{\theta}^{[k_{1}]}\), \(\hat{S}_{\theta}^{[k_{2}]}\) and \(\hat{R}_{\theta}^{[k_{3}]}\) to be trainable convolution operators, where \(k_{1}\), \(k_{2}\) and \(k_{3}\) denote the kernel sizes, and we impose symmetry on \(\hat{A}_{\theta}^{[k_{1}]}\) and \(\hat{R}_{\theta}^{[k_{3}]}\) and skew-symmetry on \(\hat{S}_{\theta}^{[k_{2}]}\). Furthermore, we let \(\hat{\mathcal{H}}_{\theta}\) and \(\hat{\mathcal{V}}_{\theta}\) be two separate neural networks that take input vectors of length \(M\), the number of spatial discretization points, and outputs a scalar. The neural network \(\hat{f}_{\theta}\) can take input vectors representing both \(u\), \(x\) and \(t\) and outputs a vector of length \(M\).

The full pseudo-Hamiltonian neural network model for PDEs is then given by

\[\hat{g}_{\theta}(u,x,t)=(\hat{A}_{\theta}^{[k_{1}]})^{-1}\Big{(}\hat{S}_{ \theta}^{[k_{2}]}\nabla\hat{\mathcal{H}}_{\theta}(u)-\hat{R}_{\theta}^{[k_{3}] }\nabla\hat{\mathcal{V}}_{\theta}(u)+k_{4}\hat{f}_{\theta}(u,x,t)\Big{)}, \tag{17}\]

where we also have introduced \(k_{4}\), which should be \(1\) or \(0\) depending on whether or not we want to learn a force term. Given a set of \(N\) training points \(\{(u^{j_{n}},u^{j_{n}+1},t^{j_{n}})\}_{n=1}^{N}\) varying across time and different stochastic realizations of initial conditions, we let the loss function be defined as

\[\mathcal{L}_{g_{\theta}}(\{(u^{j_{n}},u^{j_{n}+1},t^{j_{n}})\}_{n=1}^{N})= \frac{1}{N}\sum_{n=1}^{N}\bigg{|}\frac{u^{j_{n}+1}-u^{j_{n}}}{\Delta t}-\hat{g} _{\theta}\Big{(}\frac{u^{j_{n}}+u^{j_{n}+1}}{2},x,\frac{t^{j_{n}}+t^{j_{n}+1}}{2 }\Big{)}\bigg{|}^{2}, \tag{18}\]
if the implicit midpoint integrator is used. For the experiments in Section 5.3 we use the fourth-order symmetric integrator introduced in [23], and the loss function is amended accordingly.

### Implementation

PHNN is comprised of up to six trainable models; the framework is very flexible and assumptions may be imposed so one or more of the parts do not have to be learned. Moreover, careful considerations should be made on how to best model the different parts. In the following, we explain how we have set up the models in our code.

#### 4.1.1 Modelling \(\mathcal{H}\) and \(\mathcal{V}\)

The networks \(\hat{\mathcal{H}}_{\theta}\) and \(\hat{\mathcal{V}}_{\theta}\) consist of one convolutional layer with kernel size two followed by linear layers corresponding to convolutional layers with kernel size one, and then in the last layer performs a summation of the \(M\) inputs to one scalar. Following each of the first two layers, we apply the \(\tanh\) activation function. To impose the periodic boundary conditions, we pad the input to the convolutional layer by adding \(u(P)=u(0)\) at the end of the array of the discretized \(u\). A similar technique was suggested in [44], although they use a kernel of size three on the first convolutional layer. We opt to have a smaller filter, since kernel size two is sufficient to learn the forward difference approximation of the first derivative in the integrals, which in turn is sufficient to obtain second order approximation of the resulting variational derivative; this is shown for the KdV-Burgers equation in Example 2. If we want to be able to learn derivatives of order two in the integral, we would need kernel size three, and to pad the input with one element on each side. If we want to learn derivatives of order three or four, or if we want to learn third- or fourth-order approximations of the derivatives, we would need the kernel size of the convolutional layer to be five, and to pad the input by two elements on each side. This adjustment can easily be made in our code. For the examples in this paper, we would not gain anything by increasing the kernel size, because we only have up to first derivatives in the integrals and because the training data is generated using second order spatial discretizations. On the other hand, having a kernel of size two simplifies the learning and may facilitate superior performance over a model that does not rely on a pseudo-Hamiltonian structure and have to approximate up to third derivatives by convolutional neural networks.

#### 4.1.2 Modelling \(A\), \(S\), \(R\) and \(f\)

In (17), \(k=(k_{1},k_{2},k_{3},k_{4})\) are hyperparameters that determine the expressiveness of the model. Setting \(k_{1}=k_{2}=k_{3}=M\) and \(k_{4}=1\) means that we can approximate the general system (16), while setting \(k_{1}=k_{3}=1\), \(k_{2}=3\) and \(k_{4}=0\) would be sufficient to learn a model for the discretized KdV-Burgers system (8). In fact, if we set \(k_{3}=3\), the skew-symmetric operator \(\hat{S}_{\theta}^{[k_{2}]}\) is uniquely defined up to a constant, since we require \(w_{0}=0\) and \(w_{1}=-w_{-1}\). Moreover, since this constant would only amount to a scaling between the operator and the discrete variational derivative it is applied to, \(\hat{S}_{\theta}^{[k_{2}]}\) does not have to be trained in this case; determining \(w_{1}\) would just lead to a scaling of the second-order approximation of the first derivative in space that could be compensated by a scaling of \(\mathcal{H}\). Similarly, if the kernel size of \(A\) or \(R\) is 3, we could set \(w_{0}=1\) and learn a single parameter \(w_{1}=w_{-1}\) for each of these when training the model. This corresponds to learning a linear combination of the identity and the second-order approximation of the second derivative in space.
We model \(f\) by the neural network \(\hat{f}_{\theta}\) that may take either of the variables \(u\), \(x\) and \(t\) as input. This has three linear layers, i.e. convolutional layers with kernel size one, with the \(\tanh\) activation function after each of the first two. If \(\hat{f}_{\theta}\) depends on \(x\), periodicity on the domain \([0,P]\) is imposed in a similar fashion as suggested in [58, 43, 42] for hard-constraining periodic boundary conditions in physics-informed neural networks and DeepONet. That is, we replace the input \(x\) by the first two Fourier basis functions, \(\sin\left(\frac{2\pi}{P}x\right)\) and \(\cos\left(\frac{2\pi}{P}x\right)\), which is sufficient for expressing any \(x\)-dependent periodic function.

In the numerical experiments of the next section, we do not consider systems where \(A\) and \(R\) are anything other than linear combinations of the identity and the spatial second derivative, or \(S\) is anything other than the first derivative in space. Thus we set \(k=[3,3,3,1]\) in our most general model, which is expressive enough to learn all the systems we consider. We also consider what we call _informed_ PHNN models where we assume to have prior knowledge of the operators, affecting \(k\), and also what variables \(f\) depends on.

#### 4.1.3 Leakage of constant

If \(R\) is the identity, or a linear combination of the identity with differential operators, the separation between the dissipation term and the external force term in (16) is at best unique up to a constant, which means that there may be leakage of a constant between the two last terms of the PHNN model. Hence, we must make some assumptions about these terms to separate them as desired. If we want the external force term to be small, we may use regularization and penalize large values of \(\|\hat{f}_{\theta}\|\) during training. The option to do this is implemented in the phlearn package. However, for the numerical experiments in the next section, we have instead opted to assume that the dissipative term should be zero for the zero solution, and thus correct the two terms in question after training is done so that it adheres to this without changing the full model. That is, if we have the model

\[\hat{g}_{\theta}^{\text{pre}}(u,x,t)=(\hat{A}_{\theta}^{[k_{1}]})^{-1}\Big{(} \hat{S}_{\theta}^{[k_{2}]}\nabla\hat{\mathcal{H}}_{\theta}(u)-\hat{R}_{\theta} ^{[k_{3}]}\nabla\hat{\mathcal{V}}_{\theta}^{\text{pre}}(u)+k_{4}\hat{f}_{ \theta}^{\text{pre}}(u,x,t)\Big{)}, \tag{19}\]

when the last training step is performed, we set

\[\nabla\hat{\mathcal{V}}_{\theta}(u) :=\nabla\hat{\mathcal{V}}_{\theta}^{\text{pre}}(u)-k_{4}\nabla \hat{\mathcal{V}}_{\theta}^{\text{pre}}(0),\] \[\hat{f}_{\theta}(u,x,t) :=\hat{f}_{\theta}^{\text{pre}}(u,x,t)-\hat{R}_{\theta}^{[k_{3}]} \nabla\hat{\mathcal{V}}_{\theta}^{\text{pre}}(0)\]

to get our final model (17), which is equivalent to (19). Then we may remove the dissipation or external forces from the model simply by setting \(k_{3}=0\) or \(k_{4}=0\). Note, however, that this correction may not work as expected if the zero solution is far outside the domain of the training data, since the neural network \(\hat{\mathcal{V}}_{\theta}^{\text{pre}}\) like most neural networks generally extrapolates poorly. In that case, regularization is to be preferred.

#### 4.1.4 Algorithms

We refer to Algorithm 1 and Algorithm 2 for the training of the PHNN and baseline models, respectively.

## 5 Numerical experiments

In this section, we test how PHNN models perform on a variety of problems with different properties. In the experiments, we compare PHNNs to a baseline that do not have the pseudo-Hamiltonian
**Data:** Observations \(D=\{(t_{1},\vec{x}^{1},\vec{u}^{1}),\ldots,(t_{N},\vec{x}^{N},\vec{u}^{N}\}\)

**Data:** Number of epochs \(K\)

**Data:** Batch size \(M_{b}\)

**Data:** Initial CNN \(\hat{H}_{\theta},\hat{V}_{\theta}\)

**Data:** Initial DNN \(\hat{f}_{\theta}\)

**Data:** Matrices \(\hat{A}_{\theta}^{[k_{1}]}\), \(\hat{S}_{\theta}^{[k_{2}]}\) and \(\hat{R}_{\theta}^{[k_{3}]}\)

**Data:**\(g_{\theta}\) defined in (17)

**Data:** Loss function \(\mathcal{L}_{g_{\theta}}\) defined in (18)

**Result:** Parameters \(\theta\) for \(g_{\theta}\)

**for** \(k\) _in_ \(1\ldots K\)**do**

**for** _batch in Batches_ **do**

\(B:=\{(u^{j_{m}},u^{j_{m}+1},t^{j_{m}})\}_{m=1}^{M_{b}}\leftarrow\) DrawRandomBatch\((D,M_{b})\);

Step using \(\mathcal{L}_{g_{\theta}}(B)\) and \(\nabla_{\theta}\mathcal{L}_{g_{\theta}}(B)\)

**end**

**end**

**Algorithm 1:** The training phase of the PHNN algorithm

**Data:** Observations \(D=\{(t_{1},\vec{x}^{1},\vec{u}^{1}),\ldots,(t_{N},\vec{x}^{N},\vec{u}^{N}\}\)

**Data:** Number of epochs \(K\)

**Data:** Batch size \(M_{b}\)

**Data:** Initial CNN \(g_{\theta}\)

**Data:** Loss function \(\mathcal{L}_{g_{\theta}}\) defined in (18)

**Result:** Parameters \(\theta\) for \(g_{\theta}\)

**for** _k in_\(1\ldots K\)**do**

**for** _batch in Batches_ **do**

\(B:=\{(u^{j_{m}},u^{j_{m}+1},t^{j_{m}})\}_{m=1}^{M_{b}}\leftarrow\) DrawRandomBatch\((D,M_{b})\);

Step using \(\mathcal{L}_{g_{\theta}}(B)\) and \(\nabla_{\theta}\mathcal{L}_{g_{\theta}}(B)\)

**end**

**end**

**Algorithm 2:** The training phase of the baseline algorithm
structure but is otherwise as similar as possible to the PHNNs and is trained in the same way. We test either two or three PHNNs for each problem, in addition to a baseline model. The models we test on all problems are:

* PHNN (general): A PHNN model with kernel sizes \(k=[3,3,3,1]\) and an \(\hat{f}_{\theta}\) that depends on \(u\), \(x\) and t;
* PHNN (informed): A PHNN model where the operators \(A\), \(S\) and \(R\) are known a priori and \(\hat{f}_{\theta}\) depends only on the variable(s) which \(f\) depend on;
* Baseline: A model consisting of one neural network that takes \(u\), \(x\) and \(t\) as input, where the output is of the same dimension as \(u\). The network consists of two parts: First a five layer deep neural network with a hidden dimension of 20 and the tanh activation function, then a convolutional layer with kernel size five and activation function tanh, followed by two additional layers with hidden dimension 100 and the tanh activation function. The first five layers are meant to approximate any non-linear function (say \(u\mapsto\frac{1}{2}u^{2}\) in the case of Burgers' equation), while the convolutional layer is supposed to represent a finite-difference approximation of the spatial derivatives. The baseline model needs a kernel size of five to approximate the third and fourth derivatives present in the first resp. last example we consider.

### The KdV-Burgers equation

Consider again the KdV-Burgers equation from examples 1 and 2, but with external forces:

\[u_{t}+\eta(\frac{1}{2}u^{2})_{x}-\nu u_{xx}-\gamma^{2}u_{xxx}=f(x,t). \tag{20}\]

We consider the domain \([0,P]\) and assume periodic solutions \(u(0,t)=u(P,t)\), with \(P=20\). We let \(\eta=6\), \(\nu=0.3\), \(\gamma=1\) and \(f(x,t)=\sin\left(\frac{2\pi x}{P}\right)\). We generate training data from initial conditions

\[u(x,0)=2\sum_{l=1}^{2}c_{l}^{2}\,\mathrm{sech}^{2}\left(c_{l}\Big{(}\big{(}x+ \frac{P}{2}-d_{l}P\big{)}\bmod P-\frac{P}{2}\Big{)}\right), \tag{21}\]

where \(c_{1},c_{2}\) and \(d_{1},d_{2}\) are randomly drawn from the uniform distributions \(\mathcal{U}(\frac{1}{2},2)\) and \(\mathcal{U}(0,1)\) respectively. That is, the initial states are two waves of height \(2c_{1}^{2}\) and \(2c_{2}^{2}\) centered at \(d_{1}P\) and \(d_{2}P\), with periodicity imposed.

We compare the general and informed PHNN models to the baseline model on (20) with

\[f(x,t)=\frac{3}{5}\sin\big{(}\frac{4\pi}{P}x-t\big{)}. \tag{22}\]

Ten models were trained of each type using training sets consisting of 410 states, obtained from integrating 10 random initial states of the form (21) and evaluating the solution at every time step \(\Delta t=0.05\) until \(t=2\). The PHNN models were trained for 20 000 epochs, while the baseline models were trained for 50 000 epochs. At every epoch, the models were validated by integrating them to time \(t=2\) starting at three random initial states and calculating the average mean squared error (MSE) from these. The model with the lowest validation score after the last epoch was saved as the final model.
The PHNN models are quite sensitive to variations in the initialization of the learnable parameters of the model, an issue we discuss further in Section 6.1. Hence we get a large average MSE from these models when applying them on 10 different initial states, as evident in Table 1. However, for certain initial conditions the PHNN models are consistently outperforming the baseline model; see figures 1 and 2. Moreover, the best PHNN models perform well on all the test cases; of the 30 models trained, 10 of each type, the seven models with the lowest average MSE on 10 test sets are all PHNN models. Thus it would be advisable to run several PHNN models with different initializations of the neural networks and disregard those models who behave vastly different from the others.

Figure 2 gives a demonstration of one of the main qualitative features of the PHNN models: we can remove the force and dissipation from the model and still get an accurate solution of the system without these. In this figure, we have also extracted the external force part from the baseline model by

\[\hat{f}_{\theta}(x,t):=\hat{g}_{\theta}(u,x,t)-\hat{g}_{\theta}(u,0,0). \tag{23}\]

This works here, since the external force is independent of the solution states and the integrals are zero when \(u=0\), but it would not be an option in general. Moreover, we note that there is no way to separate the conservation and dissipation terms of the baseline model. This can however be done with the PHNN models, so that we also have a model for the energy-preserving KdV equation.

\begin{table}
\begin{tabular}{l r r} \hline \hline Model type & mean & std \\ \hline PHNN (general) & 9.60e+00 & 1.42e+01 \\ PHNN (informed) & 3.37e+01 & 3.45e+01 \\ Baseline & 1.90e+00 & 1.62e+00 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Mean and standard deviation of the MSE at \(t=2\), for 10 models of each type trained on the KdV–Burgers equation with the external force (22).

Figure 1: Predictions of the forced KdV–Burgers system system (20) with force (22), obtained from the best of 10 models of each model type, as evaluated by the mean MSE at \(t=2\) on predictions from 10 random initial states.
When the external force is explicitly dependent on time, we are generally not able to learn a model that is accurate beyond the temporal domain of the training data. For autonomous systems, we can make due with less training data. Consider thus instead of (22) the external force

\[f(x)=\frac{3}{5}\sin\big{(}\frac{4\pi}{P}x\big{)}. \tag{24}\]

We train models for (20) with this \(f\), where now we have training sets consisting of \(60\) states, from solutions obtained at every time step \(\Delta t=0.1\) from \(t=0\) to \(t=0.5\). As can be seen in figures 3 and 4, the PHNN models perform better than the baseline model, and especially so with increasing time. The worst-performing baseline models become unstable before the final test time \(t=4\). We also see that the most general PHNN model struggles to correctly separate the external force from the viscosity term with this amount of training data. This is however not an issue when the model is informed that the force is purely dependent on the spatial variable.

Figure 2: Solutions of the various learned models of the forced KdV–Burgers system (20) with \(f\) given by (22) at time \(t=2\). The line and the shaded area, barely visible in these plots, are the mean resp. standard deviation of \(10\) models of each type. The dashed black line is the ground truth. _Upper row:_ The original system (20) that the models are trained on. _Second row:_ The learned force approximating \(f\) in (20). _Third row:_ Predictions with the force \(f\) removed from the models. _Lower row:_ Predictions with the external force and the dissipation term removed from the models.
### The forced BBM equation

The Benjamin-Bona-Mahony (BBM) equation was introduced as an improvement on the KdV equation for modelling waves on a shallow surface [48, 3]. We consider this equation with a time- and state-dependent source term:

\[u_{t}-u_{xxt}+u_{x}+uu_{x}=f(u,t), \tag{25}\]

which can be written on the form (16) with \(A=1-\frac{\partial^{2}}{\partial x^{2}}\), \(S=\frac{\partial}{\partial x}\) and \(R=0\). This requires

\[\mathcal{H}=\frac{1}{2}\int_{\Omega}\Big{(}u^{2}+\frac{1}{3}u^{3}\Big{)}\, \mathrm{d}x.\]

As for the KdV-Burgers equation, we train the model on a forced system starting with a two-soliton initial condition. In this case, the initial states are given by

\[u(x,0)=3\sum_{l=1}^{2}(c_{l}-1)\operatorname{sech}^{2}\bigg{(}\frac{1}{2}\sqrt {1-\frac{1}{c_{l}}}\Big{(}\big{(}x+\frac{P}{2}-d_{l}P\big{)}\bmod P-\frac{P}{2 }\Big{)}\bigg{)}, \tag{26}\]

i.e. two waves of amplitude \(3(c_{1}-1)\) and \(3(c_{2}-1)\) centered at \(d_{1}P\) and \(d_{2}P\), where \(c_{1},c_{2}\) and \(d_{1},d_{2}\) are randomly drawn from the uniform distributions \(\mathcal{U}(1,4)\) and \(\mathcal{U}(0,1)\) respectively, and with periodicity imposed on \(\Omega=[0,P]\). We set \(P=50\) in the numerical experiments. Furthermore, we let

\[f(u,t)=\frac{1}{10}\sin(t)u. \tag{27}\]

In addition to the three models described in the introduction of this section, we also test a model that is identical to the most general PHNN model except that it does not include a dissipation term. We do this because it is not clearly defined whether or how (27) should be separated into a term that is constantly dissipative and one that is not. The most general model does learn that the system has a non-zero dissipative term; however, this term added to the learned force term is close

Figure 3: Predictions of the forced KdV–Burgers system (20) with \(f\) given by (24) obtained from the best of 10 models of each model type, as evaluated by the mean MSE at \(t=0.5\) on predictions from 10 random initial states.
to the ground truth force term. This is due to a leakage of a term \(\alpha u\) for some random constant \(\alpha\) between the terms, similar to the constant leakage described in Section 4.1.3, so that we learn an approximated integral \(\hat{\mathcal{V}}_{\theta}=\alpha\frac{\Delta x}{2}\sum_{i=0}^{M}u_{i}^{2}\) with \(\hat{R}_{\theta}^{[3]}=I\) and a corresponding external force \(\hat{f}_{\theta}=(\frac{1}{20}\sin{(t)}+\alpha)u\). This leakage could be combatted by regularization, i.e. by penalizing the mean absolute value of the dissipation term. We have not done that in the numerical experiments, but instead opted to also learn a model without the dissipation term and compare to this.

For every type of model, we trained 10 distinct models with random initializations for a total of 20 000 epochs for the PHNN models and 50 000 epochs for the baseline model. We used training data comprising 260 states, obtained from integrating 10 randomly drawn initial states with time step \(\Delta t=0.4\) from time \(t=0\) to time \(t=10\). A validation score at each epoch was generated by integrating the models to time \(t=1\) starting at three initial states and calculating the mean MSE, and the model with the lowest validation score was kept. After training was done, we integrated the models starting from 10 new arbitrary initial conditions to determine the average MSE at \(t=10\)

Figure 4: Solutions of the various learned models of the forced KdV–Burgers system with \(f\) given by (24). The line and the shaded area is the mean resp. standard deviation of predictions at \(t=4\) of 10 models of each type. The dashed black line is the ground truth. _Upper row:_ The original system. _Second row:_ The learned force approximating \(f\) in (24). _Third row:_ Predictions with the external force \(f\) removed from the models. _Lower row:_ Predictions with the force and the dissipation term removed.
By the error of the models as reported in Table 2, we see that all PHNN models perform better than the baseline model. Interestingly, the average MSE is lowest for the PHNN model where \(A\) and \(S\) has to be learned but \(R\) is known to be zero. However, the best PHNN model of all those that we trained is one of those informed of these operators, as seen in Figure 5. Note that the baseline model cannot be expected to learn a perfect model for this example, since the discrete approximation of \(A^{-1}S=(1-\frac{\partial^{2}}{\partial x^{2}})^{-1}\frac{\partial}{\partial x}\) used when generating the training data is a discrete convolution operator with kernel size bigger than five. The PHNN models are more stable and behave especially well with increasing time compared to the baseline model. Beyond \(t=10\), the accuracy of all models quickly deteriorates. We attribute this to the poor extrapolation abilities of neural networks; the models are not able to learn how the time-dependent \(f\) behaves beyond the temporal domain \(t\in[0,10]\) of the training data. Figure 6 shows the external forces learned by the PHNNs, and how well the models predict the system with these forces removed. For the general PHNN model, we have in this case also removed the dissipation term.

### The Perona-Malik equation

In addition to modelling physical systems, PDEs can be used for image restoration and denoising. For instance, if the heat equation is applied to a greyscale digital image, where the state \(u\) gives the intensity of each pixel, it will smooth out the image with increasing time. The Perona-Malik equation for so-called anisotropic diffusion is designed to smooth out noise but not the edges of an image [49]. Several variations of the equation exist. We consider the one-dimensional case, with a space-dependent force term, given by

\[u_{t}+\left(\frac{u_{x}}{1+u_{x}^{2}}\right)_{x}=f(x). \tag{28}\]

This is a PDE of the type (16) with \(A=I\), \(S=0\) and \(R=I\), and

\[\mathcal{V}[u]=\frac{1}{2}\int_{\Omega}\ln(1+u_{x}^{2})\,\mathrm{d}x.\]

Note that the equation can be written on the form \(u_{t}=\frac{\partial}{\partial x}\phi[u]+f(x)\) for \(\phi[u]=-\frac{u_{x}}{1+u_{x}^{2}}\), but this \(\phi[u]\) is not the variational derivative of any integral. We consider (28) on the domain \([0,P]\) with \(P=6\), and set

\[f(x)=10\sin\left(\frac{4\pi}{P}x\right) \tag{29}\]

\begin{table}
\begin{tabular}{l r r} \hline \hline Model type & mean & std \\ \hline PHNN (general) & 9.23e-01 & 1.04e+00 \\ PHNN (no diss. term) & 1.43e-01 & 9.91e-02 \\ PHNN (informed) & 2.22e-01 & 5.33e-01 \\ Baseline & 2.72e+00 & 3.67e+00 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Mean and standard deviation of the MSE at \(t=10\), for 10 models of each type trained on the BBM equation with an external force.
for the following experiments. The initial conditions are given by

\[u(x,0)=a-\sum_{l=1}^{2}\Big{(}h_{l}\Big{(}\tanh\big{(}b(x-d_{l})\big{)}-\tanh \big{(}b(x-P+d_{l})\big{)}\Big{)}\Big{)}+c\sin^{2}\left(r\pi x\right)\sin\left(s \pi x\right) \tag{30}\]

where \(a\in\mathcal{U}(-5,5)\), \(b\in\mathcal{U}(20,40)\), \(c\in\mathcal{U}(0.05,0.15)\), \(d_{l}\in\mathcal{U}(0.3,3)\), \(h_{l}\in\mathcal{U}(0.5,1.5)\), \(r\in\mathcal{U}(0.5,3)\) and \(s\in\mathcal{U}(10,20)\).

We trained 10 models of each type for 10 000 epochs, on 20 pairs of data at time \(t=0\) and \(t=0.02\). This corresponds to the original noisy image and and an image where the noise is almost completely removed, as judged by visual inspection. Because the step between these states is quite large, a high-order integrator is required to get an accurate approximation of the time-derivative. Indeed, models trained with the second-order implicit midpoint method fails to remove the noise as fast or accurately as the ground truth (28). Thus we use instead the fourth-order symmetric Runge-Kutta method (SRK4) introduced in [23]. This requires roughly four times the computational cost per epoch as using the midpoint method, but gives a considerably improved performance, as demonstrated in Figure 7.

Table 3 and Figure 8 reports the result of applying the learned models on an original noisy state (30) with \(a=1\), \(b=30\), \(c=0.15\), \(d_{1}=1\), \(d_{2}=2\), \(h_{1}=h_{2}=1\), \(r=2\) and \(s=15\). Interestingly, the general PHNN model performs better than the informed one. Moreover, the PHNN models perform better when the kernel size of the first convolutional layer of \(\hat{\mathcal{V}}_{\theta}\) is three instead of two. This indicates that the model does not learn the Perona-Malik equation but rather a different PDE

Figure 5: Predictions of the forced BBM system (25) obtained from the best of 10 models of each model type, as evaluated by the mean MSE at \(t=10\) on predictions from 10 random initial states.
that denoises the image. This may be as expected when we only train on initial states and end states. An odd-numbered filter size is the norm when convolutional neural networks are used for imaging tasks, since this helps to maintain spatial symmetry, and the improved performance with the same number of filters.

Figure 6: Solutions of the learned BBM system obtained from the different models. The line and the shaded area is the mean resp. standard deviation of predictions at \(t=9\) of 10 models of each type. The dashed black line is the ground truth. _Upper row:_ The original system (25) that the models are trained on. _Middle row:_ The learned force approximating \(f\) in (25). _Lower row:_ Predictions with the force \(f\) removed from the models.

Figure 7: The result at different times from integrating the mean of five general PHNN models trained on the Perona–Malik system (28) using two different integration schemes in the training. _Upper row:_ The second-order midpoint method. _Lower row:_ The fourth-order symmetric method SRK4.
a kernel of size three in \(\hat{\mathcal{V}}_{\theta}\) can perhaps be related to this.

### The Cahn-Hilliard equation

The Cahn-Hilliard equation was originally developed for describing phase separation [7], but has applications also in image analysis, and specifically image inpainting [6, 54]. Machine learning of pattern-forming PDEs, which include the Cahn-Hilliard and Allen-Cahn equations, has been studied in [59]. Results on applying PHNN to the Allen-Cahn equation is included in our GitHub repository. However, here we only consider the Cahn-Hilliard equation, with an external force, given by

\[u_{t}-(\nu u+\alpha u^{3}+\mu u_{xx})_{xx}=f(u,x). \tag{31}\]

\begin{table}
\begin{tabular}{l r r} \hline \hline Model type & mean & std \\ \hline PHNN (general) & 5.36e-04 & 1.58e-04 \\ PHNN (informed) & 4.14e-03 & 1.03e-03 \\ Baseline & 3.09e-03 & 3.08e-04 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Mean and standard deviation of the MSE at end time \(t=0.02\), for 10 models of each type trained on the Perona–Malik equation with an external force, and evaluated on predictions from 10 random initial states.

Figure 8: Perona–Malik at time \(t=0.02\), models and exact. The line plot is the average of 10 models of each type, while the shaded region indicates the standard deviation. _Upper row:_ The original system (31) that the models are trained on. _Middle row:_ The learned force approximating (29). _Lower row:_ Predictions with the force \(f\) removed from the models.
This is a dissipative PDE if the external force is zero, and it can be written on the form (16) with \(A=I\), \(S=0\) and \(R=-\frac{\partial^{2}}{\partial x^{2}}\), and

\[\mathcal{V}[u]=\frac{1}{2}\int_{\Omega}\left(\nu u^{2}+\frac{1}{2}\alpha u^{4} -\mu u_{x}^{2}\right)\,dx.\]

In the experiments, we have used \(\nu=-1\), \(\alpha=1\) and \(\mu=-\frac{1}{1000}\), and

\[f(u,x)=\begin{cases}30u&\text{if }0.3<x<0.7,\\ 0&\text{otherwise.}\end{cases}\]

The initial conditions of the training data are

\[u(x,0)=\sum_{l=1}^{2}\left(a_{l}\sin\left(c_{l}\frac{2\pi}{P}x\right)+b_{l} \cos\left(d_{l}\frac{2\pi}{P}x\right)\right) \tag{32}\]

on the domain \([0,P]\) with \(P=1\), where \(a_{l}\), \(b_{l}\), \(c_{l}\) and \(d_{l}\) are random parameters from the uniform distributions \(\mathcal{U}(0,\frac{1}{5})\), \(\mathcal{U}(0,\frac{1}{20})\), \(\mathcal{U}(1,6)\) and \(\mathcal{U}(1,6)\), respectively.

In addition to the models described in the introduction of this section, we also train a "lean" model, with \(k=[1,0,3,1]\) but no prior knowledge of how \(R\) looks. For each model type, 10 randomly initialized models were trained for 20 000 (for the PHNN models) or 50 000 (for the baseline models) epochs on different randomly drawn data sets consisting of a total of 300 states, at times \(t=0\), \(t=0.004\) and \(t=0.008\). At each epoch, the model was evaluated by comparing to the ground truth solution of three states at \(t=0.008\), and the model with the lowest MSE on this validation set was kept. The resulting 10 models were then evaluated on 10 random initial conditions and the mean MSE in the last time step was calculated from this. The mean and standard deviation from the 10 models of each type are given in Table 4, and the prediction of the model of each type with the lowest mean MSE is shown in Figure 9. In Figure 10 we give the results of the average of all models, including the standard deviation. Here we also show the learned external force and the prediction when this is removed from the model. The initial state of the plots in figures 9 and 10 is (32) with \(a_{1}=0.1\), \(a_{2}=0.06\), \(b_{1}=0.01\), \(b_{2}=0.02\), \(c_{1}=2\), \(c_{2}=5\), \(d_{1}=1\), \(b_{2}=2\).

We see from Figure 9 that the most general PHNN model may model the system moderately well, but it is highly sensitive to variations in the training data and the initialization of the neural networks in the model; from Figure (10) and Table 4 we see that this model may produce unstable predictions. In any case, the PHNN models struggle to learn the external force of this problem accurately without knowing \(R\), which we see by comparing the predictions of _PHNN (lean)_ and _PHNN (informed)_ in Figure 10, where the difference between the models is that the former has to learn an approximation \(\hat{R}_{\theta}^{[2]}\) of \(R\) and is not informed that \(f\) is not explicitly time-dependent.

## 6 Analysis of the models and further work

Here we provide some preliminary analysis of the PHNN models, which lays the groundwork for further analysis to be performed in the future.
### Stability with respect to initial neural network

The training of neural networks is often observed to be quite sensitive to the intial guesses for the weights and biases of the network. Here we test this sensitivity for both the PHNN model and the baseline model on the KdV-Burgers experiment in Section 5.1. We keep the training data fixed and re-generate the initial weights for the neural networks and rerun the training procedure described in algorithms 1 and 2. In Figure 11 we plot the solution at the final time together with the standard deviation and the pointwise maximum and minimum values for both the baseline model and the PHNN approach, where the standard deviation and maximum/minimum is computed across an ensemble of different initial weights for the deep neural network. In Figure 12 we plot the \(L^{2}\) error at the final time step against the exact solution for varying number of epochs, where the shaded areas represent the maximum and minimum values of an ensemble of varying initial weights of the neural network.

### Spatial discretization and training data

We will strive to develop PHNN further to make the models discretization invariant. For now, we settle with noting that this is already a property of our model in certain cases; a sufficiently well-trained informed PHNN model will be discretization invariant if the involved integrals do not depend on derivatives. Of the examples considered in this paper, that applies to the BBM equation, the inviscid Burgers' equation, and the Cahn-Hilliard equation if \(\mu=0\) in (31). Figure 13 shows how the learned BBM system can be discretized and integrated on spatial grids different from where

\begin{table}
\begin{tabular}{l r r} \hline \hline Model type & mean & std \\ \hline PHNN (general) & 1.18e+00 & 9.49e-01 \\ PHNN (lean) & 2.51e-01 & 1.42e-01 \\ PHNN (informed) & 7.01e-05 & 5.69e-05 \\ Baseline & 2.85e-01 & 7.09e-02 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Mean and standard deviation of the MSE at \(t=0.02\), for 10 models of each type trained on the forced Cahn–Hilliard system (31) and evaluated on predictions from 10 random initial states.

Figure 9: Predictions of the forced Cahn–Hilliard system obtained from the best of 10 models of each model type, as evaluated by the mean MSE at \(t=0.02\) on predictions from 10 random initial states.
there was training data.

For the experiments in the Section 5, we generated training data using first and second order finite difference operators to approximate the spatial derivatives. Then we trained our models on the same spatial grid, thus making it possible to learn convolution operators of kernel size two or three that perfectly captured the operators in the data. In a real-world scenario, this is unrealistic, as we would have to deal with discretization of a continuous system in space, as well as in time. We chose to disregard this issue in the experiments, to give a clearer comparison between PHNN and the baseline model not clouded by the error from the spatial discretization that would affect both. However, we have also tested our models on data generated on a spatial grid of four times as many discretization points, so that the data is generated from a more accurate approximation of the differential operators than what is possible to capture by the convolution operators. As we see in Figure 14, the PHNN models do appear to work well even with the introduction of approximation error in the spatial discretization. A more thorough study of this issue is required to gain a good understanding of how to best handle the spatial discretization.

### Proof of convergence in the idealized case

In this section we show a simplified error estimate for learning the right hand side of an ODE. Consider thus a model ODE of the form

Figure 10: Mean and standard deviation of predictions at \(t=0.02\) obtained from integrating 10 models of each type, for the Cahn–Hilliard problem (31). The dashed black line is the ground truth. _Upper row:_ The original system (31) that the models are trained on. _Middle row:_ The learned force approximating \(f\) in (31). _Lower row:_ Predictions with the force \(f\) removed from the models.
\[\begin{cases}\dot{u}(t)&=g(u(t),t)\\ u(0)&=u_{0},\end{cases} \tag{33}\]

where \(u:[0,T)\rightarrow\mathbb{R}^{M}\) and \(g:\mathbb{R}^{M}\rightarrow\mathbb{R}^{M}\). Note that the spatially discretized equation (13) can be cast into this form. In a certain sense, both the baseline model and the PHNN model tries to identify \(g\) by minimizing the \(L^{p}\)-norm of the observations \(u(t^{j})\) and the predictions \(u_{\theta}(t^{j})\). On a high level, this gives us a sequence \(u_{\theta}\to u\) in \(L^{p}\), but as is well-known, this would not be enough to conclude anything about the convergence of \(g_{\theta}\to g\), since \(L^{p}\) convergence in general does not imply convergence of the the derivatives. However, by utilizing the fact that we have a certain control over the discretized temporal derivatives in the learning phase, we can show that the \(g_{\theta}\to g\) in

Figure 11: Stability comparison of the baseline model and the general PHNN model. We retrain the models 20 times and compute the pointwise mean (plotted as a solid line) together with the standard deviation (plotted as error bars) and the pointwise maximum and minimum value (plotted as the shaded area).
the same \(L^{p}\) norm, _provided the training loss is small enough_. The following theorem makes this precise.

**Theorem 1**.: _Let \(\Delta t>0\), and \(g,\tilde{g}:\mathbb{R}^{M}\to\mathbb{R}^{M}\). Assume that \(u:[0,T)\to\mathbb{R}^{M}\) solves (33) and that

Figure 12: Convergence of the solution with respect to the number of training epochs. Here, the shaded area represents the maximum and minimum errors obtained, and the solid middle line represents the mean value across different initial weights.

Figure 13: The solution at \(t=10\) obtained from models of the BBM equation (25) with \(f=0\), learned from data discretized on \(M+1=101\) equidistributed points on the domain \([0,50]\). The \(M\) above the plots indicates the number of equidistributed discretization points used in the integration. The PHNN model was trained for 10 000 epochs, while the baseline model was trained for 100 000 epochs, on 20 pairs of states at time \(t=0\) and \(t=0.4\), with initial states (26).
\(\tilde{u}^{1},\ldots,\tilde{u}^{N}\in\mathbb{R}^{M}\) obey1

Footnote 1: This is essentially saying they are obtained during training for the loss function.

\[\frac{\tilde{u}^{j+1}-u^{j}}{\Delta t}=\tilde{g}(u^{j})\qquad\text{ for }j=0,\ldots,N-1. \tag{34}\]

_Then,_

\[\left(\Delta t\sum_{j=1}^{N-1}\left(g(u(t^{j}),t^{j})-\tilde{g}(u^{ j},t^{j})\right)^{p}\right)^{1/p}\leq\frac{1}{\Delta t}\left(\sum_{j=1}^{N} \Delta t\left|u^{j}-\tilde{u}^{j}\right|^{p}\right)^{1/p}+C_{g}\Delta t.\]

Proof.: Define \(u^{0},\ldots,u^{N}\in\mathbb{R}^{M}\) as

\[u^{j}:=u(t^{n})\qquad j=1,\ldots,N. \tag{35}\]

By a Taylor expansion, we have

\[\frac{u^{j+1}-u^{j}}{\Delta t}=g(u^{n})+\left(\left[\frac{\partial g }{\partial u}(u(\xi),\xi)\right]g(u(\xi,\xi))+\frac{\partial g}{\partial t}( u(\xi),\xi)\right)\Delta t\qquad\xi\in[t^{j},t^{j+1}],j=0,\ldots,N-1.\]

Figure 14: Solution obtained from models of the KdV–Burgers equation (8), i.e. without a force term, learned from 410 training states, with 10 different initial conditions and points equidistributed in time between \(t=0\) and \(t=2\). _Upper row:_ Models trained on data generated on a spatial grid with \(M=100\), same as used for training. _Lower row:_ Models trained on data generated on a spatial grid with \(M=400\) and then projected down on a grid with \(M=100\).
Hence, we get

\[\left(\sum_{j=0}^{N-1}\Delta t\left|g(u(t^{j}),t^{j})-\tilde{g}(u^{j},t^{j}) \right|^{p}\right)^{1/p}\leq\left(\sum_{j=1}^{N}\Delta t\left|\frac{u^{j}-u^{j- 1}}{\Delta t}-\frac{\tilde{u}^{j}-u^{j-1}}{\Delta t}\right|^{p}\right)^{1/p}+C_ {g}\Delta t.\]

We furthermore have

\[\left(\sum_{j=1}^{N}\Delta t\left|\frac{u^{j}-u^{j-1}}{\Delta t}- \frac{\tilde{u}^{j}-u^{j-1}}{\Delta t}\right|^{p}\right)^{1/p} =\left(\sum_{j=1}^{N}\Delta t\left|\frac{u^{j}-\tilde{u}^{j}}{ \Delta t}-\frac{u^{j-1}-u^{j-1}}{\Delta t}\right|^{p}\right)^{1/p}\] \[=\left(\sum_{j=1}^{N}\Delta t\left|\frac{u^{j}-\tilde{u}^{j}}{ \Delta t}\right|^{p}\right)^{1/p}\] \[=\frac{1}{\Delta t}\left(\sum_{j=1}^{N}\Delta t\left|u^{j}- \tilde{u}^{j}\right|^{p}\right)^{1/p}.\qed\]

**Remark 1**.: We note that (34) means that the sequence \(\{\tilde{u}^{j}\}_{j}\) is the learning data obtained by either the baseline or PHNN algorithm using the forward Euler integrator.

**Remark 2**.: In the above theorem, \(\left(\sum_{i=1}^{N}\Delta t\left|u^{j}-\tilde{u}^{i}\right|^{p}\right)^{1/p}\) is proportional to the training loss. In other words, the error in the approximation of \(g\) we get is bounded by \(1/\Delta t\cdot(\text{Training loss})\). Hence, to achieve an accuracy \(\epsilon\) in \(g\), we need to train to a training loss of \(\epsilon\Delta t\).

## 7 Conclusions

One of the advantages of PHNN is that it facilitates incorporating prior knowledge and assumptions into the models. The advantage of this is evident from the experiments in Section 5; the informed PHNN model performs consistently very well. We envision that our models can be used in an iterative process where you start by the most general model and as you learn more about the system from this, priors can be imposed.

As discussed in Section 6.1, the PHNN models are highly sensitive to variations in the initialized parameters of the neural networks. By the numerical results, the general PHNN model in particular seems to be more sensitive to this than the baseline model. However, the best trained PHNN models outperform the baseline models across the board. In a practical setting, where the ground truth is missing, we could train a number of models with different initialization of the neural networks and disregard those that deviate greatly from the others.

The aim of this paper has been to introduce a new method, demonstrate some of its advantages and share the code for the interested reader to study and develop further. It is the intent of the authors to also continue the work on these models. For one, we will be doing further analysis to address the issues raised in the previous section and improve the training of the models under various conditions. Secondly, we would like to extend the code to also work on higher-dimensional PDEs, and consider more advanced problems. One of the most promising uses of the methodology may be on image denoising and inpainting, motivated by the results of sections 5.3 and 5.4. Lastly,
the pseudo-Hamiltonian formulation could be used with other machine learning models than neural networks. Building on [30], we will develop methods for identifying analytic terms for one or several or all of the parts of the pseudo-Hamiltonian model (17), and compare the performance to existing system identification methods like [51, 53, 34]. We are especially intrigued by the possibility to identify the integrals of (16), while the external forces might be best modelled by a neural network.

### Acknowledgement

This work was supported by the research project PRAI (Prediction of Riser-response by Artificial Intelligence) financed by the Research Council of Norway, project no. 308832, with Equinor, BP, Subsea7, Kongsberg Maritime and Aker Solutions. The authors thank Katarzyna Michalowska and Signe Riemer-Sorensen for helpful comments on the manuscript, Eivind Bohn for help with coding issues, and Benjamin Tapley for both.
