# Energy-Dissipative Evolutionary Deep Operator Neural Networks

Jiahao Zhang111These two authors contributed equally to this work. Shiheng
Zhang222These two authors contributed equally to this work. Jie Shen
<EMAIL_ADDRESS>Guang Lin<EMAIL_ADDRESS>Department of Mathematics,
Purdue University, 150 N. University Street, West Lafayette, IN 47907-2067,
USA  School of Mechanical Engineering, Purdue University, 585 Purdue Mall,
West Lafayette, IN 47907-2067, USA

###### Abstract

Energy-Dissipative Evolutionary Deep Operator Neural Network is an operator
learning neural network. It is designed to seek numerical solutions for a
class of partial differential equations instead of a single partial
differential equation, such as partial differential equations with different
parameters or different initial conditions. The network consists of two sub-
networks, the Branch net, and the Trunk net. For an objective operator
$\mathcal{G}$, the Branch net encodes different input functions $u$ at the
same number of sensors $y_{i},i=1,2,\cdots,m$, and the Trunk net evaluates the
output function at any location. By minimizing the error between the evaluated
output $q$ and the expected output $\mathcal{G}(u)(y)$, DeepONet generates a
good approximation of the operator $\mathcal{G}$. In order to preserve
essential physical properties of PDEs, such as the Energy Dissipation Law, we
adopt a scalar auxiliary variable approach to generate the minimization
problem. It introduces a modified energy and enables unconditional energy
dissipation law in the discrete level. By taking the parameter as a function
of the time $t$ variable, this network can predict the accurate solution at
any further time with feeding data only at the initial state. The data needed
can be generated by the initial conditions, which are readily available. In
order to validate the accuracy and efficiency of our neural networks, we
provide numerical simulations of several partial differential equations,
including heat equations, parametric heat equations, and Allen-Cahn equations.

###### keywords:

Operator Learning
Evolutionary Neural Networks
Energy Dissipative
Parametric equation
Scalar auxiliary variable
Deep learning

††journal: Journal of Computational Physics

##  1 Introduction

Operator learning is a popular and challenging problem with potential
applications across various disciplines. The opportunity to learn an operator
over a domain in Euclidean spaces[1] and Banach spaces[2] opens a new class of
problems in neural network design with generalized applicability. In
application to solve partial differential equations(PDEs), operator learning
has the potential to predict accurate solutions for the PDE by acquiring
extensive prior knowledge [3, 4, 5, 6, 7, 8, 9, 10, 11]. In a recent
paper[12], Lu, Jin, and Karniadakis proposed an operator learning method with
some deep operator networks, named as DeepONets. It is based on the universal
approximation theorem [13, 14, 15]. The goal of this neural network is to
learn an operator instead of a single function, which is usually the solution
of a PDE. For any operator $\mathcal{G}$ on a domain $\Omega$, we can define
$\mathcal{G}$ as a mapping from $\Omega^{*}\rightarrow\Omega^{*}$ with
$\mathcal{G}(u)(y)\in R$ for any $y\in\Omega$. $\mathcal{G}(u)(y)$ is the
expected output of the neural network, which is usually a real number. The
objective of the training is to obtain an approximation of $\mathcal{G}$,
where we need to represent operators and functions in a discrete form. In
practice, it is very common to represent a continuous function or operator by
the values evaluated at finite and enough locations
$\\{x_{1},x_{2},\cdots,x_{m}\\}$, which is called “sensors” in DeepONet. The
network takes $[u(x_{1}),u(x_{2}),\cdots,u(x_{m})]$ and $y$ as the input. The
loss function is the difference between the output $q$ and the expected output
$\mathcal{G}(u)(y)$. Generally, there are two kinds of DeepONet, Stacked
DeepONet, and Unstacked DeepONet. The Stacked DeepONet consists of $p$ branch
networks and one trunk network. The number of the Trunk networks of the
Unstacked DeepONet is the same as the DeepONet, but the Unstacked DeepONet
merges all the $p$ branch networks into a single one. An Unstacked DeepONet
combines two sub-networks, Branch net, and Trunk net. The Branch net encodes
the input function $u$ at some sensors,
$\\{x_{i}\in\Omega\,|\,i=1,\cdots,m\\}$. The output of the Branch net consists
of $p$ neurons, where each neuron can be seen as a scalar,
$b_{j}=b_{j}(u(x_{1}),u(x_{2}),\cdots,u(x_{m}))$, $j=1,2,\cdots,p$. The Trunk
net encodes some evaluation points $\\{y_{k}\in\Omega|k=1,\cdots,n\\}$, while
the output also consists of $p$ neurons and each neuron is a scalar
$g_{j}=g_{j}(y_{1},y_{2},\cdots,y_{n})$, $j=1,2,\cdots,p$. The evaluation
point $y_{i}$ can be arbitrary in order to obtain the loss function. The
number of neurons at the last layer of the Trunk net and the Branch net is the
same. Hence, the output of the DeepONet can be written as an inner product of
$(b_{1},b_{2},\cdots,b_{p})$ and $(g_{1},g_{2},\cdots,g_{p})$. In other words,
the relationship between the expected output and the evaluated output is
$\mathcal{G}(u)(y)\approx\sum_{j=1}^{p}b_{j}g_{j}$. The DeepONet is an
application of the Universal Approximation Theorem for Operator, which is
proposed by Chen $\&$ Chen [16]:

######  Theorem 1.1 (Universal Approximation Theorem for Operator)

Suppose that $\Omega_{1}$ is a compact set in $X$, $X$ is a Banach Space, $V$
is a compact set in $C(\Omega_{1})$, $\Omega_{2}$ is a compact set in
$\boldsymbol{R}^{d}$, $\sigma$ is a continuous non-polynomial function,
$\mathcal{G}$ is a nonlinear continuous operator, which maps $v$ into
$C(\Omega_{2})$, then for any $\epsilon>0$, there are positive integers
$M,N,m$, constants $c_{i}^{k},\zeta_{k},\xi_{ij}^{k}\in\boldsymbol{R}$, points
$\omega_{k}\in\boldsymbol{R}^{n},x_{j}\in K_{1},i=1,\cdots,M$,
$k=1,\cdots,N,j=1,\cdots,m$, such that

$\begin{gathered}\mid\mathcal{G}(u)(y)-\sum_{k=1}^{N}\sum_{i=1}^{M}c_{i}^{k}\sigma\left(\sum_{j=1}^{m}\xi_{ij}^{k}u\left(x_{j}\right)+\theta_{i}^{k}\right)\cdot\sigma\left(\omega_{k}\cdot
y+\zeta_{k}\right)\mid<\epsilon\end{gathered}$

holds for all $u\in V$ and $y\in\Omega_{2}$.

(a) Stacked EDE-DeepONet

(b) Unstacked EDE-DeepONet

Fig. 1: Energy-Dissipative Evolutionary Deep Operator Neural Network. The
yellow block represents input at sensors and the blue block represents
subnetworks. The green blocks represent the output of the subnetworks and also
the last layer of the EDE-DeepONet. The difference between the stacked and
unstacked EDE-DeepONet is the number of Branch nets. In the right minimization
problem, the energy term $r^{2}$ can be shown to be dissipative, i.e.
$(r^{n+1})^{2}\leq(r^{n})^{2}$, where
$\mathcal{J}(\gamma_{1},\gamma_{2})=\frac{1}{2}\left\|\sum_{k=1}^{p}\frac{\partial{g_{k}(W^{n}_{1})}}{\partial
W^{n}_{1}}\gamma_{1}\boldsymbol{b}_{k}(W^{n}_{2})+\sum_{k=1}^{p}g_{k}(W^{n}_{1})\frac{\partial{\boldsymbol{b}_{k}(W^{n}_{2})}}{\partial
W^{n}_{2}}\gamma_{2}-\frac{r^{n+1}}{\sqrt{E({\boldsymbol{u^{n}}})}}\mathcal{N}_{\boldsymbol{x}}({\boldsymbol{u^{n}}})\right\|_{2}^{2}$.

For any time-dependent PDE, the training data is the form of
$(u,y,\mathcal{G}(u)(y))$, where $u$ in the discrete form can be represented
as $[u(x_{1}),u(x_{2}),\cdots,u(x_{m})]$ in the neural network. In the
original paper, they used the classic FNN[17] as the baseline model. For
dynamic systems, various network architectures are used, including residual
networks[18], convolutional NNs(CNNs)[19, 20], recurrent NNs(RNNs)[21], neural
jump stochastic differential equations[22] and neural ordinary differential
equations[23]. The training performance is very promising. It predicts
accurate solutions of many nonlinear ODEs and PDEs, including the simple
dynamic system, gravity pendulum system, and diffusion-reaction system.
However, the training data need to be generated at each time step, so it is
very expensive to train the network. For a lot of initial value problems,
there is no any information of $u(x,t)$ except $t=0$. It is very natural to
raise a question: Can we learn an operator of a kind of time-dependent PDEs
with only initial conditions?
Inspired by the Evolutionary Deep Neural Network(EDNN)[24], it is more
convenient to learn an operator at a fixed time instead of an operator with
not only spatial variables but also a time variable. With the loss of
generality, we can take the time variable $t$ to be $0$ in initial value
problems. Once obtained the operator at the initial time, many traditional
numerical methods can be used to update the solution. More specifically,
assuming that the initial condition operator has been trained well, we can
consider the parameters of the Branch net and the Trunk net as a function with
respect to the time variable as shown in Figure 1. More specifically, for a
given initial value problem,

