Physics-informed neural networks and their failure modes

Before this lesson: reverse-mode automatic differentiation (lesson 01); how a time stepper moves information forward (lesson 02); the spectral Burgers solver (lesson 08); MLPs and their training (deep learning 01, deep learning 02). Helpful: the condition number (optimization 02). Time: about 90 minutes, plus the notebook.

---

Why this matters #

• A popular method. PINNs (Raissi, Perdikaris and Karniadakis, 2019) appear in thousands of papers on fluids, heat transfer, geophysics and biology. Knowing when they work saves months of effort.
• Known failures. Plain PINNs fail on simple transport problems (Krishnapriyan et al., NeurIPS 2021). A systematic review of ML for fluid-related PDEs found that 79% (60 of 76) of the papers that claim to beat a standard numerical method compare against a weak baseline (McGreivy and Hakim, Nature Machine Intelligence, 2024).
• Physics penalties everywhere. The same issues (balancing loss terms, conditioning, causality in time) appear in every loss with a physics term, for example a mass-conservation penalty added to the loss of a weather emulator, or a residual term added to a neural operator.
• Inverse problems. With sparse, noisy data and a known PDE, the PINN loss is a regularised inverse problem, close to weak-constraint 4D-Var (lesson 06).

1. Predict first #

2. The picture #

A PINN in one formula #

Choose a network $u_\theta(x, t)$ with parameters $\theta$. Pick collocation points $(x_i, t_i)$ inside the space-time domain, points $x_j$ at $t = 0$, and points on the boundary. Minimise

where $g$ is the initial condition, $\lambda_{ic}, \lambda_{bc} > 0$ are weights, and $L_{bc}$ is a squared boundary mismatch. For advection the residual is $r_\theta = \partial_t u_\theta + \beta\,\partial_x u_\theta$. The derivatives of the network with respect to its inputs come from automatic differentiation. The network is one function of all space and time at once. Observations, if you have any, add one more squared-error term.

The simplest case: one Fourier mode #

Try $u = a(t)\sin x + b(t)\cos x$. Then

so the PDE becomes the harmonic oscillator $a' = \beta b$, $b' = -\beta a$ with $a(0) = 1$, $b(0) = 0$. Its solution $a = \cos\beta t$, $b = -\sin\beta t$ gives $u = \sin(x - \beta t)$. So the PINN must learn a function that oscillates $\beta/(2\pi)$ times in $t$ on $[0, 1]$: 0.16 times at $\beta = 1$, 4.8 times at $\beta = 30$.

A solver marches; a PINN does not #

A time stepper (lesson 02) computes $u(t + \Delta t)$ from $u(t)$. Information flows forward from the initial condition, as in the physics. A PINN fits all of space-time at once. The residual at $t = 0.9$ only links values of $u$ at nearby points, and a field that knows nothing about the initial condition (a constant, for example) makes it zero. The initial condition reaches late times only through the optimiser, step by step.

Interactive figure: loading…

What goes wrong at large $\beta$. Three things happen together.

• Imbalance: $\beta u_x$ is 30 times larger at $\beta = 30$, so the residual term is about $\beta^2$ times larger and dominates the gradient. The initial-condition term, which says which solution we want, gets a small share.
• Spectral bias: tanh networks fit smooth, low-frequency functions first, and the solution oscillates fast in $t$.
• No causality: nothing makes early times converge before late times.

The typical result: the PINN fits the initial condition near $t = 0$, but after a short time it no longer follows the travelling wave (Krishnapriyan et al., Fig. 1). The network is big enough: a better training procedure makes the same network solve the problem (section 3).

3. Definitions and main results #

The PINN problem #

Take a PDE $\partial_t u = \mathcal N[u]$ for $x$ in a domain $\Omega$ and $t \in [0, T]$, with initial condition $u(\cdot, 0) = g$ and a boundary condition. Examples: advection, $\mathcal N[u] = -\beta u_x$; viscous Burgers, $\mathcal N[u] = -u u_x + \nu u_{xx}$ with viscosity $\nu > 0$. The residual of a candidate $\tilde u$ is $r = \partial_t\tilde u - \mathcal N[\tilde u]$, and the PINN loss is the formula of section 2.

Nested automatic differentiation. The network output $u_\theta(x, t)$ is a scalar, so one reverse sweep gives both $u_x$ and $u_t$ at a point (lesson 01). The second derivative $u_{xx}$ needs a second sweep, through the computation of $u_x$. Training needs $\partial L/\partial\theta$: a third sweep through everything. So the inner sweeps must be recorded as computations themselves (create_graph=True in PyTorch), and memory grows with each level. The activation must also be smooth: a ReLU network is piecewise linear, so $u_{xx} = 0$ almost everywhere, and it cannot see a viscous term.

