# Simulating Shock Solutions for Hyperbolic Conservation Laws with Physics-Informed Neural Networks

Mathilde Pascal, Régis Duvigneau, Laurent Monasse

<figure><img src="image_1.png" /></figure>

To cite this version:

Mathilde Pascal, Régis Duvigneau, Laurent Monasse. Simulating Shock Solutions for Hyperbolic Conservation Laws with Physics-Informed Neural Networks. RR-9605, Inria & Université Cote d'Azur, Sophia Antipolis, France. 2025. hal-05405130

HAL Id: hal-05405130

https://inria.hal.science/hal-05405130v1

Submitted on 8 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Inria

Simulating Shock
Solutions for Hyperbolic
Conservation Laws with
Physics-Informed Neural
Networks

Mathilde Pascal, Régis Duvigneau, Laurent Monasse

RESEARCH
REPORT

Nº 9605

December 2025

Project-Team Acumes
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

# Simulating Shock Solutions for Hyperbolic Conservation Laws with Physics-Informed Neural Networks

Mathilde Pascal*, Régis Duvigneau*, Laurent Monasse*

Project-Team Acumes

Research Report n° 9605 — December 2025 — 40 pages

**Abstract:** Physics-informed neural networks represent a new paradigm for modeling complex systems, relying on both data and physical laws (e.g., partial differential equations) for the training. This duality is particularly interesting when little data is available. However, this formulation relies on smooth representations to describe the solutions, so it is not well suited to solving problems where discontinuities occur, as can be the case with nonlinear hyperbolic conservation laws. In this work, we explore how to solve this type of problems using multiple networks, describing on the one hand the behavior of the system in regular regions and on the other hand the jump conditions between these regions. The methods developed are applied to the Burgers equation with shocks.

**Key-words:** Neural networks, hyperbolic conservation laws, discontinuity

<footnote>* Centre Inria d'Université Côte d'Azur, Acumes team</footnote>
Simulation de solutions avec chocs pour des lois de
conservation hyperboliques avec des réseaux neuronaux
informés par la physique

**Résumé** : Les réseaux neuronaux informés par la physique constituent un nouveau paradigme pour la modélisation des systèmes complexes, s’appuyant à la fois sur des données et des lois physiques (par exemple des équations aux dérivées partielles) pour leur entrainement. Cette dualité est particulièrement intéressante lorsque peu de données sont disponibles. Néanmoins, cette formulation s’appuie sur des représentations lisses pour décrire les solutions, elle est donc peu adaptée à la résolution de problèmes pour lesquels des discontinuités apparaissent, comme cela peut être le cas avec des lois de conservation hyperboliques non-linéaires. Dans ce travail, on explore comment résoudre ce type de problèmes en utilisant plusieurs réseaux, décrivant d’une part le comportement du système dans les régions régulières et d’autre part les conditions de saut entre ces régions. On applique les méthodes développées à l’équation de Burgers avec chocs.

Mots-clés : Réseaux neuronaux ; lois de conservation hyperboliques ; discontinuité
# Contents

<table><tr><td>1</td><td>Introduction</td><td>4</td></tr><tr><td>2</td><td>Physics-informed neural networks</td><td>5</td></tr><tr><td>2.1</td><td>Background</td><td>5</td></tr><tr><td>2.2</td><td>Neural networks and discontinuities</td><td>7</td></tr><tr><td>2.2.1</td><td>Universal Approximation Theorems</td><td>7</td></tr><tr><td>2.2.2</td><td>PINNs and discontinuities</td><td>7</td></tr><tr><td>2.2.3</td><td>The Non-Diffusive Neural Network method</td><td>8</td></tr><tr><td>2.2.4</td><td>A modified NDNN algorithm</td><td>10</td></tr><tr><td>3</td><td>Results</td><td>11</td></tr><tr><td>3.1</td><td>1D Burgers equation with continuous solution</td><td>12</td></tr><tr><td>3.2</td><td>Riemann problem for the Burgers equation</td><td>12</td></tr><tr><td>3.3</td><td>One shock wave, one rarefaction wave</td><td>15</td></tr><tr><td>3.4</td><td>Shock wave generation</td><td>21</td></tr><tr><td>4</td><td>Conclusion</td><td>24</td></tr><tr><td>5</td><td>Annex 1: Julia for Scientific Machine Learning</td><td>27</td></tr><tr><td>6</td><td>Annex 2: Expressivity tests</td><td>29</td></tr><tr><td>7</td><td>Annex 3: Supplementary figures for Experiment 1 (one shock wave, one rarefaction wave)</td><td>37</td></tr></table>
# 1 Introduction

Artificial neural networks are one of the main algorithms that underpin modern machine learning methods. With the increasing availability of large datasets and computing power, neural networks have made rapid progress possible in fields such as image processing [8], computer vision [28], and signal processing [20]. For some scientific and engineering applications, however, collecting enough data to train models can be impractical, if not impossible. To bypass the need for extensive datasets, Physics-Informed Neural Networks (PINNs, also sometimes referred to as Physics-Informed Machine Learning - PIML) embed physics laws into the models [24]. Such models can be used where there is little or no data available. They also have the advantage of being mesh-free and can thus offer an alternative or complementary approach to traditional numerical methods such as the Finite Element Method.

In this report, we examine the use of PINNs to solve hyperbolic conservation laws with shock solutions, using the Julia programming language. The following definitions and equations are taken from [27].

**Definition 1 (Hyperbolic conservation law [27])** A conservation law is a system of partial differential equations of the form

$$ \frac{\partial u_i}{\partial t} + \sum_{j=1}^{m} a_{ij}(x, t, u_1, \dots, u_m) \frac{\partial u_j}{\partial x} + b_i(x, t, u_1, \dots, u_m) = 0 $$

that can be written as a flux

$$ \mathbf{U}_t + \mathbf{F}(\mathbf{U})_x = 0 $$

where

$$ \mathbf{U} = \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_m \end{bmatrix}, \quad \mathbf{F}(\mathbf{U}) = \begin{bmatrix} f_1(u_1, \dots, u_m) \\ f_2(u_1, \dots, u_m) \\ \vdots \\ f_m(u_1, \dots, u_m) \end{bmatrix} $$

If, moreover, the Jacobian of the flux function $\mathbf{F}(\mathbf{U})$

$$ \mathbf{A}(\mathbf{U}) = \frac{\partial \mathbf{F}}{\partial \mathbf{U}} = \begin{bmatrix} \frac{\partial f_1}{\partial u_1} & \dots & \frac{\partial f_1}{\partial u_m} \\ \frac{\partial f_2}{\partial u_1} & \dots & \frac{\partial f_2}{\partial u_m} \\ \vdots & \vdots & \vdots \\ \frac{\partial f_m}{\partial u_1} & \dots & \frac{\partial f_m}{\partial u_m} \end{bmatrix} $$

has m real eigenvalues corresponding to m linearly independent right eigenvectors, then the conservation law is said to be hyperbolic.

Let us consider the case of scalar non-linear conservation laws where the flux $f$ depends only on $u$:

$$ \frac{\partial u}{\partial t} + \frac{\partial}{\partial x} f(u) = 0, \quad u(x, 0) = u_0(x) \qquad (1) $$

This can be written in integral form on the interval $[x_L, x_R]$ as

$$ \frac{d}{dt} \int_{x_L}^{x_R} u(x, t) dx = f(u(x_L, t)) - f(u(x_R, t)) $$
When a shock occurs, we need to take the solution in the space of distributions rather than the space of differentiable functions. In the case of a solution that is discontinuous along a line $s(t)$ and continuous elsewhere in the space-time domain, equation (1) can be rewritten as

$$ 
\begin{aligned}
f(u(x_L, t)) - f(u(x_R, t)) &= \frac{d}{dt} \int_{x_L}^{s(t)} u(x, t)dx + \frac{d}{dt} \int_{s(t)}^{x_R} u(x, t)dx \\
&= [u(s_L, t) - u(s_R, t)] \frac{ds}{dt} \\
&\quad + \int_{x_L}^{s(t)} \frac{\partial}{\partial t} u(x, t)dx + \int_{s(t)}^{x_R} \frac{\partial}{\partial t} u(x, t)dx
\end{aligned}
$$