$\left\\{\begin{aligned} &\frac{\partial u}{\partial t}=s(u)\\\
&u(x,0)=f(x),\quad x\in\omega\end{aligned}\right.$ (1)

the objective is to approximate the operator
$\mathcal{G}:u\mapsto\mathcal{G}(u)$. The input is
$([u(x_{1}),u(x_{2}),\cdots,u(x_{m})],y,\mathcal{G}(u)(y))$, where
$\\{x_{1},x_{2},\cdots,x_{m}\\}$ are the sensors and $\mathcal{G}(u)(y)=f(y)$.
The training process at the initial step is the same as the DeepONet, so we
can use the same architecture to train the initial condition operator. The
output of the Branch net can be written as
${\boldsymbol{b}}=\boldsymbol{b}(u(x_{1},0),u(x_{2},0),\cdots,u(x_{m_{1}},0))=\boldsymbol{b^{{}^{\prime}}}(x_{1},x_{2},\cdots,x_{m_{1}};W_{1})$,
where $W_{1}$ are the parameters in the Branch net. The output of the Trunk
net can be written as ${\boldsymbol{g}}=\boldsymbol{g}(y;W_{2})$, where
$W_{2}$ are the parameters in the Trunk net. Once trained well, we will regard
the parameters as a function of $t$ and $W_{1}$, $W_{2}$ as the initial
conditions of $W_{1}(t)$ and $W_{2}(t)$. By the architecture of the Unstacked
DeepONet, we can write the solution at initial time $t_{0}=0$ as

$u(x,t_{0})\approx\sum_{j=1}^{p}b_{j}g_{j}=\boldsymbol{b}^{T}\boldsymbol{g}\text{
for any given initial condition }f(x)$ (2)

We do not need any more data to obtain the approximation of $u(x,t_{1})$.
$u(x,t_{1})$ should be consistent with $W_{1}(t_{1})$ and $W_{2}(t_{1})$. With
the idea of the numerical solver for PDEs, it is easy to obtain $W_{1}(t_{1})$
and $W_{2}(t_{1})$ if $\frac{\partial W_{1}}{\partial t}$ and $\frac{\partial
W_{2}}{\partial t}$ are known. The time derivative of the solution $u$ can be
written by a chain rule:

$\frac{\partial u}{\partial t}=\frac{\partial u}{\partial W}\frac{\partial
W}{\partial t}$ (3)

where $W$ consists of $W_{1}$ and $W_{2}$. $\frac{\partial W}{\partial t}$ can
be solved by a least square problem. Once we get $\frac{\partial W}{\partial
t}$, we can use any traditional time discretization schemes to get $W^{n+1}$
with $W^{n}$.
The choice of the traditional time discretization scheme is dependent on the
specific problem. The Euler or Runge–Kutta methods are commonly used in the
evolutionary network. We are going to introduce a method with unconditional
energy dissipation, which is the Energy-Dissipative Evolutionary Deep Operator
Neural Network(EDE-DeepONet). Many kinds of PDEs are derived from basic
physical laws, such as Netwon’s Law, Conservation Law and Energy Dissipation
Law. In many areas of science and engineering, particularly in the field of
materials science, gradient flows are commonly employed in mathematical
models[25, 26, 27, 28, 29, 30, 31, 32]. When approximating the solution of a
certain PDE, it is desirable to satisfy these laws. We consider a gradient
flow problem,

$\frac{\partial{u}}{\partial t}=-\frac{\delta E}{\delta u},$ (4)

where $E$ is a certain free energy functional. Since the general explicit
Euler method does not possess the unconditionally dissipative energy
dissipation law, we applied a scalar auxiliary variable(SAV) method[33] to
generate the required least square problem. It introduces a new modified
energy and the unconditionally dissipative modified energy dissipation law is
satisfied for each iterative step. SAV method has been applied to solve plenty
of PDEs with thermodynamically consistent property. It is robust, easy to
implement and accurate to predict the solution. Introducing this method to
neural network helps us explore how to combine neural network models and
physical laws.
The objectives of this article is:

  * 1.

Designing an operator learning neural network without data except the given
information.

  * 2.

Predicting solutions of parametric PDEs after a long time period.

  * 3.

Keeping energy dissipative property of a dynamic system.

Our main contributions are:

  * 1.

Constructing an evolutionary operator learning neural network to solve PDEs.

  * 2.

Solving a kind of PDEs with different parameters in a single neural network.

  * 3.

Introducing the modified energy in the neural network and applying SAV
algorithm to keep the unconditionally modified energy dissipation law.

  * 4.

Introducing an adaptive time stepping strategy and restart strategy in order
to speed the training process.

The organization of this paper is as follows: In Section 2, we introduce the
Evolutionary Deep Operator Neural Network for a given PDE problem. In Section
3, we consider the physics law behind the gradient flow problem and apply the
SAV method to obtain the energy dissipation law. We proposed a new
architecture for neural network, EDE-DeepONet. In Section 4, we presented two
adaptive time stepping strategies, where the second one is called restart in
some cases. In Section 5, we generally introduced the architecture of the EDE-
DeepONet. In Section 6, we implement our neural network to predict solutions
of heat equations, parametric heat equations, and Allen-Cahn equations to show
the numerical results.

##  2 Evolutionary Deep Operator Neural Network

Consider a general gradient flow problem,

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}+\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u})=0$ (5)
$\displaystyle\boldsymbol{u}(\boldsymbol{x},0)=\boldsymbol{f}(\boldsymbol{x})$

where $\boldsymbol{u}\in\boldsymbol{R}^{l}$,
$\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u})$ can be written as a variational
derivative of a free energy functional $E[u(\boldsymbol{x})]$ bounded from
below, $\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u})=\frac{\delta
E}{\delta\boldsymbol{u}}$. The first step is to approximate the initial
condition operator with DeepONet.

###  2.1 Operator learning

For an operator $\mathcal{G}$,
$\mathcal{G}:\boldsymbol{u}(\boldsymbol{x})\mapsto\boldsymbol{f}(\boldsymbol{x})$,
the data feed into the DeepONet is in the form
$(\boldsymbol{u},y,\mathcal{G}(\boldsymbol{u})(y))$. It is obtained by the
given initial conditions. The branch network takes
$[\boldsymbol{u}(\boldsymbol{x_{1}}),\boldsymbol{u}(\boldsymbol{x_{2}}),\cdots,\boldsymbol{u}(\boldsymbol{x_{m}})]^{T}$
as the input, which is the numerical representation of $\boldsymbol{u}$, and
$[\boldsymbol{b}_{1},\boldsymbol{b}_{2},\cdots,\boldsymbol{b}_{p}]^{T}\in\boldsymbol{R}^{p\times
l}$, where $\boldsymbol{b}_{k}\in\boldsymbol{R}^{l}\text{ for
}k=1,2,\cdots,p$, as outputs. The trunk network takes $\boldsymbol{y}$ as the
input and $[g_{1},g_{2},\cdots,g_{p}]\in\boldsymbol{R}^{p}$ as outputs. The
Unstacked DeepONet net uses FNN as the baseline model and concatenate the
function value at sensor locations and the evaluated point together, i.e.
$[\boldsymbol{u}(\boldsymbol{x_{1}}),\boldsymbol{u}(\boldsymbol{x_{2}}),\cdots,\boldsymbol{u}(\boldsymbol{x_{m}}),\boldsymbol{y}]^{T}$.
As the equation in the Universal Approximation Theorem for Operators, we can
take the product of $\boldsymbol{h}$ and $t$, then we obtain:

$\mathcal{G}(\boldsymbol{u})({\boldsymbol{x}})\approx\sum_{k=1}^{p}g_{k}\boldsymbol{b}_{k}$
(6)

The activation functions are applied to the trunk net in the last layer. There
is no bias in this network. However, according to the theorem 1, the
generalization error can be reduced by adding bias. We also give the form with
bias $\boldsymbol{b_{0}}$:

$\mathcal{G}(\boldsymbol{u})({\boldsymbol{x}})\approx\sum_{k=1}^{p}g_{k}\boldsymbol{b}_{k}+\boldsymbol{b_{0}}$
(7)

As mentioned before, we assumed the initial condition operator has been
trained very well. We are going to find the update rule of the parameters to
evolve the neural network.

###  2.2 The evolution of parameters in the neural network

Denoting the parameters in the branch network as $W_{1}$ and the parameters in
the trunk network as $W_{2}$, $W_{1}$ and $W_{2}$ can be regarded a function
of $t$ since they change in every time step. According to the derivative’s
chain rule, we have

$\frac{\partial{\boldsymbol{u}}}{\partial
t}=\frac{\partial{\boldsymbol{u}}}{\partial
W_{1}}\frac{\partial{W_{1}}}{\partial
t}+\frac{\partial{\boldsymbol{u}}}{\partial
W_{2}}\frac{\partial{W_{2}}}{\partial t}$ (8)

Since
$\boldsymbol{u}=\sum_{k=1}^{p}g_{k}\boldsymbol{b}_{k}=\sum_{k=1}^{p}g_{k}(W_{1}(t))\boldsymbol{b}_{k}(W_{2}(t))$,
then

$\displaystyle\frac{\partial{\boldsymbol{u}}}{\partial t}$
$\displaystyle=\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}(t))}}{\partial
W_{1}}\frac{\partial{W_{1}}}{\partial
t}\boldsymbol{b}_{k}(W_{2}(t))+\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial{\boldsymbol{b}_{k}(W_{2}(t))}}{\partial
W_{2}}\frac{\partial{W_{2}}}{\partial t}$ (9)

Our objective is to obtain $\frac{\partial{W_{1}}}{\partial t}$ and
$\frac{\partial{W_{2}}}{\partial t}$, the update rule for parameters. It is
equivalent to solve a minimization problem,

$\left[\frac{\partial{W_{1}}}{\partial t};\frac{\partial{W_{2}}}{\partial
t}\right]=\text{argmin}\mathcal{J}(\gamma_{1},\gamma_{2})$ (10)

where

$\mathcal{J}(\gamma_{1},\gamma_{2})=\frac{1}{2}\left\|\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}(t))}}{\partial
W_{1}}\gamma_{1}\boldsymbol{b}_{k}(W_{2}(t))+\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial{\boldsymbol{b}_{k}(W_{2}(t))}}{\partial
W_{2}}\gamma_{2}-\mathcal{N}_{\boldsymbol{x}}({\boldsymbol{u}})\right\|_{2}^{2}$
(11)

In this article, the inner product $(a,b)$ is defined in the integral sense,
$(a,b)=\int_{\Omega}a(\boldsymbol{x})b(\boldsymbol{x})\mathrm{~{}d}\boldsymbol{x}$
and the $L_{2}$ norm is defined as $\left\lVert
a\right\rVert^{2}_{2}=\int_{\Omega}|a(\boldsymbol{x})|^{2}\mathrm{~{}d}\boldsymbol{x}$.
The minimization problem can be transformed into a linear system by the first-
order optimal condition:

$\displaystyle\frac{\partial{\mathcal{J}}}{\partial\gamma_{1}}=\int_{\Omega}\left(\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}(t))}}{\partial
W_{1}}\boldsymbol{b}_{k}(W_{2}(t))\right)^{T}\left(\gamma_{1}\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}(t))}}{\partial
W_{1}}\boldsymbol{b}_{k}(W_{2}(t))+\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial{\boldsymbol{b}_{k}(W_{2}(t))}}{\partial
W_{2}}\gamma_{2}-\mathcal{N}_{\boldsymbol{x}}({\boldsymbol{u}})\right)\mathrm{d}\boldsymbol{x}=0$
(12)
$\displaystyle\frac{\partial{\mathcal{J}}}{\partial\gamma_{2}}=\int_{\Omega}\left(\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial{\boldsymbol{b}_{k}(W_{2}(t))}}{\partial
W_{2}}\right)^{T}\left(\gamma_{1}\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}(t))}}{\partial
W_{1}}\boldsymbol{b}_{k}(W_{2}(t))+\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial{\boldsymbol{b}_{k}(W_{2}(t))}}{\partial
W_{2}}\gamma_{2}-\mathcal{N}_{\boldsymbol{x}}({\boldsymbol{u}})\right)\mathrm{d}\boldsymbol{x}=0$
(13)

In this system, the gradient with respect to $W_{1}(t)$ and $W_{2}(t)$ can be
computed by automatic differentiation at each time step. By denoting

$\displaystyle(\mathbf{J_{1}})_{i{j_{1}}}=\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}(t))}}{\partial
W^{j_{1}}_{1}}\boldsymbol{b}^{i}_{k}(W_{2}(t))$ (14)
$\displaystyle(\mathbf{J_{2}})_{i{j_{2}}}=\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial{\boldsymbol{b}^{i}_{k}(W_{2}(t))}}{\partial
W^{j_{2}}_{2}}$ (15)
$\displaystyle(\mathbf{N})_{i}=\mathcal{N}\left(\boldsymbol{u}_{\boldsymbol{x}}^{i}\right)$
(16)

where $i=1,2,\cdots,l$, $j_{1}=1,2,\cdots,N^{b}_{\text{para}}$,
$j_{2}=1,2,\cdots,N^{t}_{\text{para}}$. $N^{b}_{\text{para}}$ is the number of
parameters in Branch net and $N^{t}_{\text{para}}$ is the number of parameters
in Trunk net. $\mathbf{N}$ is generated by the DeepONet, so it can be
evaluated at any spatial point. The above integrals can be approximated by
numerical methods:

$\displaystyle\frac{1}{|\Omega|}\int_{\Omega}\left(\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}(t))}}{\partial
W_{1}}\boldsymbol{b}_{k}(W_{2}(t))\right)^{T}\left(\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}(t))}}{\partial
W_{1}}\boldsymbol{b}_{k}(W_{2}(t))\right)\mathrm{d}\boldsymbol{x}=\lim_{l\rightarrow\infty}\frac{1}{l}\mathbf{J^{T}_{1}}\mathbf{J_{1}}$
(17)
$\displaystyle\frac{1}{|\Omega|}\int_{\Omega}\left(\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial{\boldsymbol{b}_{k}(W_{2}(t))}}{\partial
W_{2}}\right)^{T}\left(\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial{\boldsymbol{b}_{k}(W_{2}(t))}}{\partial
W_{2}}\right)\mathrm{d}\boldsymbol{x}=\lim_{l\rightarrow\infty}\frac{1}{l}\mathbf{J^{T}_{2}}\mathbf{J_{2}}$
(18)
$\displaystyle\frac{1}{|\Omega|}\int_{\Omega}\left(\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}(t))}}{\partial
W_{1}}\boldsymbol{b}_{k}(W_{2}(t))\right)^{T}\left(\mathcal{N}_{\boldsymbol{x}}({\boldsymbol{u}})\right)\mathrm{d}\boldsymbol{x}=\lim_{l\rightarrow\infty}\frac{1}{l}\mathbf{J^{T}_{1}}\mathbf{N}$
(19)

By denoting $\gamma_{i}^{opt}$ as optimal values of $\gamma_{i}$, $i=1,2$, the
objective function can be reduced to

$\displaystyle\mathbf{J^{T}_{1}}\left(\gamma_{1}^{opt}\mathbf{J_{1}}+\gamma_{2}^{opt}\mathbf{J_{2}}-\mathbf{N}\right)=0$
(20)
$\displaystyle\mathbf{J^{T}_{2}}\left(\gamma_{1}^{opt}\mathbf{J_{1}}+\gamma_{2}^{opt}\mathbf{J_{2}}-\mathbf{N}\right)=0$
(21)

The feasible solutions of the above equations are the approximated time
derivatives of $W_{1}$ and $W_{2}$.

$\displaystyle\frac{dW_{1}}{dt}=\gamma_{1}^{opt}$ (22)
$\displaystyle\frac{dW_{2}}{dt}=\gamma_{2}^{opt}$ (23)

where the initial conditions $W_{1}^{0}$ and $W_{2}^{0}$ can be determined by
DeepONets for initial condition operators. The two ODEs are the updated rules
in the neural networks. The simple way to solve them is the explicit Euler
method.

$\displaystyle\frac{W^{n+1}_{1}-W^{n}_{1}}{\Delta t}=\gamma_{1}^{opt}$ (24)
$\displaystyle\frac{W^{n+1}_{2}-W^{n}_{2}}{\Delta t}=\gamma_{2}^{opt}$ (25)

The neural network can calculate the solution of given PDEs at any time step
$t_{n}$ and spatial point $\boldsymbol{x}_{i}$ by weights $W_{1}^{n}$ ,
$W_{2}^{n}$, spatial points $\boldsymbol{x}$ and initial condition
${\boldsymbol{u}}(\boldsymbol{x})$.

##  3 Energy Dissipative Evolutionary Deep Operator Neural Network

Let’s reconsider the given problem.

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}+\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u})=0$ (26)
$\displaystyle\boldsymbol{u}(\boldsymbol{x},0)=\boldsymbol{f}(\boldsymbol{x})$

where $\boldsymbol{u}\in\boldsymbol{R}^{l}$,
$\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u})$ can be written as a variational
derivative of a free energy functional $E[\boldsymbol{u}(\boldsymbol{x})]$
bounded from below, $\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u})=\frac{\delta
E}{\delta\boldsymbol{u}}$. Taking the inner product with
$\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u})$ of the first equation, we
obtain the energy dissipation property

$\frac{dE[\boldsymbol{u}(\boldsymbol{x})]}{dt}=\left(\frac{\delta
E}{\delta\boldsymbol{u}},\frac{\partial\boldsymbol{u}}{\partial
t}\right)=\left(\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}),\frac{\partial\boldsymbol{u}}{\partial
t}\right)=-\left(\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}),\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u})\right)\leq
0$ (27)

However, it is usually hard for a numerical algorithm to be efficient as well
as energy dissipative. Recently, the SAV approach [33] was introduced to
construct numerical schemes which is energy dissipative (with a modified
energy), accurate, robust and easy to implement. More precisely, assuming
$E[\boldsymbol{u}(\boldsymbol{x})]>0$, it introduces a
$r(t)=\sqrt{E[\boldsymbol{u}(\boldsymbol{x},t)]}$, and expands the gradient
flow problem as

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}=-\frac{r}{\sqrt{E(\boldsymbol{u})}}\mathcal{N}_{\boldsymbol{x}}\left(\boldsymbol{u}\right)$
(28) $\displaystyle
r_{t}=\frac{1}{2\sqrt{{E}(\boldsymbol{u})}}\left(\mathcal{N}_{\boldsymbol{x}}\left(\boldsymbol{u}\right),\frac{\partial\boldsymbol{u}}{\partial
t}\right)$

With $r(0)=\sqrt{E[\boldsymbol{u}(\boldsymbol{x},t)]}$, the above system has a
solution $r(t)\equiv\sqrt{E[\boldsymbol{u}(\boldsymbol{x},t)]}$ and
$\boldsymbol{u}$ being the solution of the original problem.

###  3.1 First order scheme

By setting $\boldsymbol{u}^{n}=\sum_{k=1}^{p}g_{k}\boldsymbol{b}_{k}$, a first
order scheme can be constructed as

$\displaystyle\frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\Delta
t}=-\frac{r^{n+1}}{\sqrt{E(\boldsymbol{u}^{n})}}\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}^{n})$
(29) $\displaystyle\frac{r^{n+1}-r^{n}}{\Delta
t}=\frac{1}{2\sqrt{E(\boldsymbol{u}^{n})}}\int_{\Omega}\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}^{n})\frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\Delta
t}dx.$

This is a coupled system of equations for $(r^{n+1},\boldsymbol{u}^{n+1})$.
But it can be easily decoupled as follows. Plugging the first equation into
the second one, we obtain:

$\frac{r^{n+1}-r^{n}}{\Delta
t}=-\frac{r^{n+1}}{2{E(\boldsymbol{u}^{n})}}\left\lVert\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}^{n})\right\rVert^{2},$
(30)

which implies

$r^{n+1}=\left(1+\frac{\Delta
t}{2{E(\boldsymbol{u}^{n})}}\left\lVert\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}^{n})\right\rVert^{2}\right)^{-1}r^{n}$
(31)

######  Theorem 3.2 (Discrete Energy Dissipation Law)

With the modified energy define above, the scheme is unconditionally energy
stable, i.e.

$(r^{n+1})^{2}-(r^{n})^{2}\leq 0.$ (32)

###### Proof 3.1

Taking the inner product of the first equation with
$\frac{r^{n+1}}{\sqrt{E(\boldsymbol{u}^{n})}}\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}^{n})$
and the second equation with $2r^{n+1}$

$\displaystyle(r^{n+1})^{2}-(r^{n})^{2}$
$\displaystyle=2r^{n+1}(r^{n+1}-r^{n})-(r^{n+1}-r^{n})^{2}$ (33)
$\displaystyle=\frac{\Delta
tr^{n+1}}{\sqrt{E(\boldsymbol{u}^{n})}}\int_{\Omega}\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}^{n})\frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\Delta
t}dx-(r^{n+1}-r^{n})^{2}$
$\displaystyle=-\left(\frac{r^{n+1}}{\sqrt{E(\boldsymbol{u}^{n})}}\right)^{2}\int_{\Omega}\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}^{n})\mathcal{N}_{\boldsymbol{x}}(\boldsymbol{u}^{n})dx-(r^{n+1}-r^{n})^{2}$
$\displaystyle\leq 0$

In order to maintain the modified energy dissipation law in the evolution
neural network, we only need to replace
$\mathcal{N}_{\boldsymbol{x}}({\boldsymbol{u}})$ by
$\frac{r^{n+1}}{\sqrt{E({\boldsymbol{u}})}}\mathcal{N}_{\boldsymbol{x}}({\boldsymbol{u}})$
in section 2. The update rule of the neural network is

$\left[\frac{\partial{W_{1}}}{\partial t};\frac{\partial{W_{2}}}{\partial
t}\right]=\text{argmin}\mathcal{J}(\gamma_{1},\gamma_{2})$ (34)

where

$\mathcal{J}(\gamma_{1},\gamma_{2})=\frac{1}{2}\left\|\sum_{k=1}^{p}\frac{\partial{g_{k}(W^{n}_{1})}}{\partial
W^{n}_{1}}\gamma_{1}\boldsymbol{b}_{k}(W^{n}_{2})+\sum_{k=1}^{p}g_{k}(W^{n}_{1})\frac{\partial{\boldsymbol{b}_{k}(W^{n}_{2})}}{\partial
W^{n}_{2}}\gamma_{2}-\frac{r^{n+1}}{\sqrt{E({\boldsymbol{u^{n}}})}}\mathcal{N}_{\boldsymbol{x}}({\boldsymbol{u^{n}}})\right\|_{2}^{2}$
(35)

The corresponding linear system of the first order optimal condition is

$\displaystyle\mathbf{J^{T}_{1}}\left(\gamma_{1}^{opt}\mathbf{J_{1}}+\gamma_{2}^{opt}\mathbf{J_{2}}-\frac{r^{n+1}}{\sqrt{E({\boldsymbol{u^{n}}})}}\mathbf{N}\right)=0$
(36)
$\displaystyle\mathbf{J^{T}_{2}}\left(\gamma_{1}^{opt}\mathbf{J_{1}}+\gamma_{2}^{opt}\mathbf{J_{2}}-\frac{r^{n+1}}{\sqrt{E({\boldsymbol{u^{n}}})}}\mathbf{N}\right)=0$
(37)

where

$\displaystyle(\mathbf{J_{1}})_{i{j_{1}}}=\sum_{k=1}^{p}\frac{\partial{g_{k}(W_{1}^{n})}}{\partial
W^{n,j_{1}}_{1}}\boldsymbol{b}^{i}_{k}(W_{2}^{n})$ (38)
$\displaystyle(\mathbf{J_{2}})_{i{j_{2}}}=\sum_{k=1}^{p}g_{k}(W^{n}_{1})\frac{\partial{\boldsymbol{b}^{i}_{k}(W^{n}_{2})}}{\partial
W^{n,j_{2}}_{2}}$ (39)
$\displaystyle(\mathbf{N})_{i}=\mathcal{N}\left(\boldsymbol{u}_{\boldsymbol{x}}^{i}\right)$
(40)

and $i=1,2,\cdots,l$, $j_{1}=1,2,\cdots,N^{b}_{\text{para}}$,
$j_{2}=1,2,\cdots,N^{t}_{\text{para}}$. $N^{b}_{\text{para}}$ is the number of
parameters in Branch net and $N^{t}_{\text{para}}$ is the number of parameters
in Trunk net. After getting $\gamma_{1}^{opt}$ and $\gamma_{2}^{opt}$,
$W^{n+1}$ can be obtained by the Forward Euler method as equation (24) and
(25).

