Physics-Informed Neural Networks
Why Embedding Physical Laws into Machine Learning Changes Everything

Why Traditional ML Struggles in Engineering

Large language models are trained on hundreds of billions of tokens of human text. The internet is saturated with language data. The statistical patterns in language are redundant, overlapping, and forgiving, if the model misses a nuance, the surrounding context usually compensates.

Now ask: how much labeled data does a pharmaceutical manufacturing process generate? A typical batch reactor might log temperature, pressure, and concentration at 10-second intervals. Run for a full year, that yields roughly 3 million data points, spread across dozens of operating variables, with process upsets, sensor drift, and changing feedstock composition corrupting a significant fraction of them. That is not internet-scale data. That is sparse, expensive, physics-constrained data.

The structural asymmetry between language and engineering data

Engineering systems are governed by differential equations

Almost every physical process in chemical, mechanical, or pharmaceutical engineering is described by an ODE or PDE. These equations encode centuries of accumulated physical understanding, not approximate statistical patterns, but precise mechanistic relationships.

Energy balance in a CSTR:

Diffusion-reaction in a catalyst pellet:

Navier-Stokes for incompressible flow:

A standard neural network trained only on data will violate these relationships. It will predict reactors with more energy out than in, fail completely outside the training range of temperatures and concentrations, and produce physically implausible results. This is not a theoretical concern, it has been observed repeatedly in industrial ML deployments.

The Three-Term Loss Function

A standard neural network learns a mapping $f_\theta: x \rightarrow y$ by minimizing the prediction error on labeled outputs:

A Physics-Informed Neural Network learns the same mapping, but adds the requirement that whatever the network predicts must also satisfy the governing differential equation. The training objective becomes:

In many engineering applications, and in the first example below, there are no labeled observations at all. The data term vanishes entirely, and the differential equation is the training signal.

A concrete example: the simplest first-order ODE

Consider:

Analytical solution: $y(x) = e^{-x}$. In the PINN approach, the network $\hat{y}(x)$ is evaluated at a set of collocation points $\{x_i\}$ scattered across the domain. At each point, the ODE residual is computed:

The physics loss forces this residual to zero across the entire domain:

The boundary condition loss pins the initial value:

No labeled $(x, y)$ pairs are needed anywhere in this training procedure. The differential equation is the training signal.

Why Autograd Wins Over Finite Differences

Enforcing $\frac{d\hat{y}}{dx} + \hat{y} = 0$ requires computing the derivative of the network's output with respect to its input. Three methods exist for this, with dramatically different tradeoffs:

When you compute $\hat{y} = f_\theta(x)$ in PyTorch with requires_grad=True on the input $x$, PyTorch constructs a computational graph recording every operation. Calling torch.autograd.grad() traverses this graph in reverse, applying the chain rule exactly, to machine floating-point precision. No mesh, no approximation, no truncation error.

The key parameter is create_graph=True on the first derivative call. Without it, PyTorch frees the intermediate graph and second derivatives become unavailable. With it, the derivative itself is a differentiable expression, enabling nested calls for $\frac{d^2\hat{y}}{dx^2}$, Laplacians, and beyond.

# Demonstrate exact automatic differentiation
# f(x) = x³ + 2x² − 5x + 1  →  f'(x) = 3x² + 4x − 5  →  f''(x) = 6x + 4

x = torch.tensor([2.0], requires_grad=True)

f = x**3 + 2*x**2 - 5*x + 1

# First derivative — create_graph=True keeps the graph alive for second pass
dfdx = torch.autograd.grad(f, x, create_graph=True)[0]

# Second derivative
d2fdx2 = torch.autograd.grad(dfdx, x)[0]

# Results at x=2: f=7, f'=15, f''=16  (all exact)
print(f.item(), dfdx.item(), d2fdx2.item())
# Output: 7.0  15.0  16.0

When the simple polynomial above is replaced by the full neural network $f_\theta(x)$, autograd can still propagate the chain rule through all layers, activation functions, and weight matrices, automatically and exactly. This is the mathematical foundation that makes PINNs possible.

Architecture, Collocation, and Training

Network architecture

For the first-order ODE problem, a compact architecture is sufficient: one input node (the scalar $x$), three hidden layers of 32 neurons each with Tanh activations, and a single output (the predicted $\hat{y}$). This yields 2,209 trainable parameters.

Why Tanh? Tanh is infinitely differentiable everywhere on $\mathbb{R}$. ReLU and its variants have zero second derivative almost everywhere, which breaks the higher-order autograd passes required for second-order ODEs and PDEs. Xavier normal initialization is used to ensure the gradients neither vanish nor explode at the start of training.

class PINN(nn.Module):
    def __init__(self, layers):
        super().__init__()
        network = []
        for i in range(len(layers) - 1):
            network.append(nn.Linear(layers[i], layers[i+1]))
            if i < len(layers) - 2:
                network.append(nn.Tanh())
        self.net = nn.Sequential(*network)
        self._init_weights()

    def _init_weights(self):
        for m in self.net:
            if isinstance(m, nn.Linear):
                nn.init.xavier_normal_(m.weight)
                nn.init.zeros_(m.bias)

    def forward(self, x):
        return self.net(x)

# Architecture: 1 → 32 → 32 → 32 → 1  (2,209 parameters)
model = PINN([1, 32, 32, 32, 1]).to(device)

Collocation points

Collocation points are the locations in the domain where the ODE residual is evaluated and penalized. They require no labels, only the $x$ coordinates. The key requirement is that requires_grad_(True) is set so that PyTorch can differentiate the network output with respect to these inputs.

The training loop and actual convergence