where $x_L < s(t) < x_R$ are fixed points on the x-axis and $u(s_L, t)$ and $u(s_R, t)$ are the limits of $u(s(t), t)$ as $x$ tends to $s(t)$ from the left and the right respectively. The speed of the discontinuity is $S := \frac{ds}{dt}$. As $\partial_t u(x, t)$ is bounded, approaching $s(t)$ from left and right gives the following condition, called the Rankine-Hugoniot (RH) condition:

$$ S[u(s_R, t) - u(s_L, t)] = f(u(x_R, t)) - f(u(x_L, t)) \quad (2) $$

In this study, we are specifically interested in entropic shocks, i.e. shocks that satisfy the Lax entropy condition:

$$ f'(u(x_L, t)) > S > f'(u(x_R, t)) \quad (3) $$

This condition expresses the fact that characteristics on the left move at a faster speed than characteristics on the right, which results in characteristics crossing. If, instead, $f'(u(x_L, t)) < f'(u(x_R, t))$, the solution is called a rarefaction shock.

Finally, if the spatial domain is bounded, boundary conditions (BC) need to be imposed in addition to the initial condition (IC). In our study we impose boundary conditions on the spatial boundaries where characteristics are incoming.

# 2 Physics-informed neural networks

## 2.1 Background

Neural networks consist of a chain of functions that takes an input and applies successive transformations to model the function underlying the input. The most elementary form of neural network is the feed-forward neural network or multilayer perceptron (MLP), where the input goes through sequentially arranged layers of neurons. Neurons in the same layer do not communicate with each other and information only flows forward, not backward, as depicted in Figure (1). For a network with $D$ hidden layers, each layer represents an activation function $\sigma_i$, which is applied component-wise to an affine transformation of the output of the previous layer:

$$ x_i = \sigma_i(W_i \cdot x_{i-1} + b_i) \quad \forall i = 1, \dots, D $$

where $n_i \in \mathbb{N}$ for $i = 0, \dots, D+1$, $x_0 \in \mathbb{R}^{n_0}$ is the input, $W_i \in \mathbb{R}^{n_i \times n_{i-1}}$ are the weights, and $b_i \in \mathbb{R}^{n_i}$ are the biases (with the bias of the output layer $b_{D+1}$ usually set to 0). The parameters of the neural network are the set $\theta = (W_i, b_i)_{1 \le i \le D}$.

The training process then consists in optimising the parameters by applying one of the algorithms described above for a number of epochs until a stop criterion is met. For data-based
<figure><img src="image_5.png" /><figcaption>Figure 1: Feedforward neural network with 1 input layer containing 2 neurons, 2 hidden layers containing 5 neurons each, and one output layer with a single neuron.</figcaption></figure>

optimisation, the objective function for the optimisation problem $min_{\theta \in \Theta} \mathcal{L}(\theta)$, where $\Theta$ is the parameter space, can be written as

$$ \mathcal{L}(\theta) = \|u(x) - u_{NN}(x)\| \approx \sum_{j=1}^{N} |u(x_j) - u_{NN}(x_j)| $$

where $u$ is the exact solution, $u_{NN}$ is the neural network solution, $x = \{x_j\}_{1 \le j \le N}$ is a set of $N$ collocation points, and $\| \cdot \|$ is usually the $L^2$ norm (or $l^2$ in the discrete space of the collocation points).

The innovation of physics-informed neural networks (PINNs) [24] is that the loss is written to reflect physics laws in terms of partial differential equations. Consider a neural network $u_{NN}$ trained to approximate the solution of the boundary-value problem

$$ \left\{ \begin{array}{ll} \partial_t u(x,t) + \partial_x f((u(x,t)) = 0 & \forall (x,t) \in \Omega \times ]0, T] \\ u(x,0) = u_0(x) & \forall x \in \Omega \\ u(x,t)|_{\Gamma} = g(x,t) & \forall (x,t) \in \Gamma \end{array} \right. \qquad (4) $$

where $f$ is a flux function and $\Gamma$ is the boundary or a subset of the boundary of the spatial domain $\Omega \subset \mathbb{R}^d, d \in \mathbb{N}$. The total loss can thus be expressed as the sum of the PDE residual ($\mathcal{L}_{PDE}$), the loss on the initial condition($\mathcal{L}_{IC}$), and the loss on the boundary conditions ($\mathcal{L}_{BC}$):

$$ \mathcal{L}(\theta) = \lambda_1 \mathcal{L}_{PDE}(\theta) + \lambda_2 \mathcal{L}_{IC}(\theta) + \lambda_3 \mathcal{L}_{BC}(\theta) \qquad (5) $$

where $\{\lambda_i\}_{i=1,2,3}$ are positive scalar weights and

$$ \mathcal{L}_{PDE}(\theta) = \sum_{(x_j, t_j) \in \Omega \times ]0, T]} |\partial_t u_{NN}(x_j, t_j) + \partial_x f((u_{NN}(x_j, t_j)))|^2 \qquad (6) $$

$$ \mathcal{L}_{IC}(\theta) = \sum_{x_j \in \Omega} |u_0(x_j) - u_{NN}(x_j, 0)|^2 \qquad (7) $$

$$ \mathcal{L}_{BC}(\theta) = \sum_{(x_j, t_j) \in \Gamma} |g(x_j, t_j) - u_{NN}(x_j, t_j)|^2 \qquad (8) $$
The only exact data needed are the initial condition and the boundary conditions that make the PDE well-defined in the domain of interest. A term can be added to take into account any data one may have on the interior of the domain, such as experimental data, data sampled by captors, etc.

## 2.2 Neural networks and discontinuities

### 2.2.1 Universal Approximation Theorems

Universal Approximation Theorems show that, under certain conditions, certain types of functions can be approximated to any degree of accuracy by neural networks with certain properties (depth, width, and/or activation function).

One of these is Cybenko's Universal Approximation theorem (1989) [6], which shows that any bounded continuous function can be approximated by a neural network with one hidden layer and a finite number of neurons, and a sigmoid as activation function.

**Theorem 1 (Cybenko Theorem)** Let $\sigma$ be any sigmoidal function, i.e. a continuously differentiable function such that $\sigma' > 0$ and

$$ \sigma(t) \to \begin{cases} 1 & \text{as } t \to +\infty \\ 0 & \text{as } t \to -\infty \end{cases} $$

Let $x \in \mathbb{R}^n$, $y_j \in \mathbb{R}^n \forall j \in 1, \dots, N$ and $\alpha_j, b_j \in \mathbb{R}$. Then sums of the form

$$ G(x) = \sum_{j=1}^{N} \alpha_j \sigma(y_j^\top x + b_j) $$

are dense in $C([0, 1]^n)$ with respect to the supremum norm.

The proof can be found in [6, 7]. The function $G$ describes a neural network with a single hidden layer with $N$ neurons.

That same year, Hornik, Stinchcombe, and White showed that multilayer feedforward neural networks are universal approximators: in particular, they can approximate any continuous function uniformly on any compact set [9].

### 2.2.2 PINNs and discontinuities

In the context of hyperbolic laws, some problems result in discontinuous solutions, even when initial conditions are continuous. Modifications are thus necessary for neural network-based approaches to be able to simulate solutions to these problems. Some approaches model shocks with neural network solutions that feature steep gradients in the region of the discontinuity. Such approaches can be successful in the sense that the error of the solution is low, although they create diffusion around the shock. One such approach is taken by Mao, Jargtap, and Karniadakis [18] in their 2020 paper on PINNs for high-speed flows. They improve the performance of the PINN model by using collocation points clustered around the discontinuity rather than points randomly selected across the domain. The approach yields better results for the Euler equations than random sampling but the output smooths out discontinuities. The clustering also relies on a priori knowledge of the position of the shock.

This inspired other approaches mixing domain decomposition methods with PINNs, such as the conservative PINN (cPINN) method [10]. The cPINN method consists in splitting the domain
into small subdomains, training separate networks on each subdomain. Interface conditions enforce flux continuity conditions (which make it possible to “stitch” the solution together) and an average of the solution along the boundary of the local neural networks sharing this boundary, which accelerates convergence. By using this domain decomposition method, the error around the discontinuity does not propagate to the whole domain.