$\displaystyle W^{n+1}_{1}=W^{n}_{1}+\gamma_{1}^{opt}{\Delta t}$ (41)
$\displaystyle W^{n+1}_{2}=W^{n}_{2}+\gamma_{2}^{opt}{\Delta t}$ (42)

##  4 Adaptive time stepping strategy and Restart strategy

One of the advantages of an unconditionally stable scheme is that the adaptive
time step can be utilized. Since the coefficient of $N_{x}$,
$\frac{r^{n+1}}{\sqrt{E^{n}}}$ should be around 1, by denoting
$\xi^{n+1}=\frac{r^{n+1}}{\sqrt{E^{n}}}$, larger $\Delta t$ is allowed when
$\xi$ is close to $1$ and the smaller $\Delta t$ is needed when $\xi$ is far
away from $1$. Thus, a simple adaptive time-stepping strategy can be described
as follows:

Algorithm 1 Adaptive time stepping strategy

1\. Set the tolerance for $\xi$ as $\epsilon_{0}$ and $\epsilon_{1}$, the
initial time step $\Delta t$, the maximum time step $\Delta t_{max}$ and the
minimum time step $\Delta t_{min}$

2\. Compute $u^{n+1}$.

3\. Compute $\xi^{n+1}=\frac{r^{n+1}}{\sqrt{E^{n}}}$.

4\. If $|1-\xi^{n+1}|>\epsilon_{0}$,

Then $\Delta t=\max(\Delta t_{min},\Delta t/2)$;

Else if $|1-\xi^{n+1}|<\epsilon_{1}$,

Then $\Delta t=\min(\Delta t_{max},2\Delta t)$.

Go to Step 2.

5\. Update time step $\Delta t$.

Another popular strategy to keep $r$ approximating the original energy $E$ is
to reset the SAV $r^{n+1}$ to be $E^{n+1}$ in some scenarios. The specific
algorithm is as following:

Algorithm 2 Restart strategy

1\. Set the tolerance for $\xi$ as $\epsilon_{2}$.

2\. Compute $u^{n+1}$.

3\. Compute $\xi^{n+1}=\frac{r^{n+1}}{\sqrt{E^{n}}}$.

4\. If $|1-\xi^{n+1}|>\epsilon_{2}$,

Then $r^{n+1}=\sqrt{E^{n+1}}$ and Go to Step 2.

5\. Go to next iteration.

The choice for $\epsilon_{0}$, $\epsilon_{1}$ should be some small tolerance,
usually $10^{-1}$ and $10^{-3}$. The choices for $\Delta t_{max}$ and $\Delta
t_{min}$ are quite dependent on $\Delta t$, usually $\Delta
t_{max}=10^{3}\times\Delta t$ and $\Delta t_{min}=10^{-3}\times\Delta t$. In
Algorithm 2, we usually take $\epsilon_{2}$ as $2\times 10^{-2}$.

##  5 Algorithm for EDE-DeepONet

A general approach to solving a time-dependent PDE with EDE-DeepONet can be
summarized in Algorithm 3.

Algorithm 3 Energy-Dissipative Evolutionary Deep Operator Neural Networks(EDE-
DeepONet)

1\. Generate input data samples in the form of $(u,y,\mathcal{G}(u)(y))$ for
the DeepONet, where $\mathcal{G}$ is the objective operator. Each specific
input function $u$ can be generated in the same sensor locations
$\\{x_{1},x_{2},\cdots,x_{m}\\}$.

2\. Feed $[u(x_{1}),u(x_{2}),\cdots,u(x_{m})])$ into the branch network and
$y\in Y$ into the trunk network. Denote the output of the DeepONet as $q$.

3\. Update the parameters in the DeepONet by minimizing a cost function, where
the cost function can be taken as the mean squared error as
$\frac{1}{|Y|}\sum_{y\in Y}\left\lVert\mathcal{G}(u)(y)-q\right\rVert^{2}$.

4\. Once the DeepONet has been trained well, solve the system of equations of
(36) and (37) to obtain $\left[\frac{\partial{W_{1}}}{\partial
t};\frac{\partial{W_{2}}}{\partial t}\right]$.

5.The value of $\left[\frac{\partial{W_{1}}}{\partial
t};\frac{\partial{W_{2}}}{\partial t}\right]$ can be obtained in the current
step. Since the parameters $W_{1}^{n}$ in the branch network, and $W_{2}^{n}$
in the trunk network are known, $W_{1}^{n+1}$ and $W_{2}^{n+1}$ for the next
step can be also obtained by the Forward Euler method or Runge-Kutta method.

6\. Repeat step 5 until the final time $T$, where $T=t_{0}+s\Delta t$, $t_{0}$
is the initial time of the given PDE, $\Delta t$ is the time step in step 5
and $s$ is the number of repeated times of step 5.

7\. Output the solution at time $T$ in the DeepONet with parameters obtained
in step 6.

##  6 Numerical Experiments

In this section, we implement EDE-DeepONet to solve heat equations, parametric
heat equations, and Allen-Cahn equations to show its performance and accuracy.

###  6.1 Example 1: Simple heat equations

To show the accuracy of the EDE-DeepONet, we start with the simple heat
equation with different initial conditions since we already have the exact
solution. A 1D heat equation system can be described by

$\displaystyle u_{t}=u_{xx}$ (43) $\displaystyle u(x,0)=f$ (44) $\displaystyle
u(0,t)=u(2,t)=0$ (45)

By the method of separation of variables, we can derive the solution to the
heat equation. If we set $f(x)=asin(\pi x)$, the solution is $u(x,t)=asin(\pi
x)e^{-\pi^{2}t}$, where $a\in[1,2]$. The corresponding energy is
$E(u)=\int_{0}^{2}\frac{1}{2}|u_{x}|^{2}dx\approx\Delta
x(\sum_{i=1}^{n}\frac{1}{2}|u_{x}(x_{i})|^{2})$. With different parameters
$a$, the above equation describes a kind of PDE. The input data samples can be
generated as $(a,x,\mathcal{G}(a)(x))$, where $\mathcal{G}(a)(x)=a\sin(\pi x)$
for specific $a$ and $x$. When generating the initial data samples, we choose
50 points from $[0,2)$ uniformly for x and 50 random values of $a$ from
$[1,2]$. The time step when updating the parameters in the neural network is
$2.5\times 10^{4}$. The number of iteration steps is 400. We compared the
different solutions with 4 different $a$, $1.0$, $1.5$, $1.8$, $2.5$ every 100
steps. Although $a=2.5$ is out of the range of training data, it still
performs well in this model. With the exact solution, we also get the error
with different $a$ as Table 1. The error is defined by
$\frac{1}{N_{x}}\sum_{k=1}^{N_{x}}({u(x_{k})-\hat{u}(x_{k})})^{2},$ where
$N_{x}=51$, $u$ is the solution obtained by EDE-DeepONet and $\hat{u}$ is the
exact solution. To illustrate the relationship between the modified energy and
the original energy, we compare $r^{2}$ and $E$ at each step as Figure 2. Both
energy are actually disspative in the EDE-DeepONet except when restart
strategy applied. The restart strategy is used to keep $r^{2}$ approaching
$E$. The modified energy is initialized when the restart strategy applied. The
restart strategy was triggered on the 370th step since the modified energy and
the original energy are offset. After that, they are on the same trajectory
again. It is clear that the modified energy approaches the original energy
before and after the restart strategy applied. In Figure 3, we give the
comparison between the exact solution and the solution obtained by EDE-
DeepONet. From this simple heat equation, we show that EDE-DeepONet correctly
predicts the solution of the PDE. The most important fact is that EDE-DeepONet
can not only predict the solution in the training subset range but also the
solution out of the training range. For instance, we take $a=2.5$ while
$a\in[1,2]$ in the training process. EDE-DeepONet shows good accuracy compared
to the exact solution as Figure 3 (a)-(d) and Table 1.

Fig. 2: The heat equation: The modified energy and original energy when
training the network. Each iteration step represents one forward step of the
PDE’s numerical solution with $\Delta t=2.5\times 10^{-4}$. In the EDE-
DeepONet, both energy are actually dissipative, with the exception of the
restart strategy. In order to maintain that $r^{2}$ approaches $E$, the
restart strategy is employed, and the modified energy is initialized on the
370th step. The modified energy and the original energy on the same trajectory
before and after the 370th step.

(a) Initial $f(x)=sin(\pi x)$

(b) $f(x)=1.5sin(\pi x)$

(c) $f(x)=1.8sin(\pi x)$

(d) $f(x)=2.5sin(\pi x)$

Fig. 3: The heat equation: The solution with 4 different initial conditions
$f(x)=a\sin(\pi x)$. The curve represents the solution obtained by the EDE-
DeepONet, and xxx represents the reference solution. The training parameter
$a$ is in the range of $[1,2)$, so we give three examples in this range. We
also present the case out of the range. It also shows accuracy in Figure
3-(d).

$\begin{array}[]{||c|ccccc||}\hline\cr\text{ Error
}&&T=0.025&T=0.05&T=0.075&T=0.1\\\ \hline\cr\hline\cr a=1.0&&1.47\times
10^{-5}&1.33\times 10^{-5}&1.32\times 10^{-5}&1.29\times 10^{-5}\\\ \hline\cr
a=1.5&&5.11\times 10^{-6}&7.05\times 10^{-6}&8.48\times 10^{-6}&9.81\times
10^{-6}\\\ \hline\cr a=1.8&&1.46\times 10^{-5}&1.69\times 10^{-5}&1.79\times
10^{-5}&1.83\times 10^{-5}\\\ \hline\cr a=2.5&&2.20\times 10^{-4}&1.34\times
10^{-4}&6.02\times 10^{-5}&1.72\times 10^{-5}\\\ \hline\cr\end{array}$

Table 1: The heat equation: The initial condition of the PDE is
$f(x)=a\sin{(\pi x)}$. The error is defined by
$\frac{1}{N_{x}}\sum_{k=1}^{N_{x}}({u(x_{k})-\hat{u}(x_{k})})^{2},$ where
$N_{x}=51$, $u$ is the solution obtained by EDE-DeepONet and $\hat{u}$ is the
exact solution.

###  6.2 Example 2: Parametric heat equations

In example 1, we take different initial conditions as our inputs. In example
2, we are going to deal with the parametric heat equations. A general
parametric heat equation in 1D can be described by

$\displaystyle u_{t}=cu_{xx}$ (46) $\displaystyle u(x,0)=sin(\pi x)$ (47)
$\displaystyle u(0,t)=u(2,t)=0$ (48)

This PDE is more complex than the PDE in Example 1 since the parameter is
inside the equation. The traditional numerical scheme needs to be run multiple
times to deal with the case with different parameters because they are
actually different equations. However, we only need to train the EDE-DeepONet
once. The training range of $c$ is chosen as $[1,2)$. We choose 50 points of
$x$ and $c$ in the same way as example 1. First, we compared the modified
energy with the original energy as Figure 4. The energy is not the same as the
first example since the energy depends on the parameter $c$. We compute the
average of the energy with different $c$ to represent the energy of the
system. This case is more complex than the first one, so it needs more
restarts during the training. Even though the modified energy oscillates when
restart strategy used, it keeps decreasing after each restart. Second, we give
the error between the solution obtained by the EDE-DeepONet and the reference
solution in Table 2, where the reference solution can be obtained explicitly
by variable separation method and the error is defined in the same way as
example 1. Third, we give the comparison between our solution and the
reference solution in Figure 5. Same as example 1, we give the predicted
solution of $c\notin[1,2]$. All of them show the good accuracy. Hence, EDE-
DeepONet can actually solve parametric PDEs.

