Learning-based solutions to nonlinear hyperbolic
PDEs: Empirical insights on generalization errors

Bilal Thonnam Thodi1,2 , Sai Venkata Ramana Ambadipudi2 , Saif Eddin Jabari1,2
1
New York University Tandon School of Engineering, Brooklyn, NY
2
New York University Abu Dhabi, Abu Dhabi, UAE{btt1, sa183, sej7}@nyu.edu

Abstract
We study learning weak solutions to nonlinear hyperbolic partial differential equations (H-PDE), which have been difficult to learn due to discontinuities in theirsolutions. We use a physics-informed variant of the Fourier Neural Operator (πFNO) to learn the weak solutions. We empirically quantify the generalization/outof-sample error of the π-FNO solver as a function of input complexity, i.e., thedistributions of initial and boundary conditions. Our testing results show that
π-FNO generalizes well to unseen initial and boundary conditions. We find thatthe generalization error grows linearly with input complexity. Further, adding
a physics-informed regularizer improved the prediction of discontinuities in thesolution. We use the Lighthill-Witham-Richards (LWR) traffic flow model as a
guiding example to illustrate the results.

1

Introduction

Hyperbolic partial differential equations (H-PDEs) arise in the study of nonlinear wave motion inapplications such as the behavior of water waves, vehicular traffic flow, and even in high-speed
physics like blast waves and sonic booms. Typically, (inviscid) H-PDEs are of the form ut + cux = 0,
where u is a conserved quantity. Wave motion in H-PDEs is characterized by finite propagation speedsc. The solutions are characteristic trajectories emanating from the datum (e.g., the initial conditions).
Nonlinear H-PDEs, where the nonlinearity results in different wave propagation speeds, e.g., c = c(u),result in characteristic lines that may cross somewhere in the domain, and at those points, the solution
is multi-valued, which means that the problem does not have a solution in the classical sense. Here,weak solutions that permit discontinuities are introduced, which are conventionally solved using finite
volume or finite difference-based numerical schemes [1, 2].
Recently, deep learning-based (DL) methods have appeared that aim to overcome the limitations ofconventional numerical solvers, namely, high computational cost, grid dependence, and knowledgeof complete initial and boundary conditions. The DL solvers have shown remarkable results forboth forward and inverse problems, especially where the PDE solutions are smooth [3, 4]. However,
DL solutions having high irregularities (e.g., discontinuities), such as those arising in the weaksolutions of H-PDEs, are only partially successful [5–8]. This is partly due to the ill-posedness of
H-PDEs, i.e., derivatives are not defined everywhere, and the H-PDE residuals do not form a correct
physics loss metric. Further, these studies lacked a systematic procedure to evaluate the out-of-sample
performance or the generalization error.
To this end, we explore Fourier Neural Operators (FNO) [4] for learning weak solutions of H-PDEs.
We propose a systematic training and testing experiment where the FNO solver is trained withsolutions of elementary input conditions and evaluated for solutions of general input conditions. We
quantify the empirical generalization error of the FNO solver as a function of the input complexity,
i.e., how the out-of-sample error grows as the distribution of input conditions becomes more general.
Machine Learning and the Physical Sciences workshop, NeurIPS 2022.

To capture physically consistent weak solutions (e.g., shocks), we propose using an integral form of
the H-PDE (written in discrete form) as the physics loss function instead of the H-PDE residual. We
demonstrate these results using vehicular traffic flow as an example.

2

Methods

Problem setting We consider the LWR traffic flow model [9, 10] as the guiding example fornonlinear scalar H-PDEs. The LWR model is a continuum description of the flow of vehicles on aroad. Consider a space-time domain Ω ⊂ R × R+ . Denote by u(x, t) : Ω → [0, umax ] the density of
traffic at position x ∈ R and time t ∈ R+ . Here density refers to the average number of vehicles perunit length. Let q(x, t) : Ω → [0, qmax ] the traffic flux, which is the number of vehicles crossing apoint per unit time. The LWR model describes the evolution of traffic density based on the principleof vehicular conservation:
ut + qx = 0; u(x, 0) = ū0 ; u(xb , t) = ūb ; (x, t) ∈ Ω;

(1)

