The 1D acoustic wave equation

Part 0 built the neural-network toolkit. Parts 1–3 turned that toolkit into a PINN engine: residual-based loss, AD for derivatives, architecture survey, training pathologies. Part 4 puts seismology into the engine. Every section in this Part is a wave equation — the underlying physics of every seismic experiment, from a single hammer-blow on a road shoulder to a 4D Sleipner survey.

The 1D acoustic wave equation

The simplest wave equation in seismology is the 1D acoustic wave equation,

with $u(x, t)$ the pressure perturbation (or the equivalent scalar wavefield in elastic settings under simplifying assumptions) and $c$ the wave-propagation speed. We drive the equation with an initial wavefield $u(x, 0) = u_0(x)$, an initial velocity $u_t(x, 0) = v_0(x)$, and boundary conditions on $u(x_{\textrm{left}}, t)$ and $u(x_{\textrm{right}}, t)$.

For this section we pick the simplest non-trivial test case:

• Domain $x \in [-1, 1]$, $t \in [0, 1]$.
• IC: $u(x, 0) = \sin(\pi x), \quad u_t(x, 0) = 0$.
• BC: $u(\pm 1, t) = 0$ (Dirichlet, "fixed walls").
• Speed $c$: a slider parameter the widget exposes.

The exact solution is the standing-wave mode

which we will use as the ground-truth reference. The PINN never sees this analytic answer during training; it sees only the PDE residual and the IC/BC samples. The relative-L² error against the analytic solution is the only metric that does not lie.

The PINN loss

The training loss is a four-term sum:

where, with $u_\theta(x, t)$ the network output,

• $\mathcal{L}${\textrm{PDE}} = \frac{1}{N_c} \sum_i \big( u{\theta,tt}(x_i, t_i) - c^2 u_{\theta,xx}(x_i, t_i) \big)^2,
• $\mathcal{L}${\textrm{IC,pos}} = \frac{1}{N{\textrm{IC}}} \sum_j \big( u_\theta(x_j, 0) - \sin(\pi x_j) \big)^2,
• $\mathcal{L}${\textrm{IC,vel}} = \frac{1}{N{\textrm{IC}}} \sum_j \big( u_{\theta,t}(x_j, 0) \big)^2,
• $\mathcal{L}${\textrm{BC}} = \frac{1}{2 N{\textrm{BC}}} \sum_k \big( u_\theta(-1, t_k)^2 + u_\theta(+1, t_k)^2 \big).

The four loss terms are weighted by NTK-balanced $\lambda$ values (§3.3 machinery: residual-Jacobian Frobenius norm, EMA-smoothed, clamped to $[10^{-2}, 10^{2}]$). This is mandatory: with naive $\lambda = 1$ on all four terms, the IC-velocity loss (a second-order term in disguise, since it requires $\partial_t u_\theta$ at $t = 0$) takes over and the PDE residual stalls.

Try it

Pick a wave speed $c$ between 0.5 and 2.0 and click ▶ Train. The widget runs 2500 epochs, then reports:

• The four per-term loss curves on a log-scale.
• Three time-slice snapshots of $u_\theta(x, t)$ at $t = 0, 0.5/c, 1/c$ — full crest, zero crossing, full trough.
• A side-by-side spacetime heatmap of the predicted wavefield versus the analytic solution.
• Relative-L² error over the $(x, t)$ grid.

What you should observe

• For $c = 1.0$ (the default) the PINN reaches relative-L² in the few-percent range. The standing wave is a low-frequency target (one spatial mode, one temporal mode); a small 2-32-32-1 Tanh MLP suffices.
• For $c = 2.0$ the temporal frequency doubles ($\cos(2\pi t)$ vs $\cos(\pi t)$). Spectral bias enters: the same architecture either takes longer to converge or saturates at higher relative-L². This is the §0.9 / §2.2 problem reasserting itself in seismic clothing.
• The IC-position term decreases first (it is a pure regression target), then the BC term, then the PDE residual catches up. The IC-velocity term is the slowest because it depends on $\partial_t u_\theta$ at $t = 0$, which the AD machinery has to expose to the gradient.
• The NTK-derived $\lambda$ values stabilise within a few hundred epochs; without NTK weighting the PDE residual stalls at $\sim 1$.

Why second-order time matters

The wave equation is second-order in time. That changes two things compared to the first-order parabolic problems (Burgers, advection) we met in earlier parts:

• The IC must specify both position and velocity. A first-order PDE needs only $u(x, 0)$; a wave equation needs $u(x, 0)$ AND $u_t(x, 0)$. The PINN loss has an extra term, and the architecture must produce a clean $\partial_t u_\theta$ at the boundary $t = 0$.
• The PDE residual involves $\partial_{tt} u_\theta$. Forward-mode AD computes Hessian diagonals; the $(x, t)$ input means the diagonal has two entries: $u_{xx}$ and $u_{tt}$. The PDE residual reads off both. Activation choice matters: Tanh is differentiable to all orders; ReLU dies at the second derivative. We are using Tanh.

From this widget to seismic-FWI

Real seismic forward modelling (and inverse problems built on top) replaces the toy IC with a Ricker wavelet source, replaces the Dirichlet BC with absorbing or PML BCs (§4.6, §4.7), goes from 1D to 2D or 3D, and uses spatially-varying $c(x, z)$ as the medium of interest. The PINN loss structure remains the same four-term recipe; the per-term machinery is what gets heavier. Section 4.2 turns to 2D, where the same construction starts to show why FDTD remains competitive (Part 5) and where PINN pays off in the inverse setting (Part 6).

References

• Raissi, M., Perdikaris, P., Karniadakis, G.E. (2019). Physics-informed neural networks: a deep learning framework for solving forward and inverse problems. JCP 378.
• Rasht-Behesht, M., Huber, C., Shukla, K., Karniadakis, G.E. (2022). Physics-informed neural networks (PINNs) for wave propagation and full-waveform inversions. J. Geophys. Res. Solid Earth 127.
• Moseley, B., Markham, A., Nissen-Meyer, T. (2020). Solving the wave equation with physics-informed deep learning. arXiv:2006.11894.