Theorem 1: for advection, a small residual gives a small error #

Let $u$ solve $u_t + \beta u_x = 0$ on $[0, 2\pi)$ with periodic boundary conditions. Let $\tilde u$ be smooth and periodic, with residual $r = \tilde u_t + \beta\tilde u_x$ and error $e = \tilde u - u$. With $\lVert\cdot\rVert$ the $L^2$ norm over $x \in [0, 2\pi)$,

Why this is true, in words. Advection only moves the error around; it never makes it larger. So the error can grow only through the residual, which acts as a source term, and at most by the size of that source.

Proof of Theorem 1

1. The PDE is linear and $u$ has zero residual, so $e_t + \beta e_x = r$.
2. Multiply by $e$ and integrate over $x$: $\tfrac12\tfrac{d}{dt}\lVert e\rVert^2 + \beta\int e\,e_x\,dx = \int e\,r\,dx$.
3. $e\,e_x = (e^2/2)_x$, and the integral of the derivative of a periodic function is zero. With Cauchy–Schwarz: $\tfrac12\tfrac{d}{dt}\lVert e\rVert^2 = \int e\,r\,dx \le \lVert e\rVert\,\lVert r\rVert$.
4. Where $\lVert e\rVert > 0$, the left side equals $\lVert e\rVert\,\tfrac{d}{dt}\lVert e\rVert$, so $\tfrac{d}{dt}\lVert e\rVert \le \lVert r\rVert$. (To include times where $\lVert e\rVert = 0$, apply step 3 to $\sqrt{\lVert e\rVert^2 + \delta}$ and let $\delta \to 0$.)
5. Integrate from $0$ to $t$. $\blacksquare$

Three consequences:

• For advection, a small loss does mean a small error. So when a PINN for advection is wrong, its loss is not small: the optimiser did not find the small loss that exists. The failure is an optimisation failure. Krishnapriyan et al. reach the same conclusion from experiments.
• The bound says how small "small" must be. If the initial condition is exact and $\lVert r(\cdot, s)\rVert = \rho$ at all times, the relative error at $t = 1$ is at most $\rho/\lVert u\rVert$, with $\lVert u\rVert = \sqrt\pi$ for $u = \sin(x - \beta t)$. A 1% guarantee needs $\rho \le 0.0177$, that is, a mean squared residual $L_r = \rho^2/(2\pi) \le 5 \times 10^{-5}$.
• For a nonlinear PDE the same steps give a Grönwall factor (next box). For a chaotic system errors grow like $e^{\lambda_1 t}$ with the largest Lyapunov exponent $\lambda_1$ (lesson 05), so the residual controls late times only weakly.

The nonlinear case: a Grönwall factor

If $\langle \mathcal N[\tilde u] - \mathcal N[u], \tilde u - u\rangle \le \ell\,\lVert\tilde u - u\rVert^2$ at all times, the same steps give

$$\lVert e(t)\rVert \le e^{\ell t}\lVert e(0)\rVert + \int_0^t e^{\ell(t - s)}\lVert r(s)\rVert\,ds .$$

For periodic viscous Burgers, $\ell = \tfrac12\max_{x,t}(-u_x)$ works, with $u$ the true solution: the term cubic in $e$ integrates to zero, and the viscous term only helps. $\ell$ is large near a steep front, so the same residual allows a larger error there. This is the idea behind the PINN error estimates of Mishra and Molinaro (2023).

Why training fails: three views #

(a) The difficulty is in the coordinates, not in the physics. Write $v(\xi, t) = u(\xi + \beta t, t)$, the solution seen by an observer who moves with speed $\beta$. Then $v_t = u_t + \beta u_x$. The map $(x, t) \mapsto (x - \beta t \bmod 2\pi,\ t)$ preserves area on the periodic domain, so, with means over the domain in place of the sums,

This does not depend on $\beta$: in the moving frame the solution is $v = g$, constant in $t$, for every $\beta$. So $\beta$ only changes how hard the solution is to represent and optimise in the coordinates $(x, t)$ that the network sees. On a grid with $N$ points in $x$, the operator $\partial_t + \beta\partial_x$ has a condition number of about $\beta N$, and the squared loss squares it to about $(\beta N)^2$ (Krishnapriyan et al., App. B, a rough estimate). For Burgers the characteristics depend on the unknown solution, so this change of frame is not available.