where ut ≡ ∂u∂t , ū0 is an initial condition, and ūb is a boundary condition. For the LWR model, one
typically prescribes a flux function q = f (u), where f (u) : [0, umax ] → [0, qmax ], which is concavein the context of traffic flow. The LWR model is given as
ut + f ′ (u)ux = 0; u(x, 0) = ū0 ; u(xb , t) = ūb ; (x, t) ∈ Ω;

(2)

We refer to (2) as the forward problem if ūb corresponds to solutions at the boundary points (xb , t) ∈∂Ω and as the inverse problem if ūb is replaced with densities at random points in the domain Ω i.e.,(xb , t) ∈ Ω. We use f (u) = uvmax (1 − u/umax ), where umax is the maximum traffic density andvmax is the maximum traffic speed. A major difficulty in solving the forward problem is handling
discontinuities when they appear in the solution due to nonlinearity in the system. The inverse
problem poses the additional challenge that the boundary condition is unknown, and the solution
needs to be inferred from partially observed measurements. We tackle both challenges using a deeplearning-based solver, discussed below.
Fourier Neural Operator solver We are interested in learning the solution operator that maps the
input function a := (ū0 , ūb ) to the weak solution u(x, t) over Ω. Let the input function be a ∈ A
and output function be u ∈ U. The problem (2) can be rephrased as one of learning an operator
G : A → U. We approximate G using the Fourier Neural Operator (FNO) of Li et al. [4], represented
by the parametric model GΘ , and given byub = GΘ (a) = Q ◦ F (L) ◦ F (L−1) ◦ · · · ◦ F (2) ◦ F (1) ◦ P (a)(3)
where {F (l) }L
l=1 is a set of Fourier operators while P and Q are projection operators. A single
Fourier operator is defined asF(z) = σ W · z + IFFT R · FFT(z)(4)
for any latent input z, where FFT and IFFT denote the Fourier transform and its inverse. Θ ={W (l) , R(l) }L
l=1 is the set of trainable parameters of the FNO operator GΘ . Our motivation for using
the FNO operator (3) is its efficient approximation in the Fourier domain and it ability to learncomplex dynamics [4]. The parameter complexity involved in learning GΘ depends on the size of R,
which is independent of the domain size |Ω|.
Physics-informed training The FNO model (3) can be trained end-to-end in a supervised learning
framework over an appropriately defined loss function. We perform physics-informed training where
the loss function has two parts − an empirical trainingP
data loss Ldata and a physics constraint loss
Lphys to emulate the PDE operator. Ldata is simply n∈N ||GΘ (a(n) ) − u(n) ||2 , where N is thenumber of samples and u(n) is the nth sample solution.
One could use the PDE residual of (2) to form the physics loss Lphys , as in the conventional physicsinformed neural networks [3]. However, (2) is not a well-posed PDE, i.e., derivatives are not defined
everywhere, especially near discontinuities and hence not a well-defined loss metric. Thus, we resort
to the integral form of (2) to form the physics loss Lphys as follows:
Lphys =u(x, t + ∆t) − u(x, t) +∆t[q(x − ∆x/2, t) − q(x + ∆x/2, t)] (x,t)
∆x
2
2

(5)

where u(x, t) and q(x, t) are defined over a discretization of Ω, and ⟨·⟩(x,t) is a concatenation operator.
Accordingly, we train two different FNO models with two different objective functions, as shown
below:(i) FNO model : min LdataΘ(ii) π-FNO model : min Ldata + λLphysΘ(6)

Data and training experiments We obtain the training and testing dataset by numerically simulating (2) for different input conditions. The numerical scheme used for generating the datasets and the
FNO training codes is described in Appendix A and Appendix C.
i2

i3

Train

i4

i5

Test

i6

i7

i8

ū0 (x)

i1

x

x

x

x

x

x

x

x

b2

b3

b4

b5

b6

b7

b8

t

t

t

t

t

t

t

t

ūb (t)b1