Lorin and Novruzi [16] suggest a non-diffusive method consisting of approximating discontinuity lines and the subdomains they define with separately defined but jointly trained neural networks. This is also a form of domain decomposition, but here it is guided by the number of entropic shock waves. The consistency of the global solution is ensured by the enforcement of the Rankine-Hugoniot condition and the discontinuity lines act as boundary conditions for the domains on either side of them. This is the method that we focus on thereafter.

A few more methods are worth mentioning. McLenny and Braga-Neto [19] develop a method using self-adaptive weights applied to each training point. These weights are trained simultaneously with the neural network parameters: the algorithm seeks to maximise the loss function with regards to the self-adaptive weights and to minimise it with regards to the network parameters. In the regions that are more difficult to learn (e.g. due to discontinuities), training points are given larger weights, forcing the minimisation algorithm to focus on these regions. Rodriguez-Torrado et al. [26] propose attention-based PINNs, which combine recurrent neural networks and attention mechanisms. Li et al. [15] address shallow water equations using a long-short term memory (LSTM) cell instead of a feed-forward neural network and enforcing $L^1$ or $L^2$ regularisation terms.

### 2.2.3 The Non-Diffusive Neural Network method

The method we choose to reproduce is the non-diffusive neural network (NDNN) method for hyperbolic conservation laws by Lorin and Novruzi [16]. The method aims specifically at modelling entropic shocks that arise in such laws.

Consider the following problem:

$$ \begin{cases} \partial_t u(x,t) + \partial_x f(u(x,t)) = 0 & (x,t) \in Q := \Omega \times [0,T] \\ u(x,0) = u_0(x) & x \in \Omega \end{cases} $$

where $\Omega \in \mathbb{R}^d$ ($d=1$ in our case) is a spatial domain, $T \in \mathbb{R}^+$, and $f: \mathbb{R} \to \mathbb{R}$ is convex and differentiable. If $\Omega$ is bounded, one may also impose boundary conditions for incoming characteristics and we denote $x_l, x_r$ the leftmost and rightmost coordinates, respectively, of the spatial domain.

The NDNN method consists in approximating shocks (i.e. discontinuous solutions) and the subdomains they define with separate neural networks, using only one loss function for the optimisation. For a problem with a single shock wave arising from the initial time, we concurrently train three networks:

1. A network $n$ to approximate the discontinuity line parametrised by $\gamma: [0, T] \to \mathbb{R}, t \mapsto \gamma(t)$;

2. Two networks denoted $N^-$ and $N^+$, defined respectively over the subdomains $Q^-, Q^+$ on either side of $\gamma(t)$.

The sets of weights and biases of each network are respectively denoted $\theta, \Theta^-, \Theta^+$, and we write $\Theta := (\Theta^-, \Theta^+, \theta)$. The loss function integrates terms for the PDE residual, the initial condition, the Rankine-Hugoniot condition, and, where relevant, boundary conditions. The outputs of all three networks and their partial derivatives contribute to the loss function. The
<figure><img src="image_6.png" /><figcaption>Figure 2: Mapping collocation points onto the subdomains either side of the discontinuity line</figcaption></figure>

objective function for the optimisation algorithm is $\min_{\Theta} \mathcal{L}(\Theta)$, where $\mathcal{L}(\Theta)$ is the (weighted) sum of the terms mentioned above.

The authors of [16] suggest taking data over the domain $R := (0, 1) \times (0, T)$ (for a 1D spatial domain), and using the following transformations to map the spatial data onto both sides of the domain as depicted in Figure 2:

$$T^{-}: R \to Q^{-}, \quad (x, t) \mapsto ((n(t) - a)x + a, t) \quad (9)$$

$$T^{+}: R \to Q^{+}, \quad (x, t) \mapsto ((b - n(t))x + n(t), t) \quad (10)$$

Using this notation, we can write the loss function as the weighted sum of the PDE residual ($\mathcal{L}_{PDE}$), the loss on the initial condition($\mathcal{L}_{IC}$), and the loss on the Rankine-Hugoniot condition ($\mathcal{L}_{RH}$):

$$\mathcal{L}_{PDE}(\Theta) = \sum_{(x,t) \in R} \left( \left| (\partial_t N^{-}(\cdot) + \partial_x f(N^{-}(\cdot))) \circ T^{-}(x,t) \right|^2 + \left| (\partial_t N^{+}(\cdot) + \partial_x f(N^{+}(\cdot))) \circ T^{+}(x,t) \right|^2 \right) \quad (11)$$

$$\mathcal{L}_{RH}(\Theta) = \sum_{t \in [0, T]} \left( \begin{aligned} & |(\partial_t n(t)[N^+(n(t), t) - N^-(n(t), t)])| \\ & - [f(N^+(n(t), t)) - f(N^-(n(t), t))]^2 \end{aligned} \right) \quad (12)$$

$$\mathcal{L}_{IC}(\Theta) = \sum_{x \in [0,1]} \left( \begin{aligned} & |(N^{-}(T^{-}(x, 0)) - u_0(T^{-}(x, 0)))|^2 \\ & + |(N^{+}(T^{+}(x, 0)) - u_0(T^{+}(x, 0)))|^2 \end{aligned} \right) \quad (13)$$
<figure><img src="image_7.png" /><figcaption>Figure 3: Non-diffusive neural network architecture</figcaption></figure>

The schematic architecture of the neural network is presented in Figure 3.

The optimisation method used in [16] is the gradient descent algorithm with a learning rate (i.e. step size) of $2 \times 10^{-3}$. The experiments reported in the paper use neural networks with one or two hidden layers and the tanh activation function. The collocation points are randomly distributed over $R = (0, 1) \times (0, T)$.

### 2.2.4 A modified NDNN algorithm

We propose a modified NDNN algorithm that does not rely on the transformations $T^\pm$ and uses conventionally accepted methods to train PINNs and improve convergence. Our modifications include:

1. Explicitly including the boundary conditions in the loss function;

2. Including a term in the loss function to set the initial condition for the discontinuity line as $n(0) = 0$;

3. Using Adam as an optimisation method instead of Gradient Descent;

4. Splitting the training points alongside the discontinuity rather than mapping them onto each subdomain, as pictured in Figure 4;

5. Using mini-batch training to accelerate training: when trying to replicate the results of the first experiment, both computers available to us ran into a stack overflow error in the automatic differentiation package when using the parameters (2,500 collocation points, neural networks with two hidden layers and 20 neurons) prescribed in the paper.

Our loss function thus includes the following terms for the PDE residual, the Rankine-Hugoniot condition, the initial condition, and the boundary conditions (where $u_l$ and $u_r$ are
the values of the solution respectively on the left and right boundaries of the space-time domain, and $x_l, x_r$ the left and right extremities of the spatial domain):

$$ \mathcal{L}_{PDE}(\mathbf{\Theta}) = \sum_{(x,t) \in Q^{-}} |\partial_t N^{-}(x,t) + \partial_x f(N^{-}(x,t))|^2 + \sum_{(x,t) \in Q^{+}} |\partial_t N^{+}(x,t) + \partial_x f(N^{+}(x,t))|^2 \quad (14) $$

$$ \mathcal{L}_{RH}(\mathbf{\Theta}) = \sum_{t \in [0, T]} |(\partial_t n(t)[N^+(n(t), t) - N^-(n(t), t)] - [f(N^+(n(t), t)) - f(N^-(n(t), t))])|^2 \quad (15) $$

$$ \mathcal{L}_{IC}(\mathbf{\Theta}) = \sum_{x \in \Omega^{-}} |(N^{-}(x, 0)) - u_0(x)|^2 + \sum_{x \in \Omega^{+}} |(N^{+}(x, 0)) - u_0(x)|^2 + |n(0)|^2 \quad (16) $$

$$ \mathcal{L}_{BC}(\mathbf{\Theta}) = \sum_{t \in [0, T]} |(N^{-}(x_l, t)) - u_l(t)|^2 + |(N^{+}(x_r, t)) - u_r(t)|^2 \quad (17) $$