(b) Imbalance and spectral bias: the neural tangent kernel view. For a wide network trained with small steps, the outputs change almost linearly in the parameters. Let $J_u$ be the matrix with rows $\partial u_\theta(x_j, 0)/\partial\theta$ and $J_r$ the matrix with rows $\partial r_\theta(x_i, t_i)/\partial\theta$. The neural tangent kernel (NTK) blocks are $K_{uu} = J_uJ_u^\top$, $K_{rr} = J_rJ_r^\top$ and $K_{ur} = J_uJ_r^\top$. For the loss $\tfrac12\lVert r_\theta\rVert^2 + \tfrac{\lambda_{ic}}2\lVert u_\theta - g\rVert^2$, gradient flow in the training time $\tau$ gives

When the kernel stays almost constant, an error component along an eigenvector with eigenvalue $\kappa$ decays like $e^{-\kappa\tau}$ (Wang, Yu and Perdikaris, 2022, analyse PINNs this way). Two facts matter here.

• $J_r$ contains $\beta\,\partial u_x/\partial\theta$, so $K_{rr}$ grows like $\beta^2$ for large $\beta$, while $K_{uu}$ does not depend on $\beta$. The residual is driven down fast and the initial condition slowly. For a random one-hidden-layer tanh network with 256 units, $\operatorname{tr}K_{rr}/\operatorname{tr}K_{uu}$ is 0.64, 16 and 149 at $\beta = 1, 10, 30$ (code in section 5).
• The large eigenvalues of the kernel belong to smooth functions. This is spectral bias (Rahaman et al., 2019): gradient descent fits the low-frequency part of a target first. For the same kind of network, the initial fitting speed of $\cos 30t$ is about 300 times smaller than that of $\cos t$.

Together: the network quickly reaches fields with a small residual that ignore the initial condition (for example nearly flat ones), and it leaves them only slowly.

Interactive figure: loading…

(c) No causality. Wang, Sankaran and Perdikaris (2024) measured the residual separately at each time during training. A plain PINN often reduces the residual at late times as fast as, or faster than, at early times. But the correct late-time field depends on the correct early-time field, so this is the wrong order.

The fixes #

1. Curriculum training (Krishnapriyan et al., §5.1). Train at a small $\beta$, use the weights to start the training at a larger $\beta$, and continue in steps up to the target (for example $\beta = 1, 5, 10, 20, 30$). Each stage starts near a good solution of the next problem. Relative errors in their Table 1, plain → curriculum: $\beta = 20$: 0.750 → $9.84\times10^{-3}$; $\beta = 30$: 0.897 → $2.02\times10^{-2}$; $\beta = 40$: 0.961 → $5.33\times10^{-2}$. Same network, so its size was never the problem.

2. Time marching (sequence-to-sequence) (§5.2). Split $[0, T]$ into windows of length $\Delta t$. Train a PINN on the first window; its prediction at the end of the window is the initial condition of the next window. This is a time stepper whose "step" is a training run: it respects causality, but costs one training per window, and errors pass from window to window. For convection at $\beta = 30$ (their Table E.1) the error falls from 0.738 to 0.105 with $\Delta t = 0.1$.

3. Causal weighting (Wang, Sankaran and Perdikaris). Split the time interval into $N_t$ ordered chunks $t_1 < \dots < t_{N_t}$, and let $L_r(t_i)$ be the mean squared residual of the collocation points in chunk $i$. Replace $L_r$ by

Chunk $i$ gets a large weight only after all earlier chunks have small residuals. The weights are computed from the current losses but treated as constants: no gradient flows through $w_i$ (stop_gradient in JAX, .detach() in PyTorch). The causality parameter $\varepsilon > 0$ sets how strict this is: too small and nothing is enforced, too large and the late chunks never switch on. The paper increases $\varepsilon$ through $10^{-2}, 10^{-1}, 1, 10, 100$ and moves to the next value when $\min_i w_i > 0.99$, a useful convergence check.

Interactive figure: loading…