Figure 1: Initial and boundary conditions used for training and testing the FNO solver.
We perform systematic training and testing experiments to quantify the generalization performanceof the FNO solver. During training, the FNO solver is shown solutions of simple dynamics, forinstance, generated from step-wise initial conditions (a single vehicle queue) and one or two steppedwavelet-like boundary conditions (emulating vehicles stopping at traffic lights). The FNO solver isthen tested with solutions of complex dynamics generated from general initial conditions (multiple
vehicle queues) and multi-stepped wavelet-like boundary conditions (multiple stops at traffic lights).
The set of initial and boundary conditions used in training and testing are shown in Figure 1. Training
data consists of solutions using inputs (i0 -i3 , b0 -b2 ). Testing data are made of (i0 -i9 , b0 -b2 ) for
evaluating initial conditions, and (i0 -i3 , b0 -b8 ) for evaluating boundary conditions.
The goal is to train the FNO solver with simple solutions and assess the out-of-sample error as input
conditions become complex. For the LWR traffic flow example, the input complexity refers to (a)how the vehicles are distributed spatially at time t = 0 (i.e., ū0 ) and (b) the impact of traffic signals
at the road exit x = xmax (i.e., boundary condition ūb ). These two input factors put together can
generate complex dynamics u. Training details are summarized in Appendix B.

3

Results

Generalization error v/s input complexity The out-of-sample errors as a function of input
conditions are summarized in Figure 2. Each data point is the average mean absolute error (MAE) of
50 samples. A piece-wise linear trendline is fitted to the error plots. The trendline shows that MAE isnearly constant for input conditions seen during training (i0 -i3 , b0 -b2 ) and steadily increases for input
conditions at testing (i4 -i9 , b3 -b9 ).

4
3← Train

Test →

FNO
π-FNO

4
← Train

Test →

21
i0 i1 i2 i3 i4 i5 i6 i7 i8 i9Initial conditions(b) Inverse problem

MAE (vehs/km)

MAE (vehs/km)(a) Forward problem

b0 b1 b2 b3 b4 b5 b6 b7 b8
Boundary conditions

3← Train

Test →FNO
π-FNO← Train

Test →

2i0 i1 i2 i3 i4 i5 i6 i7 i8 i9Initial conditions

b0 b1 b2 b3 b4 b5 b6 b7 b8
Boundary conditions

Figure 2: Empirical generalization (out-of-sample) error as function of input conditions.
3

We define the generalization error as the slope of thistrendline, summarized in Table 1. For instance, Figure 2
says that MAE increases by 0.141 vehs/km (+0.12%) for
every additional traffic light at the road exit ūb using πFNO model. Similarly, an additional non-uniformity inthe initial condition ū0 increases MAE by 0.142 vehs/km(+0.12%).

Table 1: Error rates (vehs/km)
π-FNO

FNO
ForwardInverse

ū0

ūb

ū0

ūb

0.222
0.210

0.279
0.211

0.142
0.085

0.141
0.194

Figure 2 concludes that the generalization error for the testing set grows linearly with the input
complexity for both the forward and inverse problems. Also, π-FNO incurs lower error rates
compared to the FNO model.
Sample predictions The predicted and true solutions for four different input conditions are shownin Figure 3. The π-FNO solver is only shown density dynamics similar to Figure 3a at the training
stage. Figure 3b-3d are the density dynamics that π-FNO solver generalizes flawlessly. This implies
that the π-FNO solver learned to capture the traffic queuing dynamics (i.e., vehicle queue formationand dissipation) as a function of the boundary flows. Also, the inverse problem results in Figure
3c, and 3d shows that π-FNO can qualitatively recover solution without knowledge of boundary
conditions but only with sparse trajectory measurements.
(a) Forward (Train)
True(d) Inverse (Test)
True

100

0.5
0.0
1.0

x (km)(c) Inverse (Test)
True

π-FNO

π-FNO

π-FNO

0

π-FNO

u (vehs/km)

x (km)

1.0(b) Forward (Test)
True

0.5
0.0

0.0

0.1

0.0

t (hrs)0.1

0.0

t (hrs)0.1

0.0

t (hrs)0.1

t (hrs)Figure 3: Comparison of true and π-FNO predicted solutions. Sub-figure (a) is a training scenario,and (b)-(d) are test scenarios. The dotted curve is the input locations of ūb for the inverse problem.

Physics-informing on the behavior of shock solutions In Figure 4, we compare solutions (zoomedin) for a constant and a step-wise initial conditions. We see that the FNO model (physics-uninformed)
produces noisy predictions, whereas the π-FNO model (physics-informed) smoothens these artifacts.
This suggests the benefits of physics-informing in producing physically consistent solutions.
u (vehs/km) u (vehs/km)