<figure><img src="image_8.png" /><figcaption>Figure 4: Splitting the set of training points along the discontinuity line</figcaption></figure>

# 3 Results

We focus on reproducing experiments 1 and 4 from [16]. Expressivity tests were run on a laptop with a 12th Gen Intel(R) Core(R) processor running at 1.10 GHz with 16GB RAM, running Windows 11. All other tests were run on a desktop computer with an Intel(R) Xeon(R) processor (3.70GHz, 32GB RAM) running Linux Fedora 34. Tests were run with single-precision (Float32) data.
## 3.1 1D Burgers equation with continuous solution

To set up the necessary tools in Julia, we start by implementing a standard PINN method to solve the initial value problem (IVP) for the 1D inviscid Burgers equation on a space-time domain where its solution is continuous. Taking the convex flux function $f(u) = u^2/2$, the problem writes:

$$ egin{cases} 
\partial_t u + \partial_x f(u) = 0 & \forall (x, t) \in [0, 1] \times (0, 1/2\pi) \\
u(x, 0) = u_0(x) = \sin 2\pi x & \forall x \in [0, 1] \\
u(0, t) = 0, u(1, t) = 1 & \forall t \in (0, 1/2\pi)
\end{cases} 
\quad (18) $$

We use expressivity tests to tune the hyperparameters of the network: such tests feed data sampled from the exact solution space to the neural network, and the loss function only includes terms evaluating the distance between the neural network output and the exact solution at these collocation points. These tests allow to assess whether the neural network architecture is theoretically capable of approximating the solution. Details can be found in Annex 2.

Besides tuning network hyperparameters, we also establish that the best results are achieved by giving more weight to the initial condition term in the loss: of all the convex combinations tested, the best result is achieved with

$$ \mathcal{L}(\Theta) = 0.1\mathcal{L}_{PDE}(\Theta) + 0.8\mathcal{L}_{IC}(\Theta) + 0.1\mathcal{L}_{BC}(\Theta) $$

This is consistent with published findings on phase field equations [30] which suggest putting more weight on the initial conditions. The same paper also suggests mini-batching as a strategy to improve convergence, which we later employ.

Finally, we look at different optimisation strategies, combining an initial run of Adam optimisation followed by BFGS or L-BFGS optimisation and varying the step size of the Adam method and the number of epochs for each method. These methods are more commonly used for training neural networks than the classic Gradient Descent method. Note that “epochs” refer the number of passes over the whole dataset, while “iterations” are the total number of optimisation loops: if one trains a model on a dataset of N collocation points with mini-batches of size n, then #iterations = ⌊N/n⌋ × #epochs.

The summary of our findings from the expressivity tests is presented in Table 1, and the plotted results can be found in Figure 5.

## 3.2 Riemann problem for the Burgers equation

We first test our implementation of the original NDNN algorithm on a Riemann problem for the inviscid Burgers equation $\partial_t u + \partial_x(u^2/2) = 0$ on $\Omega = (-1, 1) \times (0, 1)$ with the initial condition

$$ u(x, 0) = u_0(x) = \begin{cases} u_L := 1 & \forall x < 0 \\ u_R := 0 & \forall x > 0 \end{cases} \quad (19) $$

As shown in Figure 21, the characteristics originating from the left-hand side of the domain cross the characteristics on the right-hand side. This creates a shock wave travelling at constant speed.

The exact solution is constant on either side of the discontinuity line:

$$ u(x, t) = \begin{cases} 1 & \text{if } x < \frac{t}{2} \\ 0 & \text{if } x > \frac{t}{2} \end{cases} \quad (20) $$
<figure><img src="image_9.png" /><figcaption>Figure 5: Model performance for equation (18), before and after tuning. Note the differences in scale between the naive model (left column) and the tuned model (right column)</figcaption></figure>
<table><caption>Table 1: Summary of parameters used in the expressivity tests. Note that the training for the naive model hits a stop criterion after 3 BFGS iterations, hence the short runtime.</caption><thead><tr><th></th><th>Naive model</th><th>Tuned model</th></tr></thead><tbody><tr><td>Hidden layer</td><td>2</td><td>3</td></tr><tr><td>Neurons per layer</td><td>10</td><td>15</td></tr><tr><td>Activation</td><td>tanh</td><td>tanh</td></tr><tr><td>Loss weights (PDE, IC, BC)</td><td>(1,1,1)</td><td>(0.1, 0.8, 0.1)</td></tr><tr><td>Adam step size</td><td>10<sup>-2</sup></td><td>5 × 10<sup>-2</sup></td></tr><tr><td>Adam iterations</td><td>100</td><td>500</td></tr><tr><td>(L-)BFGS</td><td>BFGS</td><td>BFGS</td></tr><tr><td>(L-BFGS) iterations</td><td>1000</td><td>1000</td></tr><tr><td>MSE</td><td>2.4848e-1</td><td>6.4164e-5</td></tr><tr><td>Training loss</td><td>2.8060e-1</td><td>7.7623e-4</td></tr><tr><td>Test loss</td><td>2.9987e-1</td><td>7.3151e-2</td></tr><tr><td>Runtime (sec)</td><td>9</td><td>1158</td></tr></tbody></table>

From the Rankine-Hugoniot condition, we obtain the speed $S$ of the shock and we can verify that the equation fulfils the Lax entropy condition $u_L > S > u_R$:

$$S = \frac{f(u_R) - f(u_L)}{u_R - u_L} = \frac{1}{2} \frac{u_R^2 - u_L^2}{u_R - u_L} = \frac{1}{2}(u_L + u_R)$$

Testing the algorithm on this simpler problem already revealed challenges and led to the adjustments we make to the original algorithm, in order to help steer the network towards the correct solution:

1. We explicitly specify the boundary conditions in the loss;

2. We add an initial loss term for the discontinuity line, which is assumed to start at (0, 0);

3. We modify the weights in the loss to give more importance to the initial condition;

4. We use Adam instead of GD. Neural network training often employs Adam followed by (L-)BFGS to fine-tune the results. However, (L-)BFGS is a full-batch method and the computers available to us were not always able to handle the full problem parameters, so we omit this fine-tuning step.

As the exact solution is constant on either side of the shock, we take shallow neural networks with the tanh activation function for $N^±$. For the network $n$ approximating the discontinuity,

<figure><img src="image_10.png" /><figcaption>Figure 6: Shock wave arising from the Riemann problem for the inviscid Burgers equation</figcaption></figure>
we take a network with one hidden layer of five neurons. The weights of the PDE residual, Rankine-Hugoniot, boundary, and initial conditions are set to (1, 1, 1, 10). The training dataset comprises 100 collocation points generated with a Latin Hypercube sampling strategy, and we train the model with 5,000 Adam iterations and an initial learning rate of $2 \times 10^{-3}$. The relatively high number of iterations is required to achieve a good fit of the discontinuity line. The results are plotted in Figure 7.

<figure><img src="image_11.png" /><figcaption>Figure 7: NDNN solution for the inviscid Burgers equation with initial condition (19)</figcaption></figure>

## 3.3 One shock wave, one rarefaction wave

Lorin and Novruzi's first experiment is a proof of concept of the method and seeks to approximate a solution with one entropic shock wave and one rarefaction wave. The problem is stated as:

$$ \partial_t u(x, t) + \partial_x f(u(x, t)) = 0 \quad \forall (x, t) \in \Omega \times [0, T] \qquad (21) $$

with $f(u) = 4u(2-u)$, $\Omega = [-4, 3]$, $T = 0.75$, and the initial condition

$$ u_0(x) = \begin{cases} 1 & x < -2 \\ 0.5 & -2 < x < 0 \\ 1.5 & x > 0 \end{cases} \qquad (22) $$

The exact solution is shown in Figure 8. A rarefaction wave issues from the discontinuity in the initial condition at $x = -2$ but this is not handled by the NDNN method, which will introduce diffusion (i.e. a regularised approximation) in the region of this discontinuity.
<figure><img src="image_12.png" /><figcaption>Figure 8: Exact solution of Experiment 1</figcaption></figure>