Fig. 4: The parametric heat equation: The modified energy and original energy
when training the network. Each iteration step represents one forward step of
the PDE’s numerical solution with $\Delta t=2.5\times 10^{-4}$. This kind of
PDEs is more complicated, so it need more restarts in the training process.
The original energy keeps decreasing and the modified energy also shows good
approximation of the original energy.

(a) $c=1.2$

(b) $c=1.5$

(c) $c=1.8$

(d) $c=2.5$

Fig. 5: The parametric heat equation: The solution with 4 different parameters
$c$. The curve represents the solution obtained by the EDE-DeepONet and xxx
represents the reference solution. The training parameter $c$ is in the range
of $[1,2)$, so we give 3 examples in this range. We also present the case out
of the range in Figure 5-(d).

$\begin{array}[]{||c|ccccc||}\hline\cr\text{ Error
}&&T=0.025&T=0.05&T=0.075&T=0.1\\\ \hline\cr\hline\cr c=1.2&&1.30\times
10^{-5}&1.43\times 10^{-5}&1.35\times 10^{-5}&1.20\times 10^{-5}\\\ \hline\cr
c=1.5&&1.35\times 10^{-5}&1.27\times 10^{-5}&9.80\times 10^{-6}&7.80\times
10^{-6}\\\ \hline\cr c=1.8&&1.17\times 10^{-5}&1.03\times 10^{-5}&7.88\times
10^{-5}&1.83\times 10^{-5}\\\ \hline\cr c=2.5&&2.20\times 10^{-4}&1.34\times
10^{-4}&6.02\times 10^{-5}&7.08\times 10^{-6}\\\ \hline\cr\end{array}$

Table 2: The parametric heat equation: The initial condition of the PDE is
$f(x)=\sin{(\pi x)}$. The error is defined by
$\frac{1}{N_{x}}\sum_{k=1}^{N_{x}}({u(x_{k})-\hat{u}(x_{k})})^{2}$, where
$N_{x}=51$, $u$ is the solution obtained by EDE-DeepONet and $\hat{u}$ is the
exact solution.

###  6.3 Example 3: Allen-Cahn equations

The energy in Examples 1 and 2 is quadratic and the right-hand side of the PDE
is linear with respect to $u$. We are going to show the result for the PDE
with more complicated energy. The Allen-Cahn equation is a kind of reaction-
diffusion equation. It is derived to describe the process of the phase
separation. It was developed to solve a problem in the material science area
and has been used to represent the moving interfaces in a phase-field model in
fluid dynamics. The Allen-Cahn equation can be treated as a gradient flow in
$L^{2}$ with some specific energy. We discussed the 1D case and 2D case as
follows:

####  6.3.1 1D case

(a) Various initial conditions:
We start with the simple case, 1D Allen-Cahn equation. It can be described by
the following equations:

$\displaystyle u_{t}=u_{xx}-g(x)$ (49) $\displaystyle u(x,0)=a\sin{\pi x}$
(50) $\displaystyle u(-1,t)=u(1,t)=0$ (51)

The corresponding Ginzburg–Landau free energy
$E[u]=\int_{0}^{1}\frac{1}{2}|u_{x}|^{2}dx+\int_{x=0}^{x=1}G(u)dx$, where
$G(u)=\frac{1}{4\epsilon^{2}}(u^{2}-1)^{2}$ and
$g(u)=G^{\prime}(u)=\frac{1}{\epsilon^{2}}u(u^{2}-1)$, $\epsilon=0.1$. The
parameter $\epsilon$ affects the width of the jump when arriving at the steady
state as the Figure 7 (c), (j), (o) and (t). In the EDE-DeepONet, we set
$\Delta t=10^{-4}$, the number of spatial points $N_{x}$ is 51 and the range
of $a$ is $[0.1,0.5]$. We also compared the modified energy and the original
energy as Figure 6. The modified energy can approximate well to the original
energy even in a much more complicated form. Then, we compared 4 different
solutions with different $a\in[0.1,0.5]$ obtained by the EDE-DeepONet and the
reference solution obtained by the SAV method in traditional numerical
computation as Figure 7. The error is shown in Table 3, where error is defined
in the same way as example 1. $a=0.6\notin[0.1,0.5)$ shows that EDE-DeepONet
can predict the solution well out of the training range. We compared the
solution with 4 different initial condition parameter $a$ every 100 steps
until the final time $T=0.04$ as Figure 7. Each row presents the solution
under the same initial condition but with different evolution time $T$. With
this example, it shows that EDE-DeepONet can deal with the PDE with a jump,
while it is hard for other neural networks.

$\begin{array}[]{||c|ccccc||}\hline\cr\text{ Error
}&&T=0.01&T=0.02&T=0.03&T=0.04\\\ \hline\cr\hline\cr a=0.1&&4.95\times
10^{-5}&2.27\times 10^{-4}&4.97\times 10^{-4}&7.40\times 10^{-4}\\\ \hline\cr
a=0.3&&2.32\times 10^{-4}&4.62\times 10^{-4}&6.66\times 10^{-4}&7.58\times
10^{-4}\\\ \hline\cr a=0.2&&1.25\times 10^{-4}&3.46\times 10^{-4}&4.86\times
10^{-3}&7.15\times 10^{-4}\\\ \hline\cr a=0.4&&2.80\times 10^{-4}&4.83\times
10^{-4}&6.35\times 10^{-4}&7.09\times 10^{-4}\\\ \hline\cr a=0.6&&6.45\times
10^{-4}&5.02\times 10^{-4}&4.90\times 10^{-4}&6.15\times 10^{-4}\\\
\hline\cr\end{array}$

Table 3: 1D Allen-Cahn equation: The initial condition of the Allen-Cahn
equation is $f(x)=a\sin{(\pi x)}$. The error is defined by
$\frac{1}{N_{x}}\sum_{k=1}^{N_{x}}({u(x_{k})-\hat{u}(x_{k})})^{2}$, where
$N_{x}=51$, $u$ is the solution obtained by EDE-DeepONet and $\hat{u}$ is the
reference solution. Fig. 6: 1D Allen-Cahn equation: The modified energy and
original energy when training the network are shown above. Each iteration step
represents one forward step of the PDE’s numerical solution with $\Delta
t=10^{-4}$. The modified energy shows the same trends as the original energy.

(a) $a=0.1$, $T=0$

(b) $a=0.1$, $T=0.01$

(c) $a=0.1$, $T=0.02$

(d) $a=0.1$, $T=0.03$

(e) $a=0.1$, $T=0.04$

(f) $a=0.3$, $T=0$

(g) $a=0.3$, $T=0.01$

(h) $a=0.3$, $T=0.02$

(i) $a=0.3$, $T=0.03$

(j) $a=0.3$, $T=0.04$

(k) $a=0.4$, $T=0$

(l) $a=0.4$, $T=0.01$

(m) $a=0.4$, $T=0.02$

(n) $a=0.4$, $T=0.03$

(o) $a=0.4$, $T=0.04$

(p) $a=0.6$, $T=0$

(q) $a=0.6$, $T=0.01$

(r) $a=0.6$, $T=0.02$

(s) $a=0.6$, $T=0.03$

(t) $a=0.6$, $T=0.04$

Fig. 7: 1d Allen-Cahn equation: The solution for 1d Allen-Cahn equation with 4
different initial conditions $f(x)=a\sin{\pi x}$. The curve represents the
solution obtained by our model, and xxx represents the reference solution. We
draw the figure for every 100 steps. The range of $a$ is $[0.1,0.5]$. We also
compare the solution with $a\notin[0.1,0.5]$. All the figures show the trends
of the phase separation.

(b) Various thickness of the interface:
Heuristically, $\epsilon$ represents the thickness of the interface in the
phase separation process. We are able to obtain a sharp interface when
$\epsilon\rightarrow 0$ with evolving in time. Each theoretical and numerical
analysis of the limit makes a difference in the purpose of the understanding
of the equation, cf. e.g. [34, 35]. We take $\epsilon$ as a training
parameter. The problem can be described as:

$\displaystyle u_{t}=u_{xx}-\frac{1}{\epsilon^{2}}(u^{3}-u)$ (52)
$\displaystyle u(-1,t)=u(1,t)=0$ (53)

Since the training sample contains the parameter $\epsilon$, we can not use
the same initial condition as the last example. We use spectral methods for a
few steps with initial condition $u(x,0)=0.4\sin{(\pi x)}$. The training
sample is generated based on the numerical solution of $u_{\epsilon}(x,0.02)$.
We randomly select 50 different $\epsilon$ from $[0.1,0.2]$. We set the
learning rate as $\Delta t=10^{-4}$ and apply the adaptive time stepping
strategy. We obtain the predicted solution after 400 iterations with different
$\epsilon$. The rest setting is the same as the last example. The solution
with different $\epsilon$ is shown in Figure 8. As $\epsilon$ goes smaller,
the interface is sharper. Besides, the range of the training parameter is
$(0.1,0.2)$. We are also able to obtain the solution out of the above range.
EDE-DeepONet can track the limit of $\epsilon$ in only one training process,
EDE-DeepONet can track the limit of $\epsilon$ in only one training process,
whereas other traditional numerical methods hardly make it.

Fig. 8: 1d Allen-Cahn equation: Solutions with different thickness of the
interface at the same final time. The curve represents the solution obtained
by EDE-DeepONet. xxx represents the reference solution.

####  6.3.2 2D case

The 2D case Allen-Cahn equation is even more complex. The problem can be
described as follows:

$\displaystyle u_{t}=\Delta u-g(u)$ (54) $\displaystyle u(x,y,0)=asin(\pi
x)sin(\pi y)$ (55) $\displaystyle u(-1,y,t)=u(1,y,t)=u(x,-1,t)=u(x,1,t)=0$
(56)

The corresponding Ginzburg–Landau free energy
$E[u]=\int_{-1}^{1}\int_{-1}^{1}\frac{1}{2}(|u_{x}|^{2}+|u_{y}|^{2})dxdy+\int_{-1}^{1}\int_{-1}^{1}G(u)dx$,
where $G(u)=\frac{1}{4\epsilon^{2}}(u^{2}-1)^{2}$ and
$g(u)=G^{\prime}(u)=\frac{1}{\epsilon^{2}}u(u^{2}-1)$. Usually, we take
$\epsilon=0.1$. In the training process, we take $\Delta t=2\times 10^{-4}$.
The number of spatial points is $51\times 51$ and the number of training
parameters $a$ is 20. The way to choose $a\in(0.1,0.4)$ and $x$ is the same as
in example 1. We first compared the exact solution and the solution obtained
by EDE-DeepONet with initial condition $f(x,y)=0.2sin(\pi x)sin(\pi y)$, where
the exact solution is obtained by the traditional SAV method. EDE-DeepONet
predicts the solution correctly based on Table 4 and Figure 9. Then in order
to show its accuracy, we draw Figure 10 with more parameters. All the examples
show good trends to separate. The case $a=0.4$ is out of the training range,
but it still approaches the exact solution.