t = 2 (s)

True
FNO

90
80

True
π-FNO

90
80
400

500

x (m)

400

500

x (m)

400

500

t = 0 (s)u (vehs/km) u (vehs/km)(a) Example 1
t = 1 (s)t = 0 (s)x (m)(b) Example 2
t = 1 (s)t = 2 (s)

80
40

True
FNO

80
40

True
π-FNO

400

500

x (m)

400

500

x (m)

400

500

x (m)

Figure 4: Solution profiles as a function of x

4

Summary and Discussion

We explored the Fourier Neural Operator (FNO) in learning the weak solutions of scalar non-linear
hyperbolic partial differential equations (H-PDEs), which has seen limited success in the deep4

learning-based computational literature. We focused on quantifying the generalization error of the
FNO solver as a function of input complexity, taking vehicular traffic flow as an example. We foundthat the FNO solver can be trained using simple solutions and can easily generalize to complex inputs
with an acceptable error tolerance − the out-of-sample errors grew linearly with input complexity.
We also showed the benefits of physics-informing in predicting physically consistent solutions, e.g.,correct shock behavior. To the authors’ knowledge, this is the first empirical study on generalization
capabilities of learning-based solvers for non-linear H-PDEs.
A limitation of the current solver is that it requires a regular grid-like computational domain, partly
due to the Fourier Transform operator. Our efforts continue to extend these solvers to irregulargraph-like computational domains, e.g., to solve traffic flow on a city network.

5

Broader Impact

The modeling and control of dynamical systems such as road traffic, water supply, and communication
networks are fundamental in functioning large-scale urban cities, which contribute to one-third of
the global carbon footprint. The techniques developed in this study aid in building low-resource
computational tools for controlling these dynamical systems, which are often modeled as partial
differential equations. The data-driven learning paradigms studied in this work help advance thedevelopment of edge computing infrastructures such as connected vehicles, smart personal gadgets,and even augmented/virtual reality applications.

References[1] Randall J LeVeque. Numerical methods for conservation laws, volume 214. Springer, 1992.
[2] G B Whitham. Linear and nonlinear waves. John Wiley & Sons, New York, NY, 1974.
[3] Maziar Raissi, Paris Perdikaris, and George 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.
[4] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik
Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametricpartial differential equations. In International Conference on Learning Representations, 2021.
URL https://openreview.net/forum?id=c8P9NQVtmnO.
[5] Ruben Rodriguez-Torrado, Pablo Ruiz, Luis Cueto-Felgueroso, Michael Cerny Green, Tyler
Friesen, Sebastien Matringe, and Julian Togelius. Physics-informed attention-based neuralnetwork for hyperbolic partial differential equations: application to the buckley–leverett problem.
Scientific Reports, 12(1):7557, 2022. doi: 10.1038/s41598-022-11058-2. URL https://doi.
org/10.1038/s41598-022-11058-2.
[6] Xiaoping Zhang, Tao Cheng, and Lili Ju. Implicit form neural network for learning scalar
hyperbolic conservation laws. In Joan Bruna, Jan Hesthaven, and Lenka Zdeborova, editors,
Proceedings of the 2nd Mathematical and Scientific Machine Learning Conference, volume 145
of Proceedings of Machine Learning Research, pages 1082–1098. PMLR, 16–19 Aug 2022.
URL https://proceedings.mlr.press/v145/zhang22a.html.
[7] Ravi G. Patel, Indu Manickam, Nathaniel A. Trask, Mitchell A. Wood, Myoungkyu Lee, Ignacio
Tomas, and Eric C. Cyr. Thermodynamically consistent physics-informed neural networks for
hyperbolic systems. Journal of Computational Physics, 449:110754, 2022. ISSN 0021-9991.
doi: https://doi.org/10.1016/j.jcp.2021.110754. URL https://www.sciencedirect.com/
science/article/pii/S0021999121006495.
[8] Ameya D. Jagtap, Ehsan Kharazmi, and George Em Karniadakis. Conservative physicsinformed neural networks on discrete domains for conservation laws: Applications to forward
and inverse problems. Computer Methods in Applied Mechanics and Engineering, 365:113028,
2020. ISSN 0045-7825. doi: https://doi.org/10.1016/j.cma.2020.113028. URL https://www.
sciencedirect.com/science/article/pii/S0045782520302127.
5