Having established that the algorithm works on a simple problem, we run expressivity tests to confirm that the network architecture suggested in the paper (networks with two hidden layers and 20 neurons per layer) is appropriate. We train one network for the discontinuity line and two networks for each side of the domain with data from the exact solution. We test different network architectures and step sizes for the Adam method, which is more tolerant to variations in the step size. By contrast, the GD method can return NaN loss values if the step size was not carefully selected, hinting at an explosion of the gradient. The details of the tests are available in Annex 2. They confirm that the network size used in [16] is broadly appropriate for the problem at hand.

We implement the NDNN algorithm with the parameters used in [16], with the only difference that we use a mini-batch procedure to accelerate calculations. We study the following methods:

1. A direct PINN method, i.e. a single neural network and a loss function embedding the PDE residual, the initial conditions, and the boundary conditions;

2. The original NDNN algorithm without explicitly stating the boundary conditions (hereafter referred to as the “naive NDNN method”);

3. The original NDNN algorithm with explicit boundary conditions and an additional initial condition on the discontinuity line;

4. The modified NDNN algorithm presented above.

We test these methods with Gradient Descent and Adam optimisers, except for the direct PINN method which we only evaluate with an Adam optimiser. We also test them on a shifted domain starting at $t = 0.1$ instead of $t = 0$: this gives a continuous initial condition on the left-hand side of the domain, which may be easier for the neural network $N^-$ to approximate. The loss function is expressed as $\mathcal{L}_{\Theta} = 0.3\mathcal{L}_{PDE}(\Theta) + 0.3\mathcal{L}_{RH}(\Theta) + 0.3\mathcal{L}_{BC}(\Theta) + \mathcal{L}_{IC}(\Theta)$.

As expected, the direct PINN solution in Figure 9 smooths out the discontinuity and does not manage to approach the constant solution in $Q^+$ nor the rarefaction wave in $Q^-$ convincingly. Similarly, the naive NDNN approach does not approximate the exact solution well, as boundary conditions and an initial condition for the discontinuity line are missing to guarantee a unique
<figure><img src="image_13.png" /><figcaption>Figure 9: Direct PINN solution for problem (21) with initial condition (22): single neural network with two hidden layers, 20 neurons per layer, and tanh activation. We train the model with 2,000 Adam iterations on mini-batches of 50 points with a learning rate of $10^{-4}$. The weights on the PDE residual, the boundary conditions, and the initial condition are (0.1, 0.1, 0.8). The training set contains 2,500 points inside the domain and 50 on each boundary.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10: Original NDNN algorithm: Solutions (top row) and errors (bottom row) with the GD (left column) and Adam (right column) optimisers for problem (21) with initial condition (22)</figcaption></figure>

solution. The tests on the shifted domain are not consistently better than those on the full domain. Even when they were better, the improvement remains marginal. The solution and
<figure><img src="image_15.png" /><figcaption>Figure 11: Modified NDNN algorithm: Solutions (top row) and errors (bottom row) with the GD (left column) and Adam (right column) optimisers for problem (21) with initial condition (22)</figcaption></figure>

error plots for the naive NDNN algorithm and the tests on the shifted domain can be found in the supplementary figures in Annex 3 (Figures 25-29).

The original and the modified methods exhibit similar patterns, as seen in Figures 10 and 11: they both manage to approach the exact solution closely when Adam is used as the optimiser. However, the error is more pronounced when using Gradient Descent, as can be clearly seen when taking time slices (Figure 12): the GD-based models fail to approximate the discontinuity line at later times, and they produce a very smooth function on the left hand-side domain. The time slices also show that even the Adam method produces diffusion around the rarefaction wave, where the exact solution is continuous but not differentiable. Notice that at later times, oscillations around the exact solution start to appear.

Table 2 presents a summary of the tests, their accuracy measured as the mean standard error $(\frac{1}{N}\sum_{i=1}^{N}|u_{NN}(x_i, t_i) - u(x_i, t_i)|^2)$, and the training runtime. As observed on the plots, the naive NDNN solution and the direct PINN both perform poorly. The original and modified versions of the NN method on the full domain get closer to the exact solution with Adam than with Gradient Descent. Tests on the shifted domain are slightly better than those on the full domain for the modified NDNN method, but much worse for the NDNN method. Finally, one noticeable advantage of the modified NDNN method, especially when using Adam, is its shorter runtime compared to the original algorithm. This is expected, as it uses half as many collocation points as the original NDNN method, but this also shows that comparable accuracy can be achieved twice as quickly with half as many points.

The table also includes tests with different network architectures. A smaller network runs faster, but does not achieve the same accuracy, and a larger network achieves good accuracy but with longer runtimes.
<figure><img src="image_16.png" /><figcaption>Figure 12: Solutions at given times for problem (21) with initial condition (22)</figcaption></figure>
<table><caption>Table 2: Results on replicating experiment 1. Training consists of 2,000 epochs of the chosen optimiser on mini-batches of 50 points. The runtime is obtained with the @timed macro in Julia.</caption><thead><tr><th>Test</th><th>Optimiser</th><th>Domain</th><th>MSE</th><th>Runtime (s)</th></tr></thead><tbody><tr><td colspan="5">2 hidden layers, 20 neurons per layer</td></tr><tr><td>Direct PINN</td><td>Adam</td><td>full</td><td>2.0487e-2</td><td></td></tr><tr><td>Naive NDNN</td><td>GD</td><td>full</td><td>1.6399e-2</td><td>8328</td></tr><tr><td>Naive NDNN</td><td>GD</td><td>shifted</td><td>2.1924e-2</td><td>8579</td></tr><tr><td>Naive NDNN</td><td>Adam</td><td>full</td><td>9.7956e-3</td><td>3207</td></tr><tr><td>Naive NDNN</td><td>Adam</td><td>shifted</td><td>4.5270e-2</td><td>3252</td></tr><tr><td>NDNN</td><td>GD</td><td>full</td><td>5.3338e-3</td><td>8500</td></tr><tr><td>NDNN</td><td>GD</td><td>shifted</td><td>1.8833e-2</td><td>8350</td></tr><tr><td>NDNN</td><td>Adam</td><td>full</td><td>4.3554e-4</td><td>3413</td></tr><tr><td>NDNN</td><td>Adam</td><td>shifted</td><td>1.6878e-1</td><td>3247</td></tr><tr><td>Modified NDNN</td><td>GD</td><td>full</td><td>5.5993e-3</td><td>7476</td></tr><tr><td>Modified NDNN</td><td>GD</td><td>shifted</td><td>4.3440e-3</td><td>6227</td></tr><tr><td>Modified NDNN</td><td>Adam</td><td>full</td><td>3.2802e-4</td><td>1645</td></tr><tr><td>Modified NDNN</td><td>Adam</td><td>shifted</td><td>6.1978e-5</td><td>1166</td></tr><tr><td colspan="5">1 hidden layer, 20 neurons per layer</td></tr><tr><td>NDNN</td><td>Adam</td><td>full</td><td>3.635494e-3</td><td>1395</td></tr><tr><td>Modified NDNN</td><td>Adam</td><td>full</td><td>3.0816828e-3</td><td>1056</td></tr><tr><td colspan="5">3 hidden layers, 10 neurons per layer</td></tr><tr><td>NDNN</td><td>Adam</td><td>full</td><td>2.0306492e-4</td><td>3642</td></tr><tr><td>Modified NDNN</td><td>Adam</td><td>full</td><td>2.843106e-4</td><td>2796</td></tr></tbody></table>
## 3.4 Shock wave generation

We test the original and the modified NDNN methods on a simplified version of experiment 4 of [16], which looks at a situation with one entropic shock wave at $t = 0$ and a second entropic shock wave appearing at $t = 0.6$ in the left hand-side domain. The problem is the inviscid Burgers equation in $\Omega \times (0, T) = (-1, 2) \times (0, 0.6)$ with the initial condition

$$u_0(x) = \begin{cases} \frac{3}{4} - \tanh 2x & -1 < x < 1.5 \\ -\frac{1}{2} & 1.5 < x < 2 \end{cases} \qquad (23)$$