$\begin{array}[]{||c|cccc||}\hline\cr\text{ Error }&&T=0.01&T=0.02&T=0.03\\\
\hline\cr\hline\cr a=0.15&&1.23\times 10^{-4}&6.53\times 10^{-4}&2.75\times
10^{-3}\\\ \hline\cr a=0.2&&2.24\times 10^{-4}&1.10\times 10^{-3}&4.04\times
10^{-3}\\\ \hline\cr a=0.3&&4.28\times 10^{-4}&1.84\times 10^{-3}&5.76\times
10^{-3}\\\ \hline\cr a=0.35&&5.31\times 10^{-4}&2.17\times 10^{-3}&6.25\times
10^{-3}\\\ \hline\cr a=0.4&&6.94\times 10^{-4}&2.71\times 10^{-3}&7.22\times
10^{-3}\\\ \hline\cr\end{array}$

Table 4: 2D Allen-Cahn equation: The initial condition of the 2D Allen-Cahn
equation is $f(x,y)=a\sin{(\pi x)}\sin{(\pi y)}$. The error is defined by
$\frac{1}{N_{x}N_{y}}\sum_{k=1}^{N_{x}}\sum_{j=1}^{N_{y}}({u(x_{k},y_{j})-\hat{u}(x_{k},y_{j})})^{2},$
where $N_{x}=N_{y}=51$, $u$ is the solution obtained by EDE-DeepONet and
$\hat{u}$ is the reference solution.

(a) $a=0.2$, $T=0$

(b) $a=0.2$, $T=0.01$

(c) $a=0.2$, $T=0.02$

(d) $a=0.2$, $T=0.03$

(e) $a=0.2$, $T=0$

(f) $a=0.2$, $T=0.01$

(g) $a=0.2$, $T=0.02$

(h) $a=0.2$, $T=0.03$

Fig. 9: 2D Allen-Cahn equation: (a)-(d) represents the reference solution of
the 2D Allen-Cahn equation with initial condition $f(x,y)=0.3\sin{(\pi
x)}\sin{(\pi y)}$. (e)-(h) is the solution obtained by the EDE-DeepONet.

(a) $a=0.15$, $T=0$

(b) $a=0.15$, $T=0.01$

(c) $a=0.15$, $T=0.02$

(d) $a=0.15$, $T=0.03$

(e) $a=0.35$, $T=0$

(f) $a=0.35$, $T=0.01$

(g) $a=0.35$, $T=0.02$

(h) $a=0.35$, $T=0.03$

Fig. 10: 2D Allen-Cahn equation: The solution of 2D Allen-Cahn equation with 4
different initial conditions $f(x,y)=a\sin{\pi x}\sin{\pi y}$. The training
parameter $a\in[0.1,0.4]$. We draw three figures where $a$ is in the training
range and one figure where $a$ is out of the training range. All the figures
show the phase separation trends according to the reference solution. As $a$
is further away from the training range, the error tends to be larger.

##  7 Concluding Remarks

In this paper, we provide a new neural network architecture to solve
parametric PDEs with different initial conditions, while maintaining the
energy dissipative of dynamic systems. We first introduce the energy
dissipative law of dynamic systems to the DeepONet. We also introduce an
adaptive time stepping strategy and restart strategy. With our experiments,
both above strategies help keep the modified energy approaching the original
energy. To avoid much cost of training the DeepONet, we evolve the neural
network based on Euler methods. In this article, we adopt the SAV method to
solve gradient flow problems. With this successful attempt, more work could be
done. For example, we can consider a general Wasserstein gradient flow
problem. We are only adopting the basic architecture of the DeepONet. The more
advanced architecture is compatible to our work. It may further improve the
accuracy of EDE-DeepONet.

## Acknowledgments

SJ and SZ gratefully acknowledge the support of NSF DMS-1720442 and AFOSR
FA9550-20-1-0309. GL and ZZ gratefully acknowledge the support of the National
Science Foundation (DMS-1555072, DMS-2053746, and DMS-2134209), Brookhaven
National Laboratory Subcontract 382247, and U.S. Department of Energy (DOE)
Office of Science Advanced Scientific Computing Research program DE-SC0021142
and DE-SC0023161.

## References

  * Kovachki et al. [2021] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, A. Anandkumar,  Neural operator: Learning maps between function spaces,  arXiv preprint arXiv:2108.08481 (2021).
  * Li et al. [2020] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar,  Neural operator: Graph kernel network for partial differential equations,  arXiv preprint arXiv:2003.03485 (2020).
  * Khoo et al. [2021] Y. Khoo, J. Lu, L. Ying,  Solving parametric pde problems with artificial neural networks,  European Journal of Applied Mathematics 32 (2021) 421–435.
  * Bhattacharya et al. [2020] K. Bhattacharya, B. Hosseini, N. B. Kovachki, A. M. Stuart,  Model reduction and neural networks for parametric pdes,  arXiv preprint arXiv:2005.03180 (2020).
  * Nelsen and Stuart [2021] N. H. Nelsen, A. M. Stuart,  The random feature model for input-output maps between banach spaces,  SIAM Journal on Scientific Computing 43 (2021) A3212–A3243.
  * Li et al. [2020] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar,  Fourier neural operator for parametric partial differential equations,  arXiv preprint arXiv:2010.08895 (2020).
  * Patel et al. [2021] R. G. Patel, N. A. Trask, M. A. Wood, E. C. Cyr,  A physics-informed operator regression framework for extracting data-driven continuum models,  Computer Methods in Applied Mechanics and Engineering 373 (2021) 113500\.
  * Opschoor et al. [2020] J. A. Opschoor, C. Schwab, J. Zech,  Deep learning in high dimension: Relu network expression rates for bayesian pde inversion,  SAM Research Report 2020 (2020).
  * Schwab and Zech [2019] C. Schwab, J. Zech,  Deep learning in high dimension: Neural network expression rates for generalized polynomial chaos expansions in uq,  Analysis and Applications 17 (2019) 19–55.
  * O’Leary-Roseberry et al. [2022] T. O’Leary-Roseberry, U. Villa, P. Chen, O. Ghattas,  Derivative-informed projected neural networks for high-dimensional parametric maps governed by pdes,  Computer Methods in Applied Mechanics and Engineering 388 (2022) 114199\.
  * Wu and Xiu [2020] K. Wu, D. Xiu,  Data-driven deep learning of partial differential equations in modal space,  Journal of Computational Physics 408 (2020) 109307.
  * Lu et al. [2021] L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis,  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators,  Nature Machine Intelligence 3 (2021) 218–229.
  * Cybenko [1989] G. Cybenko,  Approximation by superpositions of a sigmoidal function,  Mathematics of control, signals and systems 2 (1989) 303–314.
  * Hornik [1991] K. Hornik,  Approximation capabilities of multilayer feedforward networks,  Neural networks 4 (1991) 251–257.
  * Hornik et al. [1989] K. Hornik, M. Stinchcombe, H. White,  Multilayer feedforward networks are universal approximators,  Neural networks 2 (1989) 359–366.
  * Chen and Chen [1995] T. Chen, H. Chen,  Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems,  IEEE Transactions on Neural Networks 6 (1995) 911–917.
  * Raissi et al. [2018] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Multistep neural networks for data-driven discovery of nonlinear dynamical systems,  arXiv preprint arXiv:1801.01236 (2018).
  * Qin et al. [2021] T. Qin, Z. Chen, J. D. Jakeman, D. Xiu,  Deep learning of parameterized equations with applications to uncertainty quantification,  International Journal for Uncertainty Quantification 11 (2021).
  * Winovich et al. [2019] N. Winovich, K. Ramani, G. Lin,  Convpde-uq: Convolutional neural networks with quantified uncertainty for heterogeneous elliptic partial differential equations on varied domains,  Journal of Computational Physics 394 (2019) 263–279.
  * Zhu et al. [2019] Y. Zhu, N. Zabaras, P.-S. Koutsourelakis, P. Perdikaris,  Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data,  Journal of Computational Physics 394 (2019) 56–81.
  * del Águila Ferrandis et al. [2021] J. del Águila Ferrandis, M. S. Triantafyllou, C. Chryssostomidis, G. E. Karniadakis,  Learning functionals via lstm neural networks for predicting vessel dynamics in extreme sea states,  Proceedings of the Royal Society A 477 (2021) 20190897.
  * Jia and Benson [2019] J. Jia, A. R. Benson,  Neural jump stochastic differential equations,  Advances in Neural Information Processing Systems 32 (2019).
  * Chen et al. [2018] T. Chen, Y. Rubanova, J. Bettencourt, D. Duvenaud,  Neural ordinary differential equations, in ‘advances in neural information processing systems’,  La Jolla (2018).
  * Du and Zaki [2021] Y. Du, T. A. Zaki,  Evolutional deep neural network,  Physical Review E 104 (2021) 045303.
  * Allen and Cahn [1979] S. M. Allen, J. W. Cahn,  A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening,  Acta metallurgica 27 (1979) 1085–1095.
  * Anderson et al. [1998] D. M. Anderson, G. B. McFadden, A. A. Wheeler,  Diffuse-interface methods in fluid mechanics,  Annual review of fluid mechanics 30 (1998) 139–165.
  * Cahn and Hilliard [1958] J. W. Cahn, J. E. Hilliard,  Free energy of a nonuniform system. i. interfacial free energy,  The Journal of chemical physics 28 (1958) 258–267.
  * Doi et al. [1988] M. Doi, S. F. Edwards, S. F. Edwards, The theory of polymer dynamics, volume 73, oxford university press, 1988.
  * Elder et al. [2002] K. Elder, M. Katakowski, M. Haataja, M. Grant,  Modeling elasticity in crystal growth,  Physical review letters 88 (2002) 245701.
  * Gurtin et al. [1996] M. E. Gurtin, D. Polignone, J. Vinals,  Two-phase binary fluids and immiscible fluids described by an order parameter,  Mathematical Models and Methods in Applied Sciences 6 (1996) 815–831.
  * Leslie [1979] F. M. Leslie,  Theory of flow phenomena in liquid crystals,  in: Advances in liquid crystals, volume 4, Elsevier, 1979, pp. 1–81.
  * Yue et al. [2004] P. Yue, J. J. Feng, C. Liu, J. Shen,  A diffuse-interface method for simulating two-phase flows of complex fluids,  Journal of Fluid Mechanics 515 (2004) 293–317.
  * Shen et al. [2018] J. Shen, J. Xu, J. Yang,  The scalar auxiliary variable (sav) approach for gradient flows,  Journal of Computational Physics 353 (2018) 407–416.
  * Caginalp and Chen [1998] G. Caginalp, X. Chen,  Convergence of the phase field model to its sharp interface limits,  European Journal of Applied Mathematics 9 (1998) 417–445.
  * Chen et al. [2006] X. Chen, G. Caginalp, C. Eck,  A rapidly converging phase field model,  Discrete & Continuous Dynamical Systems 15 (2006) 1017.


# Evaluating Pedagogical Incentives in Undergraduate Computing: A Mixed
Methods Approach Using Learning Analytics
††thanks: This project received funding from UCL ChangeMakers.

Laura J. Johnston  Department of Statistical Science
University College London
London, UK
<EMAIL_ADDRESS>Takoua Jendoubi  Department of Statistical
Science
University College London
London, UK
<EMAIL_ADDRESS>
###### Abstract

In the context of higher education’s evolving dynamics post-COVID-19, this
paper assesses the impact of new pedagogical incentives implemented in a
first-year undergraduate computing module at University College London. We
employ a mixed methods approach, combining learning analytics with qualitative
data, to evaluate the effectiveness of these incentives on increasing student
engagement.

