🌊 Damped Harmonic Oscillator
True Solution
Analytical/Numerical Ground Truth
Standard Neural Network
Data-Only Training
Physics-Informed NN
Physics + Data Training
Ready to begin training...
Advanced multi-system demonstration with exact automatic differentiation for physical parameter discovery and scientific machine learning
Analytical/Numerical Ground Truth
Data-Only Training
Physics + Data Training
Let \(\mathcal{D}=\{(t_i,y_i)\}_{i=1}^N\) be (possibly scarce) observations of a state variable and let the governing law be expressed as a differential operator \(\mathcal{F}[u](t;\lambda)=0\), with unknown parameters \(\lambda\). A neural ansatz \(u_\theta(t)\) approximates the state, and the residual at a collocation point \(s_j\) is defined by
Training solves a regularized, multi-objective least-squares problem that balances data fidelity, physics residuals, and constraints (e.g., initial/boundary conditions):
Derivatives required by \(\mathcal{F}\) are computed exactly via automatic differentiation. In this demo, first- and second-order time derivatives \(\dot u_\theta,\, \ddot u_\theta\) are obtained by an order‑2 forward‑mode AD recurrence, while gradients of \(\mathcal{L}\) with respect to network weights are accumulated by reverse‑mode backpropagation. Parameters \(\lambda\) enter the residual analytically, yielding closed‑form physics gradients \(\partial \mathcal{L}/\partial \lambda\) without finite differences. Optimization uses Adam with light gradient clipping.
When \(\lambda\) is unknown, PINNs perform system identification by minimizing physics residuals. For the damped oscillator with unknown \((\omega,\zeta)\), the residual is
For the nonlinear pendulum with unknown \((g/L,\beta)\), the residual incorporates the sine nonlinearity:
The gradients \(\partial f_\theta/\partial \omega\), \(\partial f_\theta/\partial \zeta\), \(\partial f_\theta/\partial (g/L)\), and \(\partial f_\theta/\partial \beta\) have analytic forms, enabling stable updates of \(\lambda\) alongside \(\theta\). The nonlinear sine term requires automatic differentiation to compute \(\cos(u_\theta(t))\) contributions to parameter gradients. Similar constructions apply to the Lotka–Volterra rates \(\alpha,\beta,\gamma,\delta\).
For the heat equation with unknown thermal diffusivity \(\alpha\) and domain length \(L\), the residual is
The gradient \(\partial f_\theta/\partial \alpha = -\partial^2 u_\theta/\partial x^2\) is readily computed, while \(\partial f_\theta/\partial L\) involves the boundary condition dependencies. The spatiotemporal nature requires careful collocation point distribution across the \((x,t)\) domain to ensure adequate physics coverage.
The heat equation presents unique challenges as a parabolic PDE requiring 2D spatiotemporal neural approximation. The network \(u_\theta(x,t)\) operates in the 2D domain \((x,t) \in [0,L] \times [0,T]\) and must satisfy:
Physical Interpretation: This models heat diffusion in a rod with initial sinusoidal temperature distribution. The 2D nature comes from the spatiotemporal coupling - the network must learn how temperature varies across both space \(x\) and time \(t\) simultaneously. The thermal diffusivity \(\alpha\) controls the rate of heat spreading, while the boundary conditions represent fixed temperatures at the rod ends. The analytical solution \(u(x,t) = \sin(\pi x/L) \exp(-\alpha (\pi/L)^2 t)\) shows exponential temporal decay with spatial sinusoidal structure.
2D PINN Implementation: The physics loss samples collocation points uniformly across the 2D \((x,t)\) domain and computes the residual using automatic differentiation. The network architecture takes both spatial coordinate \(x\) and temporal coordinate \(t\) as inputs, learning the complete 2D solution field \(u(x,t)\). This is fundamentally different from 1D temporal ODEs, as the network must capture spatial patterns that evolve over time. Boundary conditions are enforced through penalty terms weighted by \(\beta\).
2D Learning Dynamics: The PINN must simultaneously learn the spatial diffusion pattern and temporal decay rate across the full 2D domain. The network learns to map any point \((x,t)\) to the corresponding temperature \(u(x,t)\), requiring it to understand both the spatial distribution at any given time and how that distribution evolves temporally. This 2D spatiotemporal learning is significantly more complex than 1D temporal problems.
Collocation points \(\{s_j\}\) probe the operator across the domain; sampling can be uniform, stratified, or adaptive. Initial and boundary conditions can be enforced softly via the penalty \(\mathcal{B}[u_\theta]\) or hard‑coded through constrained ansatzes. Proper scaling of \((\alpha,\beta)\), non‑dimensionalization, and normalization of states and time are often crucial for stable training.
The five exemplar systems span linear and nonlinear dynamics, single‑ and multi‑state settings, and temporal vs. spatiotemporal PDEs:
The pendulum model demonstrates true nonlinear dynamics where the small-angle approximation \(\sin(\theta) \approx \theta\) is not used, making it suitable for larger amplitude oscillations. The PINN must learn to handle the nonlinear \(\sin(\theta)\) term through automatic differentiation.
Heat Equation Physics: This is a 2D spatiotemporal problem where the heat equation models thermal diffusion in a rod of length \(L\) with thermal diffusivity \(\alpha\). The neural network \(u_\theta(x,t)\) learns the complete solution field across both spatial coordinate \(x\) and time \(t\). We solve on the domain \([0, L] \times [0, T]\) with homogeneous boundary conditions \(u(0,t) = u(L,t) = 0\) and initial condition \(u(x,0) = \sin(\pi x/L)\). The analytical solution exhibits exponential decay: \(u(x,t) = \sin(\pi x/L) \exp(-\alpha (\pi/L)^2 t)\).
Spatiotemporal PINN Challenges: Unlike 1D temporal ODEs, the heat equation requires the network \(u_\theta(x,t)\) to capture spatial and temporal dependencies simultaneously across a 2D domain. The physics residual \(R_\theta = \partial u_\theta/\partial t - \alpha \partial^2 u_\theta/\partial x^2\) involves mixed derivatives computed via automatic differentiation. Boundary conditions can be enforced through penalty terms or by architectural constraints. The time-dependent visualization reveals how the initial sinusoidal temperature profile diffuses and decays according to the thermal physics in the 2D \((x,t)\) space.
Wave Equation Theory: The 1D wave equation \(\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}\) describes propagating disturbances with wave speed \(c\). This hyperbolic PDE arises from Newton's second law applied to continuous media under tension, where \(u(x,t)\) represents displacement, \(x \in [0,L]\) is spatial position, and \(t \geq 0\) is time. The wave speed \(c = \sqrt{T/\rho}\) depends on tension \(T\) and linear density \(\rho\). For our system with \(L = 1\), \(c = 1.5\), and boundary conditions \(u(0,t) = u(L,t) = 0\), the analytical solution is \(u(x,t) = \sin(\pi x/L) \cos(c\pi t/L)\), exhibiting standing wave behavior with temporal period \(T = 2L/(c\pi) \approx 1.33\).
Mathematical Structure: The wave equation exhibits several fundamental properties: (1) Hyperbolicity: Characteristic curves \(x \pm ct = \text{const}\) determine information propagation at finite speed \(c\), unlike parabolic diffusion equations with infinite propagation speed. (2) Energy Conservation: Total energy \(E = \frac{1}{2}\int_0^L \left[(\partial u/\partial t)^2 + c^2(\partial u/\partial x)^2\right] dx\) remains constant in time, representing kinetic plus potential energy density. (3) D'Alembert Decomposition: General solutions decompose as \(u(x,t) = f(x-ct) + g(x+ct)\) representing rightward and leftward traveling waves. (4) Superposition Principle: Linear combinations of solutions are also solutions, enabling Fourier modal analysis.
PINN Computational Challenges: Wave equations pose unique difficulties for neural network approximation: (1) Second-order Time Derivatives: Computing \(\partial^2 u_\theta/\partial t^2\) requires two sequential automatic differentiation passes, creating deeper computational graphs with vanishing/exploding gradient risks. (2) Oscillatory Solutions: High-frequency temporal oscillations \(\cos(c\pi t/L)\) challenge neural networks' spectral bias toward low frequencies, requiring specialized architectures like Fourier features or curriculum learning. (3) Parameter Sensitivity: The quadratic coupling \(c^2\) in the PDE makes parameter gradients \(\partial r/\partial c = -2c \partial^2 u/\partial x^2\) scale nonlinearly with wave speed, demanding careful learning rate tuning. (4) Boundary Reflection: Standing wave patterns from boundary reflections create spatially varying solution complexity requiring adaptive collocation strategies.
Physics-Informed Training Methodology: Our implementation addresses these challenges through: (1) Multi-scale Collocation Sampling: Strategic point placement with 25% near high spatial gradient regions (\(x = L/4, 3L/4\) where \(|\partial u/\partial x|\) peaks), 25% near temporal inflection points (\(t = T/4, T/2, 3T/4\) where \(|\partial^2 u/\partial t^2|\) is large), and 15% in corners for boundary effects. (2) Energy Conservation Penalty: Additional loss term \(\lambda_E \sum_{i=1}^{N_t} |E(t_i) - E_0|^2\) penalizes energy violations, where \(E_0\) is initial energy and \(N_t\) time slices sample energy evolution. (3) Curriculum Learning: Progressive training from low-frequency (Phase 0: using 1/3 of collocation points) to mid-frequency (Phase 1: using 2/3) to full-frequency (Phase 2: all points) components with automatic transitions at epochs 300/600. (4) Enhanced Architecture: Residual connections \(u_\theta = \text{MLP}(x,t) + W_{\text{res},x} x + W_{\text{res},t} t + W_{\text{res},xt} xt\) improve gradient flow for oscillatory solutions.
Network Architecture: For single‑state ODEs, a network \(u_\theta(t)\) suffices. For multi‑state dynamics like Lotka–Volterra, one may use either a multi‑output network or separate surrogates for each state. For the 2D spatiotemporal PDEs (heat and wave equations), the network \(u_\theta(x,t)\) takes both spatial and temporal coordinates as input, learning the full solution field \(u(x,t)\) across the 2D \((x,t)\) domain. The wave equation requires enhanced architectures with residual connections to handle oscillatory solutions and second-order time derivatives.
PINN optimization is non‑convex. Useful practices include curriculum strategies for \(\alpha\), balanced batching of data and collocation points, early stopping on a held‑out physics set, gradient clipping, and diagnostics that monitor each loss component. Sensible units and non‑dimensionalization often improve conditioning and identifiability.
Physical Derivation: The 1D wave equation emerges from Newton's second law applied to a continuous string under tension \(T\) with linear mass density \(\rho\). Consider a small string element of length \(\Delta x\): the net vertical force from tension components is \(T[\sin\theta(x+\Delta x, t) - \sin\theta(x,t)] \approx T u_x(x+\Delta x, t) - T u_x(x,t) = T u_{xx} \Delta x\) for small angles. Newton's second law \(F = ma\) gives \(\rho \Delta x \cdot u_{tt} = T u_{xx} \Delta x\), yielding the wave equation \(\frac{\partial^2 u}{\partial t^2} = \frac{T}{\rho} \frac{\partial^2 u}{\partial x^2} = c^2 \frac{\partial^2 u}{\partial x^2}\) where the wave speed \(c = \sqrt{T/\rho}\) has units of velocity [m/s].
General Solution Structure: The d'Alembert solution for infinite domains is \(u(x,t) = f(x-ct) + g(x+ct)\) representing rightward and leftward traveling waves. For our bounded domain \([0,L]\) with fixed endpoints \(u(0,t) = u(L,t) = 0\), the method of separation of variables yields the modal expansion:
The coefficients \(A_n, B_n\) are determined by initial conditions \(u(x,0) = f(x)\) and \(u_t(x,0) = g(x)\). For our specific case with \(u(x,0) = \sin(\pi x/L)\) and \(u_t(x,0) = 0\), only the fundamental mode \(n=1\) is excited with \(A_1 = 1\), giving the exact solution \(u(x,t) = \sin(\pi x/L) \cos(c\pi t/L)\).
Energy Analysis: Wave energy has kinetic and potential components. The kinetic energy density is \(\rho u_t^2/2\), representing the motion of string elements. The potential energy density is \(T u_x^2/2\), representing the elastic energy stored in string stretching. The total energy per unit length is:
where we've used \(c^2 = T/\rho\). Total energy \(E(t) = \int_0^L e(x,t) dx\) is conserved: \(\frac{dE}{dt} = \int_0^L \frac{\partial e}{\partial t} dx = \int_0^L [u_t u_{tt} + c^2 u_x u_{xt}] dx = 0\) after integration by parts and applying the wave equation \(u_{tt} = c^2 u_{xx}\).
Characteristic Analysis: The wave equation is hyperbolic with characteristic curves \(x - ct = \xi\) and \(x + ct = \eta\). These represent the paths along which information propagates at speed \(c\). The canonical form in characteristic coordinates \((\xi, \eta)\) is \(\frac{\partial^2 u}{\partial \xi \partial \eta} = 0\), whose general solution is \(u = F(\xi) + G(\eta) = F(x-ct) + G(x+ct)\). This geometric interpretation explains why wave disturbances maintain their shape while traveling and why the domain of dependence for point \((x_0, t_0)\) is the interval \([x_0 - ct_0, x_0 + ct_0]\).
Boundary Condition Analysis: For fixed endpoints \(u(0,t) = u(L,t) = 0\), the boundary conditions impose constraints on the traveling wave components: \(F(-ct) + G(ct) = 0\) and \(F(L-ct) + G(L+ct) = 0\). These lead to the requirement that \(F\) and \(G\) are odd, periodic extensions with period \(2L\), creating the standing wave pattern \(u(x,t) = \sum_{n} [A_n \cos(\omega_n t) + B_n \sin(\omega_n t)] \sin(k_n x)\) with frequencies \(\omega_n = nc\pi/L\) and wavenumbers \(k_n = n\pi/L\).
Numerical and PINN Considerations: The wave equation's hyperbolic character creates unique computational challenges: (1) CFL Condition: Explicit finite difference schemes require \(\Delta t \leq \Delta x/c\) for stability, (2) Dispersion Errors: Discrete schemes introduce artificial frequency-dependent wave speeds, (3) Boundary Reflections: Improper boundary treatment can create spurious reflections. For PINNs, these translate to: (1) Temporal Resolution: Networks must resolve oscillations with period \(T = 2L/(c\pi) \approx 1.33\), (2) Gradient Pathologies: Second-order time derivatives create deeper AD graphs, (3) Standing Wave Patterns: Spatial nodes at boundaries and antinodes at \(x = L/2\) require careful collocation sampling to capture the full solution structure.
Theoretical Foundation: Wave equations represent hyperbolic PDEs with finite propagation speed, fundamentally different from parabolic (heat) or elliptic (Laplace) equations. The mathematical challenges arise from: (1) Characteristic Method: Information travels along characteristic curves \(dx/dt = \pm c\), creating directional dependencies that neural networks must capture. (2) Conservation Laws: Energy conservation \(\frac{\partial}{\partial t}\left[\frac{1}{2}(u_t^2 + c^2 u_x^2)\right] = 0\) and momentum conservation impose physical constraints on solutions. (3) Dispersion Relations: Fourier modes \(u_k(t) = A_k \cos(\omega_k t) + B_k \sin(\omega_k t)\) with \(\omega_k = ck\pi/L\) create multi-scale temporal dynamics spanning frequency ranges from DC to Nyquist. (4) Boundary Conditions: Dirichlet conditions \(u(0,t) = u(L,t) = 0\) enforce standing wave patterns with spatial nodes at boundaries and antinodes at \(x = L/2\).
PINN-Specific Mathematical Challenges: Beyond standard numerical methods, PINNs face unique difficulties: (1) Automatic Differentiation Depth: Computing \(\partial^2 u_\theta/\partial t^2\) requires AD through potentially deep networks, amplifying numerical errors and gradient pathologies. (2) Loss Landscape Complexity: The physics residual \(R = u_{tt} - c^2 u_{xx}\) creates non-convex optimization with multiple local minima corresponding to different frequency modes. (3) Parameter Coupling: The quadratic term \(c^2\) in the PDE makes parameter identifiability challenging - small changes in \(c\) create large changes in high-frequency content. (4) Spectral Bias: Neural networks preferentially learn low-frequency components, while wave solutions require accurate high-frequency representation for proper oscillatory behavior.
Research Background: Wave equations pose fundamental challenges for PINNs due to their second-order temporal nature and oscillatory solutions. This implementation incorporates 24 specific mathematical insights from recent literature to address these challenges, including energy conservation penalties based on Hamiltonian mechanics and frequency-aware training strategies inspired by spectral methods.
Core Challenges Identified:
Implemented Mathematical Solutions:
Mathematical Foundation - Parameter Gradient Analysis: For the wave equation \(\partial^2 u/\partial t^2 = c^2 \partial^2 u/\partial x^2\), the physics residual is \(r(x,t; c) = u_{tt} - c^2 u_{xx}\). The parameter gradient \(\frac{\partial L}{\partial c} = \frac{2}{N} \sum_{i=1}^N r_i \frac{\partial r_i}{\partial c}\) involves \(\frac{\partial r}{\partial c} = -2c u_{xx}\). This creates a quadratic coupling where the gradient magnitude scales as \(\|\nabla_c L\| \propto c \|u_{xx}\|_2\), making learning rates critically dependent on wave speed. The second-order spatial derivatives \(u_{xx}\) exhibit high variance near boundary layers and oscillatory regions, requiring specialized variance reduction through multi-scale sampling and gradient clipping with adaptive thresholds \(\tau = \min(1.0, \sigma_0/\sigma_t)\) where \(\sigma_t\) is current gradient variance.
Energy Conservation Mathematics: Wave energy \(E(t) = \frac{1}{2}\int_0^L [u_t^2 + c^2 u_x^2] dx\) should remain constant (\(\frac{dE}{dt} = 0\)) for conservative systems. Our energy penalty \(\mathcal{L}_E = \lambda_E \sum_{j=1}^{N_t} |E(t_j) - E_0|^2\) with \(\lambda_E = 0.1\) enforces this constraint by sampling energy at \(N_t = 5\) time slices across the domain \([0,T]\). The energy density integrand \(e(x,t) = u_t^2 + c^2 u_x^2\) combines kinetic and potential energy, where numerical integration uses spatial discretization with \(N_x = 10\) points and trapezoidal rule. This approach is mathematically justified by Noether's theorem: temporal translation symmetry of the Lagrangian \(\mathcal{L} = \frac{1}{2}(u_t^2 - c^2 u_x^2)\) leads to energy conservation as the corresponding conserved quantity.
Validation & Theoretical Basis: These insights combine classical wave theory with modern PINN methodology. The mathematical framework draws from: (1) Hamiltonian Mechanics: Energy conservation penalties from variational principles, (2) Spectral Theory: Frequency-aware loss weighting based on eigenfunction analysis of the wave operator \(\mathcal{L} = \partial^2/\partial t^2 - c^2 \partial^2/\partial x^2\), (3) Characteristic Theory: Collocation sampling aligned with characteristic directions \(x \pm ct = \text{const}\), and (4) Automatic Differentiation Theory: Gradient flow analysis for second-order derivative computations in deep networks. Implementation validation references include "Physics-Informed Neural Networks (PINNs) for Wave Propagation and Full Waveform Inversions" (American Geophysical Union) and "Energy-Conserving Neural Networks for Hamiltonian Systems" (Journal of Computational Physics).
ADNode Architecture: This implementation uses a custom automatic differentiation system based on dual numbers and forward-mode AD. Each variable is represented as an \(\texttt{ADNode}\) containing both the function value and its derivatives up to second order:
where \(v\) is the function value, \(\dot{v} = \frac{dv}{dt}\), and \(\ddot{v} = \frac{d^2v}{dt^2}\) for temporal problems, or \(\dot{v} = \frac{\partial v}{\partial x}\) and \(\ddot{v} = \frac{\partial^2 v}{\partial x^2}\) for spatial derivatives.
Forward-Mode Recurrence Rules: The AD system implements exact differentiation rules for neural network operations:
Network Forward Pass with AD: For a neural network \(u_\theta(t) = \sigma(W_L \sigma(W_{L-1} \cdots \sigma(W_1 t + b_1) + b_{L-1}) + b_L)\), the forward pass simultaneously computes:
Spatiotemporal AD for Heat Equation: For the 2D heat equation, the network \(u_\theta(x,t)\) requires mixed partial derivatives. The implementation uses separate AD passes:
Code Implementation Highlights: The \(\texttt{forward}\) method performs forward propagation where each layer operation maintains derivative tracking. The network architecture uses tanh activations, and each operation follows the AD chain rule. For example, the tanh activation implements:
Network Forward Pass Implementation: The \(\texttt{AutoDiffNN}\) class implements a single hidden layer network: \(u_\theta(t) = \mathbf{W}_2^T \tanh(\mathbf{W}_1 t + \mathbf{b}_1) + b_2\). The \(\texttt{forward}\) method returns an ADNode output that contains both the function value and its computational graph. Methods like \(\texttt{computeFirstDerivative}\) and \(\texttt{computeSecondDerivative}\) then call \(\texttt{ADNode.d1d2}\) to compute exact derivatives via the chain rule.
Physics Residual Computation: The physics loss directly uses AD outputs without finite differences. For the heat equation:
This provides machine-precision derivatives, avoiding numerical errors from finite difference approximations and enabling stable parameter identification even for stiff systems.