We crop the spatial domain to $\Omega_c = (-1, 1)$ to only study the situation where a shock wave appears at an unspecified time $t^*$. Our initial spatial subdomains are $\Omega_0^- = [-1, 0.5]$, $\Omega_0^+ =$

<figure><img src="image_17.png" /><figcaption>Figure 13: Results from Experiment 4 on $\Omega = (-1, 1)$</figcaption></figure>
[0.5, 1], and we use two space-time-dependent neural networks $N^-$, $N^+$ on these subdomains, as well as one time-dependent network $n$ to approximate the shock. Each network has one hidden layer of 30 neurons, the training dataset comprises 900 points inside the domain, 100 points along the initial condition, and 50 points on each boundary. We put a weight of 2 on the initial condition term in the loss, and we train the model using Adam, with an initial learning rate of $2 \times 10^{-3}$, 2,000 epochs, and mini-batches of 50 points. Otherwise, the implementation is identical to the previous example: before the shock appears, there should be no jump along the discontinuity line, i.e. $f(N^+(n(t), t)) - f(N^-(n(t), t)) \simeq 0$, and so the output of the network $n$ approximating the discontinuity line should not be visible on the plots.

On this example, the original algorithm approximates the solution better than the modified algorithm (for reference, see Figures 7 and 8 in [17]). The NDNN algorithm puts the discontinuity roughly in the correct region, especially as $t$ nears $T = 0.6$. Meanwhile, the modified algorithm produces a discontinuity line that is almost constant over time. Figure 15e also shows that the jump across the discontinuity is near 0 at all times, which shows that the discontinuity line is placed in a region where the solution is actually smooth.

A potential explanation for this underperformance for the modified model is that, in this instance, the shock is not a line that roughly divides the domain in two equal halves; instead, it starts 3/4 of the way into the spatial domain. At each iteration, there are thus many more training points in the space-time domain left of the shock than on the right of the shock. When considering the fact that we are training with a mini-batch method with 50 points per batch, there may just be too few training points in the right-hand side of the domain to steer the neural networks towards the solution.

<figure><img src="image_18.png" /><figcaption>(a) Distribution of a minibatch of 50 data points at model initialisation without hard constraints on the initial condition</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Distribution of a minibatch of 50 data points at model initialisation with initial condition embedded in the model</figcaption></figure>

<caption>Figure 14: Impact on training point distribution of embedding initial conditions for the discontinuity lines</caption>

We repeat the tests on the full problem, taking two time-dependent neural networks $n_1, n_2$ initiated at 0.5 and 1.5 respectively to approximate each discontinuity. These two networks define three domains on which we train space-time dependent networks $N_1, N_2, N_3$. The original NDNN method is implemented as in the previous example but with 2,500 epochs. The modified method needs to be adapted: $n_1$ and $n_2$ may be initiated too close to each other and either cross each other, run very close together, or be near the edge (or even outside) of the domain, in which case we may end up having few or no training points in one of the subdomains, as illustrated in Figure 14a. To avoid this issue, we enforce the initial conditions on $n_1$ and $n_2$ as a hard constraint and we consider $\bar{n}_i = n_i(0) + t n_i(t)$ with $i \in \{1, 2\}$, $n_1(0) = 0.5$, $n_2(0) = 1.5$ (Figure 14b). We also train the modified model using 2,700 training points divided into larger
mini-batches of 300 points each and 5,000 epochs, to have the same number of iterations (45,000) as with the original NDNN method. We increase the weight on the initial condition to 10 (from 1) and the weight on the Rankine-Hugoniot condition to 1 (from 0.3).

<figure><img src="image_20.png" /><figcaption>(a) Original NDNN solution</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) Modified NDNN solution</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(c) Solutions at initial and final time</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(d) Discontinuity lines</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(e) Flux jumps along the discontinuity line</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(f) Training losses</figcaption></figure>

<caption>Figure 15: Results on Experiment 4</caption>

As the results in Figure 15 show, these changes to the modified method do not suffice to get a correct approximation of the solution. Moreover, the increase of collocation points means we lose the time advantage of the modified method: the original method runs in 2,472 seconds, while the modified method now takes 5 times longer (12,893 seconds).
# 4 Conclusion

We reproduced the NDNN method proposed by Lorin and Novruzi [16] and suggested a modified NDNN method for faster convergence with fewer collocation points. Our modified method improved the original algorithm for the first experiment of the paper, where the solution exhibits one shock wave arising at time $t = 0$. However, we failed to achieve similar results on experiment 4, where a shock wave is generated from a smooth initial condition. Other approaches could be tested to help place the discontinuity lines in the correct position. One option would be to steer the shock towards areas with steep gradients by adding the following term in the total loss:

$$ \sum_{i \in \#Q} (x_i - n(t_i))^2 |\nabla u_{NN}(x_i, t_i)|^2 $$

For each point $(x_i, t_i)$ in the training dataset, if the gradient in that point is high, then the network should look for parameters that bring $x_i$ and $n(t_i)$ closer. One could also implement the initial conditions as hard constraints for all networks (not just the ones approximating the discontinuity lines). This would place the characteristics at $t = 0$ exactly in the right place and should improve the ability of the network to place the shock in the correct place.

Our tests confirmed that the NDNN method can simulate entropic shock solutions without creating diffusion (i.e. smoothing out the solution) around the shock. Further work can look at testing our modified method on other cases studied by Lorin and Novruzi, such as the interaction of two shock waves and the introduction of a domain decomposition method. The authors have since published a multidimensional extension of the NDNN method (only the pre-print is available at the time of writing) [17].
References

[1] Weightinitializers.jl: Api reference. https://lux.csail.mit.edu/stable/api/Building_Blocks/WeightInitializers#Main-Functions.

[2] Zygote.jl: Limitations. https://fluxml.ai/Zygote.jl/stable/limitations/#Second-derivatives-1.

[3] Atilim Gunes Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: a survey. *Journal of Machine Learning Research*, 18(153):1–43, 2018.

[4] Edward Berman and Jacob Ginesin. The state of julia for scientific machine learning. 2024. arXiv:2410.10908.

[5] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. Julia: A fresh approach to numerical computing. *SIAM Review*, 59(1):65–98, 2017.

[6] G Cybenko. Approximation by superpositions of a sigmoidal function. *Math. Control Signals Systems*, 2(4):303–314, December 1989.

[7] Bruno Després. *Neural networks and numerical analysis*. De Gruyter Series in Applied and Numerical Mathematics. De Gruyter, Berlin, Germany, August 2022.

[8] M. Egmont-Petersen, D. de Ridder, and H. Handels. Image processing with neural networks—a review. *Pattern Recognition*, 35(10):2279–2301, 2002.

[9] Kurt Hornik, Maxwell Stinchcombe, and Halbert White. Multilayer feedforward networks are universal approximators. *Neural Netw.*, 2(5):359–366, January 1989.

[10] Ameya D. Jagtap, Ehsan Kharazmi, and George Em Karniadakis. Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems. *Computer Methods in Applied Mechanics and Engineering*, 365:113028, June 2020.

[11] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. 2017. arXiv:1412.6980.

[12] Elham Kiyani, Khemraj Shukla, Jorge F. Urbán, Jérôme Darbon, and George Em Karniadakis. Optimizing the optimizer for physics-informed neural networks and kolmogorov-arnold networks. 2025.

[13] I. E. Lagaris, A. Likas, and D. I. Fotiadis. Artificial Neural Networks for Solving Ordinary and Partial Differential Equations. *IEEE Transactions on Neural Networks*, 9(5):987–1000, September 1998. arXiv:physics/9705023.

[14] Kuangdai Leng and Jeyan Thiyagalingam. On the compatibility between neural networks and partial differential equations for physics-informed learning. 2023.

[15] Yanling Li, Qianxing Sun, Junfang Wei, and Chunyan Huang. An Improved PINN Algorithm for Shallow Water Equations Driven by Deep Learning. *Symmetry*, 16(10):1376, October 2024.

[16] Emmanuel Lorin and Arian Novruzi. Non-diffusive neural network method for hyperbolic conservation laws. *Journal of Computational Physics*, 513:113161, 2024.
[17] Emmanuel Lorin and Arian Novruzi. A diffusion-free multi-layer neural-network method for n-dimensional nonlinear hyperbolic systems. 07 2025. eprint.

[18] Zhiping Mao, Ameya D. Jagtap, and George Em Karniadakis. Physics-informed neural networks for high-speed flows. *Computer Methods in Applied Mechanics and Engineering*, 360:112789, 2020.

[19] Levi McClenny and Ulisses Braga-Neto. Self-Adaptive Physics-Informed Neural Networks using a Soft Attention Mechanism. *Journal of Computational Physics*, 474:111722, February 2023. arXiv:2009.04544 [cs].

[20] A S Miller, B H Blott, and T K Hames. Review of neural network applications in medical imaging and signal processing. *Med. Biol. Eng. Comput.*, 30(5):449-464, September 1992.

[21] Jorge Nocedal. Updating quasi-newton matrices with limited storage. *Math. Comput.*, 35(151):773, July 1980.

[22] Avik Pal. Lux: Explicit Parameterization of Deep Neural Networks in Julia. https://zenodo.org/records/17842754, April 2023.

[23] Guofei Pang, Lu Lu, and George Em Karniadakis. fPINNs: Fractional Physics-Informed Neural Networks. *SIAM Journal on Scientific Computing*, 41(4):A2603-A2626, 2019. https://doi.org/10.1137/18M1229845.

[24] M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational Physics*, 378:686-707, 2019.

[25] Prajit Ramachandran, Barret Zoph, and Quoc V. Le. Searching for activation functions. 2017. arXiv:1710.05941.

[26] Ruben Rodriguez-Torrado, Pablo Ruiz, Luis Cueto-Felgueroso, Michael Cerny Green, Tyler Friesen, Sebastien Matringe, and Julian Togelius. Physics-informed attention-based neural network for solving non-linear partial differential equations. May 2021. arXiv:2105.07898.

[27] Eleuterio F. Toro. *Notions on Hyperbolic Partial Differential Equations*, pages 41-85. Springer Berlin Heidelberg, 1997. in: Riemann Solvers and Numerical Methods for Fluid Dynamics, https://doi.org/10.1007/978-3-662-03490-3_2.

[28] Athanasios Voulodimos, Nikolaos Doulamis, Anastasios Doulamis, and Eftychios Protopadakis. Deep learning for computer vision: A brief review. *Computational Intelligence and Neuroscience*, 2018(1):7068349, 2018.

[29] Frames White, Miles Lubin, and Guillaume Dalle. Juliadiff: Differentiation tools in julia.

[30] Colby L. Wight and Jia Zhao. Solving Allen-Cahn and Cahn-Hilliard Equations using the Adaptive Physics Informed Neural Networks. July 2020. arXiv:2007.04542.
# 5 Annex 1: Julia for Scientific Machine Learning

Julia is a general-purpose, multi-paradigm language designed to be as easy to use as dynamic languages such as Python, while offering performance close to statically compiled languages such as C [5]. Publicly released in 2012, it includes features such as multiple dispatch (the ability to define functions allowing different argument types), just-in-time compiling (which improves performance), and package management (which largely avoid conflicts between packages by automatically taking dependencies into account and selecting package versions accordingly).

While Julia is less fully-fledged than older, more established languages, the SciML suite offers robust tools for scientific computing and machine learning that enable the implementation of PINNs.

## Deep Learning packages

The foremost packages for deep learning are Lux.jl and Flux.jl. Flux was the first package developed but is, at the time of writing, being replaced by Lux to fix issues with precision, randomness, performance, and parameter manipulation [22].

Packages exist to implement PINNs “out of the box”, such as NeuralPDE.jl, but we chose to not use these, in order to retain full control over network parameters, algorithm implementation, etc.

## Automatic Differentiation

Automatic differentiation packages are available for forward-mode and reverse-mode differentiation. ForwardDiff.jl, Zygote.jl, Enzyme.jl are the most notable ones. Other packages support numeric (FiniteDifferences.jl, FiniteDiff.jl), symbolic (Symbolics.jl) and high-order differentiation (TaylorSeries.jl) [29].

Forward and reverse modes refer to different ways of applying the chain rule for composed functions: forward differentiation goes from the innermost to the outermost function in the chain, while reverse differentiation goes from the outermost to the innermost function. Let $f: \mathbb{R}^m \to \mathbb{R}^n$ be the function we want to differentiate. Computationally, forward differentiation will perform $m$ evaluations and reverse differentiation, $n$ evaluations. Hence, if $m \ll n$, then forward differentiation is more efficient, and if $m \gg n$, reverse differentiation is preferable. In the context of neural networks, one may have networks with a large number of parameters, thus reverse mode is preferable to differentiate the model with respect to its parameters in the optimisation loop. For computing the partial derivatives in the PDE residual, in our case, forward differentiation is the better choice [3].

We tested combinations of four automatic differentiation packages to implement the neural network training algorithm, which involves the computation of two differentials: first, the differentials involved in the PDE are computed to calculate the PDE residual; then, the optimisation algorithm takes the derivative of the loss function according to its parameters. Note that it is not recommended to use the Zygote.jl package to compute the second derivation [2]. We used the following test case to check the accuracy of results:

$$ \nabla_{\theta}(\partial_{t} f(\theta, t)) $$

with, for $X, \theta \in \mathbb{R}^n, n \in \mathbb{N}, t \in \mathbb{R}$,

$$ f(\theta, t) = A(\theta)t = X^{\top}\theta t $$
The double derivation is easy to compute and verify: $\nabla_{\theta}(\partial_{t}f(\theta, t)) = X$. Runtimes for $n = 100$ are presented in Table 3 and led us to select Enzyme.jl as the most appropriate AutoDiff package to derive models according to their parameters, and ForwardDiff.jl to obtain the PDE residuals.

<table><caption>Table 3: Runtimes (μs) obtained with the @btime macro</caption><thead><tr><th></th><th colspan="3">Second derivation</th></tr><tr><th>First derivation</th><th>ForwardDiff</th><th>Enzyme Forward</th><th>Enzyme Reverse</th></tr></thead><tbody><tr><td>ForwardDiff</td><td>9.500</td><td>56.300</td><td>3.837</td></tr><tr><td>Zygote</td><td>68.900</td><td>incompatible</td><td>incompatible</td></tr><tr><td>Enzyme Forward</td><td>15.700</td><td>test fails</td><td>4.286</td></tr><tr><td>Enzyme Reverse</td><td>incompatible</td><td>147.500</td><td>8.267</td></tr></tbody></table>

## Optimisation

Several packages implement classic optimisation algorithms, which we use to optimise model parameters during training. Most of the following tests implement a combination adopted by much of the literature, consisting in applying a number of iterations of Adam optimisation followed by refinement with BFGS or its low-memory version L-BFGS [12].

Adam [11] is a first-order stochastic gradient descent algorithm, while the BFGS (Broyden-Fletcher-Goldfarb-Shann) algorithm is a second-order quasi-Newton method: it uses an approximation of the Hessian instead of calculating its inverse and thus is computationally cheaper to implement. However, it keeps a trace of previous iterations and thus can be memory-intensive. The L-BFGS method [21] avoids this issue by only keeping a set number of iterations in memory and replacing the oldest item with the newest at each iteration.

## Challenges

The main challenge we encountered in using Julia for this project was the limited online community of users at the time of writing. Understanding error messages and finding appropriate solutions took time. There is also a lack of variety in tutorials and detailed examples showing how to implement functions, use keyword arguments, etc. As a result, I spent considerable time reading through source code to understand issues, while answers are usually more readily available for languages such as Python with a longer history and larger support communities.

Another commonly cited issue is the limited adoption of Julia in industry circles, despite the enthusiasm in academia [4]. This impacts the amount of resources dedicated to developing new packages and functionalities.
# 6 Annex 2: Expressivity tests

## Expressivity tests for the inviscid Burgers equation in 1D with continuous solution

### Network architecture

We run tests on the number of hidden layers and of neurons per layer using a dataset of 1,000 collocation points, the tanh activation function, 100 iterations of Adam training with a learning rate of 0.1, and 1000 iterations of L-BFGS training. For this test and the following two (activation function and weight initialisation), we calculate the error as the mean square error on a different set of points, averaged over four training runs. Figure 16 shows that a balance is needed between number of layers and number of neurons: too few layers mean that the neural network cannot produce a function complex enough to approximate the exact solution, but too many layers and/or neurons result in overfitting.

<figure><img src="image_26.png" /><figcaption>Figure 16: Tests on neural network architecture</figcaption></figure>

### Activation function

From the previous experiment, we keep different combinations of number of layers and neurons per layer: 2 hidden layers of 20 neurons (501 parameters); 3 hidden layers of 15 neurons (541 parameters); and 4 hidden layers of 5 neurons (111 parameters). We compare the training results using four activation functions: the sigmoid $\sigma(x) = \frac{1}{1+e^{-x}}$; the hyperbolic tangent $\tanh(x)$ and a faster but slightly less accurate version, tanh_fast, available in the Lux.jl package; and the swish function $\text{swish}_\beta(x) = x\sigma(\beta x)$ [25]. We do not test the (thresholded) rectified linear unit (($T$)ReLU) activation functions, since their second-order derivatives are null and thus produce vanishing Hessians [14]. The mean error for all tests is of order $10^{-4}$ with a slightly better performance for the tanh-based activation functions than the sigmoid-based ones (Figure 17). As there is a potential loss of accuracy from using tanh_fast, we retain the tanh activation function for the following tests.
<figure><img src="image_27.png" /><figcaption>Figure 18: Tests on weight initialisation</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 19: Varying the weights on the different loss terms. $\lambda_2$ is plotted on the abscissa, $\lambda_3$ is given in the legend, and $\lambda_1$ is inferred as $\lambda_1 = 1 - \lambda_2 - \lambda_3$. The orange dot shows the test with unitary weights i.e. $\lambda_i = 1 \forall i \in \{1, 2, 3\}$.</figcaption></figure>

Adam

Next, we test the effect of varying the learning rate and the number of iterations of the Adam optimiser. We take again a network of 3 hidden layers and 15 neurons per layer and train the model using loss function (5) with the best combination of weights from the previous tests. We apply no iteration of the BFGS algorithm. The results in Figure 20a show that, in this instance,
a learning rate of $5 \times 10^{-2}$ and 500 epochs give the best accuracy. We also plot the evolution of the training loss in Figure 20b.

<figure><img src="image_29.png" /><figcaption>Figure 20: Tests on Adam learning rate and number of epochs</figcaption></figure>

## BFGS vs L-BFGS

We now compare the BFGS optimisation algorithm and its low-memory version, L-BFGS, in terms of runtime and final error. We use the best parameters from the previous tests: 3 hidden layers, 15 neurons per layer, tanh activation, and 500 Adam iterations with an initial learning
rate of $5 \times 10^{-2}$. We then run 1,000 iterations of (L-)BFGS optimisation. The results show that BFGS gives a slightly more accurate solution. In this instance, the runtimes of both algorithms are of the same order.

<figure><img src="image_30.png" /><figcaption>Figure 21: Error in absolute value with BFGS (top left) and L-BFGS (top right) algorithms, solutions at $t = 0$ and $t = T$ (bottom left), evolution of training loss during training (bottom right)</figcaption></figure>

## Hard constraints

In the final test, we test the effect of embedding initial and boundary conditions in the equation. This technique consists in taking a model $\bar{u}$ that is the sum of two terms: one term with no trainable parameters that expresses the initial and/or boundary conditions, and a second term that is a neural network $u_{NN}$ trained to minimise the PDE residual [13]. A multiplicative factor is applied to the second term to cancel it out at $t = 0$ and on the boundaries (trivially, $t * (x - x_l) * (x_r - x)$).

Similarly, a factor is applied to the first term so that $\bar{u}(x, 0) = u_0(x)$ and the weight of this term decreases as $t$ increases. For this factor, Lagaris et al. take a linearly decreasing factor of $(1 - t)$, Pang et al. [23] take a constant factor equal to 1, and we also test the monotonously decreasing factor $(1 - \tanh(2\pi^2 t))$. These different choices make the initial condition term evolve as depicted in Figure 22. In the case of the IVP (18) with the boundary conditions
$u(0, t) = 0$, $u(1, t) = 1 \forall t \in (0, 1/2\pi)$, the respective solutions are expressed as:

$$ \bar{u}(x, t) = u_0(x) (1 - t) + t x (1 - x) u_{NN}(x, t) \quad (24) $$

$$ \bar{u}(x, t) = u_0(x) + t x (1 - x) u_{NN}(x, t) \quad (25) $$

$$ \bar{u}(x, t) = u_0(x) (1 - \tanh(2\pi^2 t)) + t x (1 - x) u_{NN}(x, t) \quad (26) $$

<figure><img src="image_31.png" /><figcaption>Figure 22: Evolution of the different weights on the initial condition term</figcaption></figure>

The results are summarised in Table 4. We conclude that these methods work but in our case do not provide appreciable benefits in terms of accuracy or runtime.

<table><caption>Table 4: Results from enforcing hard constraints</caption><thead><tr><th>Method</th><th>MSE</th><th>runtime (sec)</th></tr></thead><tbody><tr><td>Lagaris 1998</td><td>5.37 × 10<sup>-5</sup></td><td>1,876</td></tr><tr><td>Pang 2019</td><td>8.41 × 10<sup>-5</sup></td><td>1,442</td></tr><tr><td>Lagaris modified</td><td>5.65 × 10<sup>-4</sup></td><td>2,590</td></tr></tbody></table>

## Expressivity tests for Experiment 1 (one shock wave, one rarefaction wave)

We run expressivity tests to confirm whether the network architecture suggested in [17] is as complex as necessary and as simple as possible to solve the problem in Experiment 1. We test networks with 1 or 2 hidden layers; 5, 10, 15 or 20 neurons per layer; and an Adam learning rate in $[10^{-1}, 5 \times 10^{-2}, 2 \times 10^{-2}, 10^{-2}, 5 \times 10^{-3}, 2 \times 10^{-3}, 10^{-3}]$. We run a set of tests on the network $n$ simulating the discontinuity line and another set of tests on the networks $N^{\pm}$ simulating the solution on either side of the (exact) discontinuity line. For each set of tests, we present the four networks with the smallest error (measured as the mean square error).

The results in Figures 23 and 24 show that two-layer networks with 20 neurons per layer offer the best performance, although networks with fewer neurons (10 or 15) can also approximate the solution well. Of the 56 tests performed, neural networks with only one hidden layer do not figure in the top 15.
<figure><img src="image_32.png" /><figcaption>Figure 23: Results of expressivity tests on the network approximating the discontinuity line in experiment 1. The mean square errors are, from top to bottom: 7.04 × 10⁻⁸, 7.76 × 10⁻⁸, 8.78 × 10⁻⁸, and 9.04 × 10⁻⁸.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>Figure 24: Results of expressivity tests for experiment 1. The mean square errors for the solutions on the left and right domains (bottom row) are, from left to right and top to bottom: 5.79×10⁻⁵, 6.10 × 10⁻⁵, 7.35 × 10⁻⁵, and 1.46 × 10⁻⁴.</figcaption></figure>
# 7 Annex 3: Supplementary figures for Experiment 1 (one shock wave, one rarefaction wave)

The figures below represent the neural network solutions and the error compared to the exact solution. The neural networks have two hidden layers with 20 neurons each and use the tanh activation function. The initial learning rate for Gradient Descent and the Adam algorithm is set to $2 \times 10^{-3}$. We employ a mini-batch training strategy, with batches of 50 points each.
<figure><img src="image_34.png" /><figcaption>Figure 25: Naive NDNN solution (top two rows) and error (bottom two rows)</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Figure 26: Original NDNN solution</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure 27: Original NDNN error</figcaption></figure>
<figure><img src="image_37.png" /><figcaption>Figure 28: Modified NDNN solution</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure 29: Modified NDNN error</figcaption></figure>
Inria

RESEARCH CENTRE
Centre Inria d'Université Côte d'Azur

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399