Supporting tricks (Wang et al., An Expert's Guide to Training PINNs, 2023):

• Exact periodic boundary conditions: feed $(\cos x, \sin x)$ instead of $x$. Every function of these inputs is exactly $2\pi$-periodic, so $L_{bc}$ disappears (Dong and Ni, 2021).
• Fourier features: map the inputs $z$ to $[\cos(Bz), \sin(Bz)]$ with a random matrix $B$, against spectral bias (Tancik et al., 2020).
• Loss balancing: every 100–1,000 steps set $\hat\lambda_k = \sum_j \lVert\nabla_\theta L_j\rVert / \lVert\nabla_\theta L_k\rVert$ for each loss term $k$ and update $\lambda_k$ by a moving average, so all weighted terms have the same gradient norm.
• Scaling and optimiser: non-dimensionalise; use Adam, often followed by L-BFGS (optimization 10).

With these tricks plus causal weighting, the guide reaches a relative $L^2$ error of $6.88\times10^{-4}$ for advection at speed 80 (§7.2); its ablations show that each part matters.

Burgers: steep fronts #

The residual $u_t + u u_x - \nu u_{xx}$ is nonlinear and of second order. For small $\nu$ the solution forms a steep front. The steady front $u = -U\tanh\big(Ux/(2\nu)\big)$, with $U$ the half-height of the jump, is an exact solution. It goes through 90% of its jump over a width of $4\nu\operatorname{artanh}(0.9)/U \approx 5.9\nu/U$: 0.019 for the classic PINN test case $\nu = 0.01/\pi$, $U = 1$ on $[-1, 1]$. The network must represent this steep feature (spectral bias again), and the collocation points must resolve it: uniform random points put only about 1% of the points there. Residual-based adaptive sampling adds points where $\lvert r\rvert$ is large (Lu et al., DeepXDE, 2021). As $\nu \to 0$ the front becomes a shock where the residual is not defined.

PINN or solver? #

For a forward problem with a known PDE on a simple domain, a classical solver wins. The pseudo-spectral Burgers solver of lesson 08 matches the exact Cole–Hopf solution to $2.5\times10^{-9}$ ($\nu = 0.1$, 128 grid points) in a fraction of a second on a laptop. For advection, the Fourier solver is exact up to round-off. The tuned PINN of the Expert's Guide reaches $6.88\times10^{-4}$ on 1-D advection after 200,000 Adam steps on a GPU. And a PINN solves one initial condition per training run, while a solver or a neural operator handles any new one.

Interactive figure: loading…

A PINN is a sensible choice when the problem is not a clean forward problem:

• inverse problems: unknown coefficients ($\beta$, $\nu$, a diffusivity field) or sources, estimated together with $u$ from sparse, noisy observations; the PDE acts as a regulariser (lesson 04);
• data together with a partly known PDE, an irregular domain, or no solver at hand;
• the residual as an extra loss on a neural operator (physics-informed neural operators).

The data-plus-residual loss is weak-constraint 4D-Var (lesson 06) in another form: a data misfit plus a penalty on model error. The weight of the residual term plays the role of the inverse model-error covariance $Q^{-1}$, and the network replaces the model trajectory.

Answers to section 1. (1) It fails. In Krishnapriyan et al. the relative error is close to 100% for $\beta > 10$ (0.750 at $\beta = 20$, 0.897 at $\beta = 30$; at $\beta = 10$ it is still $1.08\times10^{-2}$). The network is not too small: curriculum training of the same network reaches $2\times10^{-2}$. (2) The zero field. With $\lambda_{ic} = 1$ its loss is $0 + \tfrac12$ (the mean of $\sin^2 x$). The frozen initial condition has residual $\beta\cos x$, so $L_r = \beta^2/2 = 450$. At $\beta = 1$ both losses are $0.5$. A large $\beta$ makes "keep the initial condition" very expensive and "flatten everything" cheap.

4. Worked example #

(a) Three candidates at $\beta = 30$ ($\lambda_{ic} = 1$; means over the domain). The relative error $\lVert\tilde u - u\rVert/\lVert u\rVert$ is taken over the whole space-time domain.

candidate	$L_r$	$L_{ic}$	relative error
exact, $\sin(x - 30t)$	0	0	0
zero	0	0.5	1
frozen, $\sin x$	$\beta^2/2 = 450$	0	1.44

For the frozen candidate, step by step:

1. The mean over $x$ of $(\sin x - \sin(x - \beta t))^2$ is $1 - \cos\beta t$ (expand the square; the cross term has mean $\tfrac12\cos\beta t$).
2. Its mean over $t \in [0, 1]$ is $1 - \sin\beta/\beta$.
3. Divide by the mean of $u^2$, which is $\tfrac12$, and take the square root: $\sqrt{2(1 - \sin 30/30)} = 1.44$.

At $\beta = 1$ the same candidate has $L_r = 0.5$ and relative error $\sqrt{2(1 - \sin 1)} = 0.56$.

(b) Theorem 1 on a wrong-speed candidate. Take $\tilde u = \sin(x - 29t)$ while $\beta = 30$.

1. Its residual is $r = -29\cos(x - 29t) + 30\cos(x - 29t) = \cos(x - 29t)$, so $L_r = \tfrac12$ and $\lVert r(\cdot, s)\rVert = \sqrt\pi$ at every $s$.
2. The initial error is zero, so the bound gives $\lVert e(t)\rVert \le \sqrt\pi\,t$: a relative error of at most $t$.
3. The true error: $\lVert e(t)\rVert^2 = 2\pi(1 - \cos t)$, so the relative error is $\sqrt{2(1 - \cos t)} = 2\sin(t/2)$. At $t = 1$ it is $0.959 \le 1$.

A mean squared residual of $\tfrac12$ allows almost 100% error; for 1% you need $5\times10^{-5}$ (section 3).

(c) Causal weights by hand. Four chunks have residuals $(0.001,\ 0.02,\ 0.5,\ 2.0)$; the plain mean is 0.630. The sums of the earlier chunks are $(0,\ 0.001,\ 0.021,\ 0.521)$, so $w_i = e^{-\varepsilon\cdot\text{sum}}$:

• $\varepsilon = 1$: $w = (1,\ 0.999,\ 0.979,\ 0.594)$: a mild effect; the last chunk still trains.
• $\varepsilon = 10$: $w = (1,\ 0.990,\ 0.811,\ 0.0055)$: chunk 4 is switched off until chunk 3 is solved.
• $\varepsilon = 100$: $w = (1,\ 0.905,\ 0.122,\ 2\times10^{-23})$: even chunk 3 is nearly off.

The causal losses are $0.425$, $0.109$ and $0.020$: a larger $\varepsilon$ hides more of the late residual.

5. In code #

PyTorch is the usual tool for PINNs. The first block computes residuals by automatic differentiation with respect to the inputs, for closed-form candidates and for a network:

import math
import torch

torch.set_default_dtype(torch.float64)

def advection_residual(u_fn, x, t, beta):
    """r = u_t + beta u_x at the points (x, t), by reverse-mode AD with respect to the inputs."""
    x, t = x.detach().requires_grad_(True), t.detach().requires_grad_(True)
    u = u_fn(x, t)
    u_x, u_t = torch.autograd.grad(u, (x, t), torch.ones_like(u), create_graph=True)
    return u_t + beta * u_x

# a uniform space-time grid on [0, 2π) x [0, 1]
X, T = torch.meshgrid(torch.arange(64) * 2 * math.pi / 64, torch.linspace(0, 1, 51), indexing="ij")
x, t, beta = X.ravel(), T.ravel(), 30.0
candidates = {
    "exact  sin(x - βt)": lambda x, t: torch.sin(x - beta * t),
    "frozen sin(x)     ": lambda x, t: torch.sin(x) + 0 * t,
    "zero              ": lambda x, t: 0 * x + 0 * t,
}
for name, fn in candidates.items():
    L_r = advection_residual(fn, x, t, beta).pow(2).mean().item()
    L_ic = (fn(x[t == 0], t[t == 0]) - torch.sin(x[t == 0])).pow(2).mean().item()
    print(name, f"L_r = {L_r:.3g}   L_ic = {L_ic:.3g}")
# exact  sin(x - βt) L_r = 0   L_ic = 0
# frozen sin(x)      L_r = 450   L_ic = 0
# zero               L_r = 0   L_ic = 0.5

def burgers_residual(u_fn, x, t, nu):
    """r = u_t + u u_x - nu u_xx. The second derivative needs a second backward pass."""
    x, t = x.detach().requires_grad_(True), t.detach().requires_grad_(True)
    u = u_fn(x, t)
    u_x, u_t = torch.autograd.grad(u, (x, t), torch.ones_like(u), create_graph=True)
    u_xx = torch.autograd.grad(u_x, x, torch.ones_like(u_x), create_graph=True)[0]
    return u_t + u * u_x - nu * u_xx

nu = 0.01 / math.pi
front = lambda x, t: -torch.tanh(x / (2 * nu)) + 0 * t        # steady viscous front, U = 1
xs = torch.linspace(-1, 1, 201)
print(f"{burgers_residual(front, xs, torch.zeros(201), nu).abs().max().item():.0e}")   # 4e-15

torch.manual_seed(0)
net = torch.nn.Sequential(torch.nn.Linear(2, 32), torch.nn.Tanh(), torch.nn.Linear(32, 1))
u_net = lambda x, t: net(torch.stack([x, t], dim=-1)).squeeze(-1)
x0, t0, h = torch.tensor([0.3]), torch.tensor([0.4]), 1e-4
r_ad = burgers_residual(u_net, x0, t0, nu)
with torch.no_grad():                                          # the same residual by central differences
    u = lambda dx, dt: u_net(x0 + dx, t0 + dt)
    r_fd = ((u(0, h) - u(0, -h)) / (2 * h) + u(0, 0) * (u(h, 0) - u(-h, 0)) / (2 * h)
            - nu * (u(h, 0) - 2 * u(0, 0) + u(-h, 0)) / h**2)
print(f"{abs(r_ad - r_fd).item():.1e}")                        # 2.2e-10: AD agrees with finite differences
r_ad.pow(2).sum().backward()
print(all(p.grad is not None for p in net.parameters()))       # True: gradients reach every weight

Second, Theorem 1 on the wrong-speed candidate of worked example (b), at three times (NumPy):

import numpy as np

x = np.linspace(0, 2 * np.pi, 1024, endpoint=False)
norm = lambda f: np.sqrt(np.sum(f**2) * (2 * np.pi / len(x)))   # L2 norm over [0, 2π)
u_norm = norm(np.sin(x))                                         # sqrt(π) = 1.7725
s = np.linspace(0, 1, 2001)
r_norm = np.array([norm((30 - 29) * np.cos(x - 29 * si)) for si in s])   # ‖r(·, s)‖ = sqrt(π)
for T in [0.25, 0.5, 1.0]:
    err = norm(np.sin(x - 29 * T) - np.sin(x - 30 * T))
    bound = np.trapezoid(r_norm[s <= T], s[s <= T])
    print(T, f"{err / u_norm:.3f} <= {bound / u_norm:.3f}")
# 0.25 0.249 <= 0.250
# 0.5 0.495 <= 0.500
# 1.0 0.959 <= 1.000

Third, the two NTK facts of view (b), for a one-hidden-layer tanh network whose Jacobians we can write by hand: spectral bias for a plain regression, and the $\beta^2$ growth of the residual block at initialisation:

rng = np.random.default_rng(0)
H = 256                                     # hidden units; f = sum_j a_j tanh(w_j z + b_j) / sqrt(H)

def jac_1d(t, w, b, a):
    """Rows d f(t_i) / d(a, w, b) for a one-hidden-layer tanh network of one input."""
    th = np.tanh(np.outer(t, w) + b); s2 = 1 - th**2
    return np.hstack([th, a * s2 * t[:, None], a * s2]) / np.sqrt(H)

# (1) Spectral bias: the NTK fits cos(βt) more slowly as β grows.
w, b, a = 4 * rng.normal(size=H), rng.normal(size=H), rng.normal(size=H)
t = np.linspace(0, 1, 200)
J = jac_1d(t, w, b, a)
Theta = J @ J.T                             # the NTK on 200 points
lam_max = np.linalg.eigvalsh(Theta)[-1]
for beta in [1, 5, 10, 30]:
    y = np.cos(beta * t)
    print(beta, f"{(y @ Theta @ y) / (y @ y) / lam_max:.1e}")   # initial fitting speed, relative
# 1 9.4e-01   5 2.3e-01   10 2.0e-02   30 3.1e-03

# (2) Imbalance: the residual block of the kernel grows like β² against the initial-condition block.
w1, w2, b2, a2 = (rng.normal(size=H) for _ in range(4))
def rows_u(x, t):                           # d u / d(a, w1, w2, b) for u(x, t)
    th = np.tanh(np.outer(x, w1) + np.outer(t, w2) + b2); s2 = 1 - th**2
    return np.hstack([th, a2 * s2 * x[:, None], a2 * s2 * t[:, None], a2 * s2]) / np.sqrt(H)
def rows_r(x, t, beta):                     # d r / d(a, w1, w2, b) for r = u_t + beta u_x
    th = np.tanh(np.outer(x, w1) + np.outer(t, w2) + b2); s2 = 1 - th**2
    m, ds2 = w2 + beta * w1, -2 * th * s2   # r = sum_j a_j s2_j m_j / sqrt(H)
    return np.hstack([s2 * m, a2 * (ds2 * x[:, None] * m + beta * s2),
                      a2 * (ds2 * t[:, None] * m + s2), a2 * ds2 * m]) / np.sqrt(H)
xr, tr, x0 = 2 * np.pi * rng.random(300), rng.random(300), 2 * np.pi * rng.random(300)
J_ic = rows_u(x0, np.zeros(300))
for beta in [1.0, 10.0, 30.0]:
    J_r = rows_r(xr, tr, beta)
    print(beta, f"{(J_r**2).sum() / (J_ic**2).sum():.3g}")   # tr K_rr / tr K_uu
# 1.0 0.638   10.0 16.3   30.0 149

The ratio $y^\top\Theta y / (y^\top y\,\lambda_{\max})$ is the initial rate at which kernel gradient descent reduces the error $y$, relative to the fastest direction: about 300 times smaller for $\cos 30t$ than for $\cos t$. Both kernel blocks of the second part use 300 points, so the trace ratio compares the two loss terms fairly.

6. Where it shows up #

• Seismology. EikoNet (Smith, Azizzadenesheli and Ross, IEEE Transactions on Geoscience and Remote Sensing, 2021) trains a network on the residual of the eikonal equation $\lvert\nabla T\rvert = 1/v$ for first-arrival travel times $T$ in a 3-D velocity model $v$, and uses it for earthquake location and tomography. This is an inverse setting with a cheap residual, where a mesh-free representation helps.
• Weather and climate. As of 2026, the main ML weather models (GraphCast, ECMWF's AIFS) are trained on reanalysis data, and NeuralGCM trains through a differentiable dynamical core: none of them minimises a PDE residual (lesson 16). Residual losses appear as extra terms (physics-informed neural operators, conservation penalties), and a PINN with observation terms is a weak-constraint 4D-Var.
• Machine learning. The same spectral-bias analysis explains why coordinate networks (NeRF) use Fourier features. DeepXDE and jaxpi are widely used PINN libraries.
• Benchmarks. A fair benchmark compares with a tuned classical solver at equal accuracy or equal cost, and reports the error against a reference solution, not the training loss.

7. Common mistakes #

• Reporting the loss instead of the error. Compare with a reference solution; Theorem 1 shows how small the residual must be, and for chaotic PDEs even that is not enough.
• No create_graph=True for the input derivatives. The residual term then has no gradient. If the initial-condition term still trains, nothing crashes, but the PDE has no effect.
• ReLU for a PDE with second derivatives. $u_{xx} = 0$ almost everywhere, so the viscous term disappears. Use tanh, GELU, sine or another smooth activation.
• Unbalanced loss terms, soft periodic conditions when the exact inputs $(\cos x, \sin x)$ are available, and inputs on very different scales ($x$ in km, $t$ in s).
• Too few or fixed collocation points for a steep front. Resample every step, or add points where the residual is large.
• Gradients through the causal weights, or unsorted time chunks. The weights must be constants.
• A weak baseline (a coarse or untuned solver), and expecting generalisation: a PINN solves one problem instance; a new initial condition needs a new training run.

8. Check yourself #

1. Write the PINN loss for viscous Burgers on $[0, 2\pi)$ with periodic boundary conditions and initial condition $g$. Which derivatives does each collocation point need, and how are they computed?
   Answer

   $L = \frac1{N_r}\sum_i r_i^2 + \lambda_{ic}\frac1{N_0}\sum_j (u_\theta(x_j, 0) - g(x_j))^2$ with $r = u_t + u u_x - \nu u_{xx}$. Periodicity either by a penalty on $u(0, t) - u(2\pi, t)$ (and on $u_x$), or exactly with the inputs $(\cos x, \sin x, t)$. Each point needs $u_t$ and $u_x$ (one reverse sweep, because the output is a scalar) and $u_{xx}$ (a reverse sweep through the computation of $u_x$). Training differentiates all of this with respect to $\theta$, so the inner sweeps are recorded (`create_graph=True`).
2. For advection with $\beta = 30$ on $[0, 2\pi) \times [0, 1]$, the residual norm $\lVert r(\cdot, s)\rVert$ is about the same at all times and the initial condition is exact. What mean squared residual is enough for a relative error below 1% at $t = 1$?
   Answer

   By Theorem 1, $\lVert e(1)\rVert \le \rho$ with $\rho = \lVert r\rVert$. A relative error below 1% needs $\rho \le 0.01\sqrt\pi = 0.0177$. The mean squared residual is $\rho^2/(2\pi) = 10^{-4}\pi/(2\pi) = 5\times10^{-5}$. The answer does not depend on $\beta$.
3. Three time chunks have residuals $(0.2, 0.2, 0.2)$ and $\varepsilon = 5$. Compute the causal weights and the causal residual loss.
   Answer

   Earlier sums $(0, 0.2, 0.4)$, so $w = (1, e^{-1}, e^{-2}) = (1, 0.368, 0.135)$. The loss is $\frac13 \cdot 0.2\,(1 + 0.368 + 0.135) = 0.100$, against the plain mean $0.2$.
4. Why is the moving-frame argument of view (a) not a general fix, and what does it show?
   Answer

   For linear advection with a known $\beta$ the moving frame solves the problem exactly, so there is nothing to learn. For Burgers the characteristics depend on the unknown solution, so you cannot change coordinates in advance. It shows that the PINN loss for advection is equally well posed for every $\beta$; what gets harder is representing and optimising the solution in the coordinates $(x, t)$.
5. A colleague's Burgers PINN uses ReLU and gives a smooth but wrong solution with a low loss. What happened?
   Answer

   A ReLU network is piecewise linear in $(x, t)$, so $u_{xx} = 0$ almost everywhere (and autograd returns 0 at the kinks). The residual the network sees is $u_t + u u_x$: it trains on the inviscid equation, and the viscous term has no influence on the result.
6. When would you prefer a PINN to the spectral solver of lesson 08 for Burgers?
   Answer

   When it is not a clean forward problem: for example $\nu$ (or a forcing) is unknown and you have only sparse, noisy observations of $u$. The PINN then fits $u$ and $\nu$ together, with the residual as a regulariser. For the forward periodic problem the solver is far more accurate ($10^{-9}$) and far cheaper.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Write the PINN loss for $u_t = \mathcal N[u]$ with initial condition $g$.

   Show answer

2. How is $u_{xx}$ of a network computed for a PINN residual?

   Show answer

3. State Theorem 1 (advection, periodic) of this lesson.

   Show answer

4. Advection at $\beta = 30$, $u(x, 0) = \sin x$: what is $L_r$ of the frozen field $\sin x$?

   Show answer

5. Why does a large $\beta$ cause loss imbalance?

   Show answer

6. What is spectral bias?

   Show answer

7. Give the causal weight of time chunk $i$.

   Show answer

8. Name three fixes for the PINN failure on fast advection.

   Show answer

9. When does a PINN make sense instead of a classical solver?

   Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Why do PINNs fail on advection at high speed? Is the network too small?" — No: curriculum training of the same network works. A residual term that grows like $\beta^2$, spectral bias and no causality in time make the optimisation hard.
• "How would you fix it?" — Curriculum in $\beta$, time marching, causal weights; exact periodic inputs, Fourier features, gradient-norm loss balancing, non-dimensionalisation, L-BFGS at the end.
• "Your PINN loss is $10^{-3}$. Is the solution accurate?" — Only through a stability estimate (Theorem 1 for advection, a Grönwall factor for nonlinear PDEs, nothing useful for chaos). Check against a reference solution.
• "PINN, neural operator or solver?" — Solver for forward problems; neural operator for many instances with training data; PINN for inverse problems with sparse data and a known PDE.
• "How do PINNs relate to data assimilation?" — Data misfit plus residual penalty is weak-constraint 4D-Var, with the network as the space-time trajectory.

10. Practice #

The notebook for this lesson runs in Colab or in the browser. It uses only NumPy (no deep-learning library) and has five exercises, each with a check cell and a folded solution:

1. PINN losses of candidate fields on a grid (exact, frozen, zero), from their derivatives.
2. Check the bound of Theorem 1 for a wrong-speed and a damped candidate.
3. Spectral bias: fit $\cos\beta t$ with a linear model on random tanh features by gradient descent and compare $\beta = 1$ with $\beta = 20$.
4. A random-feature PINN for advection, solved by linear least squares: failure at $\beta = 30$ and the time-marching fix.
5. The Fourier solver of lesson 08 for the same problem: error and cost against the PINN.

11. Watch and read #

Physics Informed Neural Networks (PINNs) [Physics Informed Machine Learning]Steve Brunton · 35 min

The PINN loss, the residual by automatic differentiation, and typical uses.

Physics-Informed Neural Networks: Failure Modes and SolutionsKoren Lab | Scientific Machine Learning · 11 min

The failure on convection and the curriculum and time-marching fixes of section 3.

Fourier Features Let Networks Learn High Frequency Functions in Low Dimensional DomainsMatthew Tancik · 3 min

Spectral bias and how Fourier features remove it, from the paper's author.

Physics Informed Machine Learning: High Level Overview of AI and ML in Science and EngineeringSteve Brunton · 47 min

Where physics enters a learning problem (data, architecture, loss, optimiser); PINNs are one option among several.

• Krishnapriyan et al., Characterizing possible failure modes in physics-informed neural networks (NeurIPS 2021): §3.1 (the convection failure), §4 (loss landscapes), §5.1 (curriculum, Table 1), §5.2 (sequence-to-sequence learning), App. B (the condition-number estimate).
• Wang, Sankaran and Perdikaris, Respecting causality for training physics-informed neural networks (CMAME 2024): §2 (the causality violation), §3 (causal weights and the choice of $\varepsilon$).
• Wang et al., An Expert's Guide to Training Physics-informed Neural Networks (2023): §3 (non-dimensionalisation), §4.2 (Fourier features), §5 (causality, loss balancing, curriculum), §7.2 (advection); code in jaxpi.
• Wang, Yu and Perdikaris, When and why PINNs fail to train: a neural tangent kernel perspective (J. Comput. Phys., 2022).
• Rahaman et al., On the Spectral Bias of Neural Networks (ICML 2019); Tancik et al., Fourier Features Let Networks Learn High Frequency Functions (NeurIPS 2020).
• Raissi, Perdikaris and Karniadakis, Physics Informed Deep Learning (Part I), the preprint of the 2019 J. Comput. Phys. PINN paper.
• McGreivy and Hakim, Weak baselines and reporting biases lead to overoptimism in machine learning for fluid-related PDEs (Nature Machine Intelligence, 2024).
• Thuerey et al., Physics-based Deep Learning (free online book), part "Physical Losses": physical loss terms and a Burgers PINN.