A longitudinal overview of resource interactions is mapped through Bayesian
network analysis of Moodle activity logs from 204 students. This analysis
identifies early resource engagement as a predictive indicator of continued
engagement while also suggesting that the new incentives disproportionately
benefit highly engaged students. Focus group discussions complement this
analysis, providing insights into student perceptions of the pedagogical
changes and the module design. These qualitative findings underscore the
challenge of sustaining engagement through the new incentives and highlight
the importance of communication in blended learning environments.

Our paper introduces an interpretable and actionable model for student
engagement, which integrates objective, data-driven analysis with students’
perspectives. This model provides educators with a tool to evaluate and
improve instructional strategies. By demonstrating the effectiveness of our
mixed methods approach in capturing the intricacies of student behaviour in
digital learning environments, we underscore the model’s potential to improve
online pedagogical practices across diverse educational settings.

###### Index Terms:

educational technology, curriculum development, computer science education,
learning management systems, statistical learning

##  I Introduction

In the wake of the COVID-19 pandemic, maintaining active and sustained student
engagement has emerged as a paramount challenge in higher education. This
challenge has necessitated a critical reassessment of module designs and the
development of strategies that foster student engagement to support effective
learning. Assessing these design changes requires innovative methods to
measure their impact on student engagement and learning accurately.

Responding to this broader challenge, the Department of Statistical Science at
University College London has introduced new pedagogical incentives across its
modules. During the pandemic, the first-year computing module, STAT0004:
Introduction to Practical Statistics, adopted a flipped classroom model to
promote independent learning. Students individually navigated through lecture
notes, videos, and quizzes before applying their knowledge in live labs, where
they could seek support from the lecturer while practising their coding
skills. The module has since been updated with weekly code submissions and
weekly feedback, aiming to increase student engagement. Introducing these
pedagogical incentives underscores the need to develop an evaluation framework
to assess their impact on enhancing student engagement in the module.

A notable gap exists in the effective quantification and modelling of student
engagement, which is pivotal when evaluating the impact of pedagogical
changes. While student log data extracted from learning management systems,
such as Moodle, provide a mechanism to track online student engagement [1],
[2] caution that reliance on such data alone may overlook critical aspects of
the student experience. Recognising that student engagement is multifaceted
and challenging to capture comprehensively [3], our study adopts a mixed
methods approach integrating advanced learning analytics with focus group
discussions. The novel application of Bayesian network analysis to Moodle log
data enables the dynamic tracking of student interactions with resources over
the term, offering a temporal perspective on engagement. This method, enriched
by qualitative insights from focus group discussions, presents a holistic
model of student engagement.

Using the STAT0004 computing module as a case study, we showcase the potential
of this framework to yield actionable insights, thus equipping educators with
the tools to effectively evaluate and inform the adaptation of their teaching
strategies to enhance student engagement. The primary objectives of this
research are to: (a) assess the impact of the pedagogical changes in STAT0004;
(b) evaluate the overall effectiveness of the STAT0004 module design; and (c)
demonstrate a new framework educators can utilise for module evaluations.

We structure this paper into three main sections: II, which reviews relevant
literature; III, detailing the mixed methods approach of this study; and IV,
presenting and interpreting the findings in the context of STAT0004 while
considering their implications for educational practice and research in higher
education.

##  II Related Works

###  II-A Mixed methods in educational research

Adopting mixed methods in educational research offers a comprehensive approach
to understanding complex pedagogical phenomena [4]. By integrating
quantitative and qualitative data, researchers can capture objective data-
driven insights and subjective experiences that often elude quantitative
analyses [5]. Despite its proven effectiveness in evaluating module design
interventions [6] and assessing the influence of pedagogical changes on
student engagement [7], the adoption of mixed methods in higher education
research is surprisingly limited [8].

Furthermore, the depth and sophistication of learning analytics tools often
fall short in mixed methods studies. Traditional quantitative approaches, such
as Likert scale surveys [9], exam scores [10], and attendance and login data
[7], though helpful, only provide surface-level insights. Similarly,
qualitative methods, including focus groups [11], structured interviews [12],
and open-ended survey questions [6], while insightful for capturing student
perceptions, often lack robust quantitative support. This study aims to bridge
this gap by combining advanced learning analytics with in-depth qualitative
insights. We expect this approach to yield a comprehensive and nuanced
understanding of student engagement, paving the way for more effective
educational strategies.

###  II-B Bayesian network

A Bayesian network is a statistical model representing a set of outcomes
(nodes in the graph) and their conditional dependencies (arrows in the graph)
using a graphical structure informed by observed data. In the context of
educational research, these networks can estimate the probability of various
student outcomes, such as behaviours and academic performance. When an arrow
connects two nodes, it indicates that knowledge of one outcome provides
information that can alter the assessed likelihood of the other outcome
occurring without implying direct causation. The adaptability of Bayesian
networks across different module designs and educational settings [13]
highlights their methodological versatility. Furthermore, these networks are
visual and accessible, allowing educators of diverse statistical backgrounds
to interpret the data.

Bayesian networks have been instrumental in addressing complex educational
challenges such as building adaptive tutoring systems [14], modelling learning
styles [15], and predicting student dropout [16]. Similarly, [17] applied
dynamic Bayesian networks to track student interactions in engineering design
software. In this study, we employ a dynamic Bayesian network to analyse
longitudinal engagement patterns in Moodle log data across instructional
resources. This broadens the scope of Bayesian network applications in
educational research.

##  III Methodology

###  III-A Context

This study focuses on the STAT0004 module over the 2022–2023 academic year.
While the module spans two terms, totalling 20 weeks, we restrict attention to
the first term to expedite evaluating the new changes. This term encompasses
ten weeks, with the curriculum covering nine chapters in weeks 1–9 and
concluding with a graded assessment in week 10. Live coding labs deliver the
module, supplemented by learning resources accessed via Moodle. Most chapters
include lecture notes containing content and coding tasks, explanatory videos,
a formative quiz, and a code submission point due weekly. Though actively
encouraged, participation in the quiz and code submissions is optional, with
no direct impact on students’ grades. Feedback on the code submissions is
verbally provided during the subsequent live lab.

###  III-B Learning analytics

The objective of the learning analytics component of this study is to
quantitatively assess student engagement patterns within the STAT0004 module
to evaluate the impact of the new pedagogical incentives. To this end, the
analysis leverages a dataset extracted from the Moodle activity logs for
STAT0004, detailing the interactions of 204 students with the module’s online
resources during the first term. This dataset offers a granular view of how
students engaged with the module’s activities, such as accessing lecture
notes, watching videos, completing quizzes, and making submissions. Notably,
we cannot capture the impact of the weekly feedback in the log data.

In the dataset, each row records a student’s click on the Moodle page for the
STAT0004 module, including the time and nature of the interaction. We focused
our analysis on synchronous interactions, defined as actions occurring within
two weeks of the resource’s release date. This was designed to capture
regular, day-to-day student interactions with module resources and exclude
periods of atypical engagement, such as pre-assessment revision, to align with
the study’s intentions of evaluating active and sustained student engagement.
Synchronous actions were identified using the Time column while the Event
context and Event name columns discerned the type of resource interaction. We
converted the logs into a binary table, where each entry indicates if a
student interacted with a specific resource within the defined synchronous
period.

We fitted the Bayesian network using the bnlearn package in R [18]. The model
was constrained to allow stationary or forward connections, meaning the
network could only establish relationships between resources of the same
chapter or subsequent chapters in line with the chronological order of the
curriculum. Lecture notes from chapters 1-3 were excluded from the model as
over 95% of students accessed them, limiting their analytical value. We used
bootstrapping to enhance the reliability of the Bayesian network and mitigate
the risk of overfitting. This validation process fitted the model 100 times,
each with a different subset of the data. We retained connections in at least
50% of these iterations in the final model, ensuring it included only the most
stable and significant relationships between module resources.

The resulting network graphically represents the probabilistic dependencies
between resources, offering insights into student behaviours concurrent with
the pedagogical changes.

###  III-C Focus group

We designed the focus group methodology to complement the quantitative data,
gathering qualitative insights on students’ perceptions of the STAT0004 module
and the new pedagogical incentives. All undergraduate students from the
Department of Statistical Science were invited via email to participate,
encapsulating current and former students of the STAT0004 module. Though
initially ten students expressed interest, the focus group comprised five
participants, representing a blend of study years, genders, and nationalities,
each receiving a £20 compensation for their time.

We structured the focus group into three sections, each led by a different
project team member. The sections sought to: (1) gather general feedback on
the STAT0004 module; (2) collect specific feedback on the newly introduced
pedagogical incentives; and (3) encourage student input in future module
design changes. This structure, while incorporating feedback on the
pedagogical changes, also enabled an assessment of current student experiences
and promoted student involvement in shaping future educational strategies.
Ethical considerations, including informed consent and confidentiality, were
rigorously observed.

##  IV Results and Discussion

Fig. 1 visualises the Bayesian network analysis of the Moodle activity data.
In this graph, each node represents a specific action with a module resource
within STAT0004: accessing lecture notes, watching videos, completing quizzes,
or submitting code. The arrows between nodes represent conditional
dependencies, suggesting how engagement with one resource may relate to
engagement with another. These dependencies offer insights into possible
relationships rather than direct causation, allowing for an intuitive
understanding of patterns in student interactions with the module content.

We interpret and discuss these results in tandem with the focus group
discussions, which offer a qualitative dimension to our understanding of
student engagement in STAT0004.

Figure 1: Bayesian Network Depicting Resource Interactions Over Chapters for
the STAT0004 Module. This graph illustrates the conditional dependencies
between each module chapter’s resources (nodes). Nodes are colour–coded by
type: blue for quizzes (quiz_1 to quiz_9), green for videos (vid_1 to vid_9),
orange for lecture notes (ln_4 to ln_9), and pink for submissions (sub_1 to
sub_8), where the suffix indicates the chapter number. Directed arrows between
nodes represent probabilistic dependencies, where the engagement with one
resource may influence the likelihood of engagement with another. Notably,
lecture notes from chapters 1 to 3 are absent due to uniform high access
rates, while other missing resources were not available.

###  IV-A Weekly code submissions

The Bayesian network graph shows that the nodes representing code submissions
for Chapters 1-5, which cover the first half of the term, are disconnected.
This observation implies that participation in submissions in any of these
early weeks does not necessarily increase the likelihood of submitting work in
the following weeks alone. However, this changes after the midpoint, where 80%
of students who submitted code for Chapter 5 also submitted code for Chapter
6, and 85% of students who submitted code for Chapters 6 and 7 did the same
for Chapter 8. Additionally, in the latter half of the term, 83% of students
who watched the video for Chapter 7 also submitted their code.

This insight is supported by the focus group, who found the weekly submissions
in STAT0004 initially engaging but noted a decline in participation over time,
as these tasks lacked direct implications for grades. Participant A said,
“There’s no consequence; you don’t need to submit it. It’s just more like
holding yourself accountable.” This observation suggests that while students
initially engage with the submissions, their commitment may diminish upon
realising their optional nature. To counter this, Participant A suggested that
“…every student cares about their grades…” advocating to integrate submissions
into the grading system, with grades awarded for effort, to sustain
engagement.