[9] M. Lighthill and G. Whitham. On kinematic waves. II. A theory of traffic flow on long crowded
roads. In Royal Society of London. Series A, Mathematical and Physical Sciences, volume 229,
pages 317–345, 1955.
[10] Paul I. Richards. Shock Waves on the Highway. Operations Research, 4(1):42–51, February
1956. ISSN 0030-364X, 1526-5463.

Acknowledgments and Disclosure of Funding
This work was supported in part by the NYUAD Center for Interacting Urban Networks (CITIES),
funded by Tamkeen under the NYUAD Research Institute Award CG001, and in part by the NYUAD
Research Center on Stability, Instability, and Turbulence (SITE), funded by Tamkeen under the
NYUAD Research Institute Award CG002. The views expressed in this article are those of theauthors and do not reflect the opinions of CITIES, SITE, or their funding agencies.

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [Yes] See Section 4(c) Did you discuss any potential negative societal impacts of your work? [N/A](d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A](b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] See Appendix
A-C.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they werechosen)? [Yes] The hyperparameters and training details are discussed in Appendix B(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes] The error bars are plotted in Figure 2(d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [Yes] See Section B
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes](b) Did you mention the license of the assets? [Yes](c) Did you include any new assets either in the supplemental material or as a URL? [Yes]
See Section C(d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]6

A

Godunov scheme for LWR simulation

The LWR traffic flow model can be solved using Godunov’s numerical scheme. Let i and j be thespace and time index, denote u(i,j) as the average traffic density for each cell (i, j). The traffic
density is updated as:∆t u(i,j+1) = u(i,j) +q(i−1/2,j) − q(i+1/2,j)(7)∆xwhere q(i−1/2,j) is the cell boundary flux from cell (i − 1, j) to cell (i, j). ∆t and ∆ x denote thetemporal and spatial width. The boundary flux is given by,no(i,j−1)q(i−1/2,j) = min Qdem , Q(i,j)
sup
n
o
(i,j)
(i,j+1)q(i+1/2,j) = min Qdem , Qsup
where
(i,j)



Qdem =

f u(i,j)qmax



if u(i,j) ≤ ucr
otherwisef u(i,j)if u(i,j) > ucr
qmax
otherwise
where f (u) is the traffic flux function (also called Fundamental relation in traffic flow literature)f (u) = u(1 − u/umax )vmax .
where vmax is maximum speed, umax is maximum density and qmax is the maximum traffic flux.
The above numerical scheme is run for different sets of initial and boundary conditions shown in
Figure 1. For the inverse problem, we draw random vehicle trajectories to represent ūb . For all thetesting, we used 10 random vehicle trajectories as the input.
Q(i,j)sup =B

Hyperparameters and training details

The hyperparameters used in the study are summarized in Table 2. The training hyperparameters arechosen by an independent trial-and-error experiment. We modified the original FNO implementationfrom Li et al. [4] to implement the π-FNO model in python using the PyTorch machine learning
library; see Appendix C. The training was performed on a GPU cluster (NVIDIA Tesla V100 32GB)
and the total run time was around ∼ 45 min for a single training experiment.
Table 2: Hyper-parameters used in the study
LWR simulation

FNO model

FNO training

space dimension

1000 m# Fourier layers L

4coefficient λ2.0time dimension

600 sec# modes in x dimension

24# epochs

500

discretization size

(50 × 600)

# modes in t dimension

128batch size

128

cell width ∆x

20 m# latent width

64

learning rate

1e − 3cell width ∆t

1 sec

lifting operator P

Linear layer
with depth 128

learning rate scheduler

step-wise

max density umax

120 vehs/km

optimizer

Adam GD

lifting operator Q

2-layer FNN
with depth 128# training samples

5200# testing samples

400

max flow qmax

1800 vehs/hr

The objective function coefficient λ is obtained by training a series of π-FNO models for λ ={0, 0.05, 0.1, . . . , 0.95, 1.0}. The λ value corresponding to the least validation error is considered
optimal.

C

Datasets and codes

The datasets, codes, and pretrained models are shared at https://github.com/bilzinet/pifno
under the MIT license.
7