def ode_residual(model, x):
    y_hat = model(x)
    dy_dx = torch.autograd.grad(
        y_hat, x,
        grad_outputs=torch.ones_like(y_hat),
        create_graph=True     # required for loss to be differentiable w.r.t. weights
    )[0]
    return dy_dx + y_hat     # r(x) = dy/dx + y  →  should be zero

def total_loss(model, x_col, x_bc, y_bc):
    r = ode_residual(model, x_col)
    loss_physics = torch.mean(r**2)
    y_bc_pred = model(x_bc)
    loss_bc = torch.mean((y_bc_pred - y_bc)**2)
    return loss_physics + loss_bc

optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
for epoch in range(5000):
    optimizer.zero_grad()
    loss = total_loss(model, x_col, x_bc, y_bc)
    loss.backward()
    optimizer.step()

Training output from the notebook, showing convergence over 5,000 epochs:

Notice the characteristic pattern: the BC loss drops rapidly to near-zero within the first few hundred epochs, the network quickly learns to satisfy the single initial condition $y(0)=1$. The physics loss then dominates and drives the slow, steady improvement as the network learns to satisfy the ODE residual across the entire domain.

PINNs vs Classical Numerical Methods

Mesh-free computation

Creating a quality mesh for a complex 3D engineering geometry, a pharmaceutical tablet, a turbine blade, a reactor with internal baffles, can take days of preprocessing effort and requires specialized domain expertise. PINNs learn a continuous function over the domain without discretization, reducing geometry-to-simulation time dramatically.

Inverse problems without an outer optimization loop

In classical PDE solvers, estimating an unknown parameter (a reaction rate constant, a diffusion coefficient, a thermal conductivity) requires solving a forward problem thousands of times inside an outer optimization loop, computationally expensive and prone to getting stuck in local minima. In a PINN, the unknown parameter becomes an additional trainable variable in the same loss function and the same training loop. The forward solve and the parameter estimation happen simultaneously.

Data assimilation from sparse measurements

Adding 10 noisy experimental observations to the loss function, as a data term alongside the physics term, allows the PINN to incorporate real measurements while the physics automatically regularizes the solution between measurement points. The result is a physically consistent reconstruction that a pure data-driven model could not achieve with only 10 points.

Continuous and differentiable everywhere

A classical finite difference solution is a lookup table, values at fixed grid points. A PINN is a differentiable function defined for any $x$ in the domain. It can be evaluated at arbitrary resolution without retraining. Derivatives, integrals, and sensitivity analyses over the domain are all available analytically from the network itself.

Where PINNs Currently Struggle

Most introductory articles stop after the success stories. That would be a disservice to anyone planning to deploy PINNs in production. Here are the well-documented failure modes, with concrete mitigations where they exist.

Training instability

The PINN loss landscape is non-convex and poorly conditioned. Gradients from the physics loss and the boundary condition loss can point in conflicting directions, causing oscillations or premature plateauing. This is particularly severe when the physics loss and BC loss have different natural scales, one term dominates and the other is effectively ignored. Mitigation: adaptive loss weighting, learning rate scheduling, residual-based adaptive refinement (RAR) [Lu et al., 2021].

Stiff equations

A stiff ODE has solution components that evolve on vastly different timescales. For example, $\frac{dC}{dt} = -10^6 C$ decays to near-zero within microseconds. Standard PINN training with uniform collocation points entirely misses the fast dynamics, and the physics loss at late times provides negligible gradient information about the early transient. Mitigation: adaptive collocation density, curriculum training (gradually expanding the time horizon), domain decomposition.

Slow convergence compared to classical solvers

PINNs typically require tens of thousands of gradient descent steps to converge, orders of magnitude more computational work than a classical solver for equivalent accuracy on a simple well-posed problem:

PINNs are not faster solvers for classical forward problems. Their value is in other capabilities: inverse estimation, sparse data fusion, mesh-free handling of complex geometries.

Scaling to high dimensions and long time horizons

3D Navier-Stokes at high Reynolds numbers with complex boundary conditions remains largely unsolved by vanilla PINNs. Long time horizons compound the stiffness and spectral bias problems. High-dimensional PDEs suffer from the curse of dimensionality even in the collocation regime.

Multi-physics coupling

When heat, mass, momentum, and reaction terms are coupled, the relative magnitudes of the residuals differ by orders of magnitude. A naively assembled loss function will be dominated by the largest-magnitude residual, with the others ignored. Non-dimensionalization and learned or manually tuned loss weights are essential.

Beyond vanilla PINNs

Eight Articles on PINNs

The Jupyter notebooks for this series are available on GitHub: github.com/alishahmohammadi22/physics-ai-pinns. The notebook for this article is 01_Beginner/01_intro_to_pinn.ipynb.

Physics Becomes Part of the Network

The important idea is not that neural networks replace physics. The important idea is that physics becomes part of the neural network. For decades, engineers described systems through differential equations. For the last decade, AI described systems through data. PINNs represent the convergence of these two traditions.

In this article, we built a network that solved $dy/dx = -y$ to better than $10^{-3}$ absolute error using zero labeled data, with a 2,209-parameter network, 5,000 training epochs, and the ODE itself as the sole training signal. The BC loss vanished to $10^{-11}$ while the physics loss drove systematic improvement across the entire domain.

This convergence is already transforming pharmaceutical manufacturing (real-time parameter estimation from PAT data), computational fluid dynamics (mesh-free solvers for complex geometries), reaction engineering (simultaneous kinetics and transport identification from sparse pilot data), and scientific discovery (identifying governing equations from experimental observations). The tools are open-source, documented, and runnable on a standard laptop.

Further Reading