These findings suggest that the introduction of weekly code submissions has
successfully fostered initial active engagement among students. However, these
submissions have only maintained the interest of highly engaged students
interacting across a broad spectrum of module resources, such as watching
weekly videos. This pattern indicates that while these submissions retain the
engagement of active participants, they might not have successfully sustained
or increased engagement among those displaying lower levels of engagement.
This phenomenon could be interpreted as a novel application of the ‘Matthew
Effect’ [19], aligning with research indicating that educational interventions
often benefit already engaged or higher-achieving students more significantly
[20]. Therefore, it is imperative to implement targeted strategies that
captivate and sustain the engagement of all students, especially those at risk
of disengaging.

###  IV-B Weekly feedback

While we could not evaluate the impact of the weekly feedback through the
Bayesian network analysis, the focus group revealed concerns about the
visibility and effectiveness of this new incentive. Participant B highlighted
a communication gap: “To be honest, I’m not aware of any weekly feedback.” and
suggested making feedback channels more accessible by “…posting it on Moodle,
rather than only discussing it in the lab.”

This insight underscores that providing feedback alone is insufficient; it
must be effectively communicated and readily accessible to students. The
recommendation to employ Moodle for disseminating feedback aligns with [21],
highlighting the critical role of familiar platforms for communication to
promote student engagement. Adopting this strategy not only aims to address
the current communication gap but also supports the integration of feedback
into future quantitative evaluations for a more informed analysis of this
pedagogical incentive.

###  IV-C Module design

The focus group discussions reveal the pivotal role of early and consistent
engagement with module resources to remain aligned with the demanding
curriculum. Participants emphasised the challenges of independent study and
the necessity of self-directed learning strategies and described the module as
intensive. Participant C noted, “This is the kind of module that needs a large
amount of time to study by ourselves”, while others labelled the module super
hard. Students identified extensive reading requirements as a barrier to
engagement. Participant D shared, “It takes a long time for me to actually sit
down and read the paragraph at the beginning.”

The Bayesian network analysis extended this finding by identifying a pattern
of sustained engagement: students engaging with the quiz or video resources in
one chapter were likely to continue to engage with this same resource in
subsequent chapters. For example, 82% of students who viewed the video for
Chapter 3 also accessed the video for Chapter 4. Conversely, students who
initially did not engage with a particular resource showed a marked decrease
in subsequent participation. Only 14% of those who skipped the Chapter 1 quiz
engaged with the Chapter 2 quiz, and a mere 7% of students who avoided the
Chapter 2 quiz went on to attempt the quiz for Chapter 3.

These findings underscore the importance of initial engagement with the
resources. Extending the findings of [22], which establishes early engagement
as a predictor of future behaviour and performance, these insights demonstrate
that early engagement with a specific resource might indicate future
interaction with that same resource. This observation suggests students may
develop routines early, exhibiting reluctance to explore new resources mid-
term. Consequently, it becomes imperative to consider strategies that
encourage diversified engagement from the onset. This insight might also
suggest that interventions to identify and support students showing signs of
low initial engagement can be implemented early.

There were several further interesting observations about the module design.
From the Bayesian network analysis, quiz interactions are notably isolated,
only linking within themselves. This indicates students’ interaction with
quizzes was independent of the other resources, which might reflect quiz
completion outside the synchronous time frame or some inherently independent
nature of quizzes in the module. This aspect warrants further investigation
into how different assessment types are perceived and impact student
engagement.

Moreover, participants in the focus group recommended adding live lectures to
complement the lab sessions and resources. Participant E compared STAT0004
with a second-year computing module: “…we have a lecture for that course… I
feel like this course is better than [STAT0004].” Others echoed this,
indicating a preference for combining theoretical lectures with practical
workshops. This feedback highlights an actionable insight for enhancing module
design and is a dimension not directly captured in the analytical analysis.

###  IV-D Limitations

A notable limitation of this study lies in the sensitivity of the Bayesian
network structure to the specific definition of synchronous interactions. The
current definition considers interactions within two weeks of release, which
was chosen to model regular engagement with module resources. However, even
minor adjustments to this definition can result in changes to the structure of
the Bayesian network, indicating a lack of robustness to variations in the
definition of synchronous interactions. This sensitivity highlights the need
for careful consideration in defining engagement periods and suggests that
interpretations of the findings must account for this potential variability.

Additionally, despite offering valuable insights, the qualitative aspect of
this study is limited by the small number of participants for the focus group.
With only five participants, there’s a risk the findings may not be
representative of the student population. Future research should include a
more extensive and varied group of participants or implement more inclusive
qualitative methods, such as surveys, to enhance the reliability of the
results from the qualitative data.

##  V Conclusion

Evaluating and refining pedagogical strategies should be an ongoing process to
ensure the effectiveness of educational delivery. This paper introduces a
framework that facilitates this, blending advanced learning analytics with
qualitative insights to model student engagement. This framework’s versatility
allows its adaptation to a wide range of educational contexts, resources, and
module structures, making it a valuable tool for educators seeking to review
and improve their pedagogical approaches using a data-driven and student-
centred approach.

Key findings from the STAT0004 module include the necessity of targeted
pedagogical incentives for diverse engagement levels and the critical role of
effective communication in blended learning environments. Our analysis reveals
that early engagement with specific resources is a reliable indicator of
future participation, underscoring the importance of early and inclusive
engagement strategies.

Building upon the insights gained from the STAT0004 module, our future work
will extend the application of this framework to a broader range of modules
within the Department of Statistical Science, encompassing both computing and
non-computing subjects. This expansion aims to determine the transferability
of our findings across different educational contexts and explore how specific
aspects of module design influence student engagement.

## Acknowledgments

We want to thank Ziqing Li and Jingcheng Song for their invaluable assistance
as research support students in this project and Prof. Jim Griffin and Prof.
Ioanna Manolopoulou for their ongoing supervision and support.

## References

  * [1] A. Pardo, _Design of Learning Analytics Experiences_. Springer New York, 01 2014, pp. 15–38.
  * [2] A. Wong and S. Chong, “Modelling adult learners’ online engagement behaviour: proxy measures and its application,” _Journal of Computers in Education_ , vol. 5, no. 4, pp. 463–479, 2018. [Online]. Available: https://doi.org/10.1007/s40692-018-0123-z
  * [3] G. M. Sinatra, B. C. Heddy, and D. Lombardi, “The challenges of defining and measuring student engagement in science,” _Educational Psychologist_ , vol. 50, no. 1, pp. 1–13, 2015. [Online]. Available: https://doi.org/10.1080/00461520.2014.1002924
  * [4] Q. Nguyen, B. Rienties, and D. Whitelock, “A mixed-method study of how instructors design for learning in online and distance education,” _Journal of Learning Analytics_ , vol. 7, no. 3, pp. 64–78, Dec. 2020. [Online]. Available: https://learning-analytics.info/index.php/JLA/article/view/6807
  * [5] R. B. Johnson, A. J. Onwuegbuzie, and L. A. Turner, “Toward a definition of mixed methods research,” _Journal of Mixed Methods Research_ , vol. 1, no. 2, pp. 112–133, 2007. [Online]. Available: https://doi.org/10.1177/1558689806298224
  * [6] H. Ö. Egilsdottir, L. G. Heyn, E. A. Brembo, K. R. Byermoen, A. Moen, and H. Eide, “The value of a redesigned clinical course during covid-19 pandemic: an explorative convergent mixed-methods study,” _BMC Nursing_ , vol. 21, no. 1, p. 94, 2022. [Online]. Available: https://doi.org/10.1186/s12912-022-00872-8
  * [7] L. Lin, _Student learning and engagement in a blended environment: a mixed methods study_. IGI Global, 05 2018.
  * [8] O. Viberg, M. Hatakka, O. Bälter, and A. Mavroudi, “The current landscape of learning analytics in higher education,” _Computers in Human Behavior_ , vol. 89, pp. 98–110, 2018. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0747563218303492
  * [9] B. E. Sam Muir, Loredana Tirlea and M. Huynh, “Promoting classroom engagement through the use of an online student response system: a mixed methods analysis,” _Journal of Statistics Education_ , vol. 28, no. 1, pp. 25–31, 2020. [Online]. Available: https://doi.org/10.1080/10691898.2020.1730733
  * [10] E. M. Ahlin, “A mixed-methods evaluation of a hybrid course modality to increase student engagement and mastery of course content in undergraduate research methods classes,” _Journal of Criminal Justice Education_ , vol. 32, no. 1, pp. 22–41, 2021. [Online]. Available: https://doi.org/10.1080/10511253.2020.1831034
  * [11] K. Nomme and G. Birol, “Course redesign: an evidence-based approach,” _Canadian Journal for the Scholarship of Teaching and Learning_ , vol. 5, pp. 1–28, 11 2014.
  * [12] Z. Li and J. Li, “Using the flipped classroom to promote learner engagement for the sustainable development of language skills: a mixed-methods study,” _Sustainability_ , vol. 14, no. 10, 2022. [Online]. Available: https://www.mdpi.com/2071-1050/14/10/5983
  * [13] R. Reichenberg, “Dynamic bayesian networks in educational measurement: reviewing and advancing the state of the field,” _Applied Measurement in Education_ , vol. 31, no. 4, pp. 335–350, 2018. [Online]. Available: https://doi.org/10.1080/08957347.2018.1495217
  * [14] T. Käser, S. Klingler, A. G. Schwing, and M. Gross, “Dynamic bayesian networks for student modeling,” _IEEE Transactions on Learning Technologies_ , vol. 10, no. 4, pp. 450–462, 2017.
  * [15] C. Carmona, G. Castillo, and E. Millán, “Designing a dynamic bayesian network for modeling students’ learning styles,” in _2008 Eighth IEEE International Conference on Advanced Learning Technologies_ , 2008, pp. 346–350.
  * [16] A. I. M. Carmen Lacave and J. A. Cruz-Lemus, “Learning analytics to identify dropout factors of computer science studies through bayesian networks,” _Behaviour & Information Technology_, vol. 37, no. 10-11, pp. 993–1007, 2018. [Online]. Available: https://doi.org/10.1080/0144929X.2018.1485053
  * [17] W. Xing, C. Li, G. Chen, X. Huang, J. Chao, J. Massicotte, and C. Xie, “Automatic assessment of students’ engineering design performance using a bayesian network model,” _Journal of Educational Computing Research_ , vol. 59, no. 2, pp. 230–256, 2021. [Online]. Available: https://doi.org/10.1177/0735633120960422
  * [18] M. Scutari, “Learning bayesian networks with the bnlearn R package,” _Journal of Statistical Software_ , vol. 35, no. 3, pp. 1–22, 2010.
  * [19] D. Rigney, _The Matthew effect: how advantage begets further advantage_. Columbia University Press, 2010\. [Online]. Available: http://www.jstor.org/stable/10.7312/rign14948
  * [20] B. Otto and S. Kistner, “Is there a matthew effect in self-regulated learning and mathematical strategy application? - assessing the effects of a training program with standardized learning diaries,” _Learning and Individual Differences_ , vol. 55, pp. 75–86, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1041608017300663
  * [21] K. Betts, “Lost in translation: importance of effective communication in online education,” _Online Journal of Distance Learning Administration_ , vol. 12, 2009. [Online]. Available: https://api.semanticscholar.org/CorpusID:60957719
  * [22] H. E. H. Robert J. Summers and E. Moores, “Measures of engagement in the first three weeks of higher education predict subsequent activity and attainment in first year undergraduate students: a uk case study,” _Assessment & Evaluation in Higher Education_, vol. 46, no. 5, pp. 821–836, 2021. [Online]. Available: https://doi.org/10.1080/02602938.2020.1822282