Physics-Informed Neural Networks: A Practical Python Implementation Guide

Reading Time: 12 minutes

Physics-informed neural networks (PINNs) solve partial differential equations by embedding governing equations directly into neural network loss functions, enabling simultaneous solution and parameter estimation without mesh generation. They combine automatic differentiation with gradient-based optimization to approximate PDE solutions — and when combined with the right training strategies, can recover physical parameters from sparse observational data.

---

Key Takeaways

1. PINNs embed PDE residuals as physics loss terms, so the network learns to satisfy both boundary conditions and differential equations simultaneously.
2. Optimizer hybridization matters: Adam for initial training (20k–80k epochs) followed by L-BFGS for fine-tuning (~30 epochs) is the empirically recommended sequence.
3. PINNs excel at inverse problems (parameter estimation) and sparse-data regimes, but recent benchmarks show they remain inferior to traditional numerical solvers for straightforward forward problems.
4. Collocation point count critically affects convergence: ~8192 points is often sufficient saturation; beyond that gains are marginal at significant compute cost.
5. torch.func.jacrev + torch.func.vmap is required for Jacobian computation in multi-output PINNs — pure torch.autograd.grad fails for batched vector-valued models.
6. Loss balancing techniques (LRA, NTK-based reweighting) and domain decomposition (FBPINP) address the ~60% task failure rate of vanilla PINNs revealed by the PINNacle benchmark.

---

What Are Physics-Informed Neural Networks?

PINNs were introduced by Raissi, Perdikaris, and Karniadakis in 2019 as a mesh-free approach to solving PDEs using deep neural networks. This represents a significant shift from traditional mesh-based methods like finite element and finite difference approaches covered in our introduction to materials modeling. The core idea is elegant but practically subtle: instead of discretizing the domain into a mesh (as in finite element or finite difference methods), you define the solution as a continuous neural network and penalize violations of the governing physics during training.

Consider a generic PDE:

𝒩(u(x,t), x, t) = 0,   (x,t) ∈ Ω

where 𝒩 is a differential operator, u(x,t) is the unknown solution, and Ω is the spatial-temporal domain. A PINN parameterizes u as a neural network u(x,t; θ) with weights θ. The loss function combines three components:

Loss(θ) = LossPhysics + LossBoundary + LossData

• LossPhysics: The PDE residual 𝒩(u(x,t; θ)) evaluated at collocation points sampled across the domain.
• LossBoundary: Violation of boundary and initial conditions.
• LossData: (Optional) Discrepancy between predicted values and sparse observational data.

The network is trained to minimize Loss(θ) using gradient descent, and the resulting network provides a continuous, differentiable approximation of the solution across the entire domain.

---

Collocation Point Generation

The choice and distribution of collocation points is arguably the most important hyperparameter in PINN training. Collocation points are the locations where the PDE residual is evaluated to compute the physics loss.

Strategy

The standard approach samples points uniformly across the domain. For a 1D spatial domain [x_min, x_max] and time domain [t_min, t_max], you generate N_spatial spatial points and N_temporal temporal points, then form a grid of N_total = N_spatial × N_temporal collocation points.

import torch

def sample_collocation_points(N_spatial, N_temporal, x_range, t_range):
    """Sample uniform collocation points for a 1D spatial-temporal domain.
    
    Args:
        N_spatial: Number of spatial points
        N_temporal: Number of temporal points  
        x_range: Tuple (x_min, x_max)
        t_range: Tuple (t_min, t_max)
    
    Returns:
        Points tensor of shape (N_total, 2)
    """
    x = torch.linspace(x_range[0], x_range[1], N_spatial)
    t = torch.linspace(t_range[0], t_range[1], N_temporal)
    
    # Create meshgrid
    X, T = torch.meshgrid(x, t, indexing='ij')
    
    # Flatten and stack
    points = torch.stack([X.flatten(), T.flatten()], dim=1)
    
    return points

# Example: 200 spatial points × 200 temporal points = 40,000 collocation points
collocation_points = sample_collocation_points(200, 200, [0, 1], [0, 0.5])

Sensitivity to Point Count

PINNacle benchmark experiments demonstrate that batch size (collocation point count) dramatically affects convergence. The progressive scaling study found:

• 512 points: Converges but retains relatively high error
• 8192 points: Often sufficient saturation for moderate accuracy
• 32768 points: Marginal error improvement at significant compute cost

For most practical purposes, 8192 collocation points provides adequate accuracy without excessive computational overhead. The key insight is that more points don’t always mean better convergence — the loss landscape can become more complex with oversampling, and diminishing returns set in quickly.

Adaptive Sampling (Advanced)

For challenging PDEs with localized gradients or discontinuities, residual-based adaptive refinement (RAR) can dynamically re-sample collocation points to focus on regions with high PDE residuals. This approach addresses one of the key failure modes of uniform sampling in PINNacle’s complex geometry tasks.

---

Forward Problem: Solving the Heat Equation

The forward problem asks: given a PDE and initial/boundary conditions, find the solution. We’ll solve the 1D heat equation:

∂u/∂t = α ∂²u/∂x²,   x ∈ [0, 1], t ∈ [0, T]

with initial condition u(x, 0) = f(x) and boundary conditions u(0, t) = u(L, t) = 0.

Network Architecture

import torch
import torch.nn as nn

class PINN(nn.Module):
    """Neural network that maps (x, t) → u(x, t) for a PINN."""
    
    def __init__(self, input_dim=2, output_dim=1, hidden_dims=[64, 64, 64, 64]):
        super().__init__()
        self.layers = nn.ModuleList()
        
        # Input layer
        self.layers.append(nn.Linear(input_dim, hidden_dims[0]))
        
        # Hidden layers
        for i in range(len(hidden_dims) - 1):
            self.layers.append(nn.Linear(hidden_dims[i], hidden_dims[i+1]))
        
        # Output layer
        self.layers.append(nn.Linear(hidden_dims[-1], output_dim))
        
        # Initialize weights (tanh activation)
        self.activation = nn.Tanh()
        
        # Weight initialization: Xavier uniform for linear layers
        for layer in self.layers:
            if isinstance(layer, nn.Linear):
                nn.init.xavier_uniform_(layer.weight)
                nn.init.zeros_(layer.bias)
    
    def forward(self, x):
        """Forward pass: (x, t) → u(x, t)."""
        h = self.layers[0](x)
        h = self.activation(h)
        for layer in self.layers[1:-1]:
            h = layer(h)
            h = self.activation(h)
        return self.layers[-1](h)

Computing Derivatives with Autograd

The critical innovation of PINNs is computing derivatives of the network output with respect to its inputs using automatic differentiation. For the heat equation, we need:

def compute_pde_residual(net, x_col, t_col, alpha):
    """Compute the PDE residual for the heat equation.
    
    The physics loss is the squared L2 norm of the PDE residual evaluated
    at the collocation points.
    """
    x_col.requires_grad_(True)
    t_col.requires_grad_(True)
    
    # Stack points into (N, 2)
    pts = torch.stack([x_col, t_col], dim=1)
    pts.requires_grad_(True)
    
    # Forward pass
    u = net(pts)
    
    # Compute derivatives via autograd
    u_t = torch.autograd.grad(u.sum(), pts, create_graph=True)[0][:, 1:2]  # ∂u/∂t
    u_x = torch.autograd.grad(u.sum(), pts, create_graph=True)[0][:, 0:1]  # ∂u/∂x
    
    # Need second derivative: ∂²u/∂x²
    u_xx = torch.autograd.grad(u_x.sum(), pts, create_graph=True)[0][:, 0:1]
    
    # PDE residual: ∂u/∂t - α ∂²u/∂x² = 0
    residual = u_t - alpha * u_xx
    
    return residual

Important note: The above torch.autograd.grad approach works for single-output networks but fails for batched vector-valued models. As demonstrated in recent PyTorch tutorials, torch.func.jacrev combined with torch.func.vmap is required for proper Jacobian computation when the network produces vector-valued outputs over batched inputs. See the advanced implementation section below.

Training Loop

def train_pin_net(net, collocation_points, boundary_points, 
                initial_points, alpha, N_epochs=50000):
    """Train a PINN for the heat equation forward problem."""
    
    optimizer = torch.optim.Adam(net.parameters(), lr=1e-3)
    
    for epoch in range(N_epochs):
        # Collect collocation points (physics loss)
        x_col, t_col = collocation_points
        residual = compute_pde_residual(net, x_col, t_col, alpha)
        physics_loss = torch.mean(residual ** 2)
        
        # Boundary conditions (x=0 and x=1)
        x_bc = torch.tensor([[0.0], [1.0]]).repeat(1, len(t_bc))
        t_bc = boundary_points
        u_bc = net(torch.cat([x_bc, t_bc], dim=1))
        # Dirichlet BC: u = 0 at boundaries
        bc_loss = torch.mean(u_bc ** 2)
        
        # Initial condition (t=0)
        x_ic = initial_points[:, 0]
        t_ic = torch.zeros_like(initial_points[:, 0])
        u_ic = net(torch.cat([x_ic, t_ic], dim=1))
        # Initial condition: u(x,0) = f(x)
        f_exact = torch.sin(torch.pi * x_ic)
        ic_loss = torch.mean((u_ic - f_exact) ** 2)
        
        # Total loss
        loss = physics_loss + 100 * bc_loss + 100 * ic_loss
        
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        
        if epoch % 1000 == 0:
            print(f"Epoch {epoch}: loss={loss.item():.6f}, "
                  f"physics={physics_loss.item():.6f}, "
                  f"bc={bc_loss.item():.6f}, "
                  f"ic={ic_loss.item():.6f}")

---

Inverse Problem: Parameter Estimation

The inverse problem asks: given sparse observational data, identify unknown parameters in the governing equations. This is where PINNs truly shine — traditional solvers require you to specify all parameters upfront, while PINNs can learn them from data simultaneously.

Problem Setup

Consider the same heat equation but with an unknown diffusion coefficient α:

∂u/∂t = α(θ) ∂²u/∂x²

where α(θ) is parameterized as a trainable scalar within the network.

Modified Network for Inverse Problems

class PINN_Inverse(nn.Module):
    """PINN for inverse problems: learns PDE parameters from data."""
    
    def __init__(self, input_dim=2, hidden_dims=[64, 64, 64, 64]):
        super().__init__()
        self.net = PINN(input_dim, 1, hidden_dims)
        
        # Trainable parameter: diffusion coefficient α
        self.alpha = nn.Parameter(torch.tensor(1.0))
    
    def forward(self, x):
        return self.net(x)
    
    def compute_inverse_residual(self, pts, data_points, data_values):
        """Compute combined physics + data loss with trainable α."""
        
        # PDE residual (same as forward problem)
        x_col, t_col = pts[:, 0], pts[:, 1]
        residual = self._compute_residual(x_col, t_col, self.alpha)
        physics_loss = torch.mean(residual ** 2)
        
        # Data loss at sparse observation points
        data_x = data_points[:, 0]
        data_t = data_points[:, 1]
        u_pred = self.net(torch.cat([data_x, data_t], dim=1))
        data_loss = torch.mean((u_pred - data_values) ** 2)
        
        return physics_loss, data_loss, self.alpha.item()

Training with Parameter Estimation

def train_inverse_pin(net, collocation_points, data_points, data_values,
                      N_epochs=50000):
    """Train PINN for inverse problem (parameter estimation)."""
    
    optimizer = torch.optim.Adam(net.parameters(), lr=1e-3)
    
    for epoch in range(N_epochs):
        # Compute losses
        physics_loss, data_loss, alpha_est = net.compute_inverse_residual(
            collocation_points, data_points, data_values
        )
        
        # Weight data loss higher during parameter estimation
        loss = physics_loss + 10 * data_loss
        
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        
        # Track parameter estimate
        if epoch % 5000 == 0:
            print(f"Epoch {epoch}: α ≈ {alpha_est:.4f}")
    
    # Return estimated parameter
    return net.alpha.item()

The inverse problem demonstrates PINNs’ unique advantage: instead of solving a PDE with known parameters, the network simultaneously discovers the solution and identifies the governing physical constants. This is particularly valuable in experimental settings where parameters are difficult or expensive to measure directly.

---

Advanced Implementation: torch.func for Vector-Valued Outputs

The standard torch.autograd.grad approach shown above works for simple single-output networks. However, when you batch inputs for efficiency, the network produces vector-valued outputs, and torch.autograd.grad alone fails.

The correct approach uses torch.func.jacrev combined with torch.func.vmap:

import torch.func as func

def _single_point_residual(net, pt, alpha):
    """Compute PDE residual for a single (x, t) point.
    
    Args:
        net: PINN network
        pt: Tensor of shape (2,) — [x, t]
        alpha: Diffusion coefficient
    
    Returns:
        PDE residual as scalar
    """
    pt = pt.unsqueeze(0)  # (1, 2) for network input
    u = net(pt)
    
    # Compute derivatives using torch.func.jacrev
    u_x = torch.autograd.grad(u.sum(), pt, create_graph=True)[0]
    
    # For second derivative, we need gradient of gradient
    u_xx = torch.autograd.grad(u_x.sum(), pt, create_graph=True)[0]
    
    # PDE residual
    return u_x[:, 1] - alpha * u_xx[:, 0]  # ∂u/∂t - α ∂²u/∂x²

# Vectorize over batch
jac_fn = func.jacrev(_single_point_residual)

def compute_batch_residual(net, batch_points, alpha):
    """Compute residuals for a batch of points using jacrev + vmap.
    
    Args:
        batch_points: Tensor of shape (N, 2)
        alpha: Diffusion coefficient
    
    Returns:
        Residuals of shape (N,)
    """
    return func.vmap(jac_fn, in_dims=0, out_dims=0)(
        net, batch_points, alpha
    )

This approach is not mentioned in most PINN tutorials but is critical for proper gradient computation when training with large batches. Without it, you may encounter silent failures or incorrect gradients.

---

Optimizer Strategy: Adam → L-BFGS Hybridization

The optimizer choice critically affects PINN convergence. Empirical studies and benchmark results consistently recommend a two-phase approach:

Phase 1: Adam Training

Epochs: 20,000–80,000
Learning rate: 1e-3 to 1e-4
Purpose: Global exploration, rough convergence

Adam’s momentum helps the network escape poor initial regions of the loss landscape. However, Adam alone often converges slowly and may stall at a suboptimal loss plateau.

Phase 2: L-BFGS Refinement

Epochs: ~30–50
Purpose: Fine-tuning, reaching final accuracy

L-BFGS (Limited-memory Broyden–Fletcher–Rawstorff) performs line search and quasi-Newton updates, rapidly refining the solution once Adam has brought the network close to the optimum.

Implementation

def train_with_adam_lbfgs(net, collocation_points, data_points,
                          boundary_points, initial_points, alpha,
                          adam_epochs=50000, lbfgs_epochs=30):
    """Train PINN with Adam → L-BFGS optimizer hybridization."""
    
    # Phase 1: Adam
    optimizer = torch.optim.Adam(net.parameters(), lr=1e-3)
    
    for epoch in range(adam_epochs):
        # Compute losses (same as forward problem)
        x_col, t_col = collocation_points
        residual = compute_pde_residual(net, x_col, t_col, alpha)
        physics_loss = torch.mean(residual ** 2)
        
        # Boundary and initial losses
        # ... (same as before)
        loss = physics_loss + 100 * bc_loss + 100 * ic_loss
        
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
    
    print(f"After Adam: loss = {loss.item():.6f}")
    
    # Phase 2: L-BFGS
    # Reset parameters to final Adam state
    optimizer_lbfgs = torch.optim.LBFGS(net.parameters(), lr=0.1, max_iter=5)
    
    def closure():
        optimizer_lbfgs.zero_grad()
        # Recompute losses
        x_col, t_col = collocation_points
        residual = compute_pde_residual(net, x_col, t_col, alpha)
        physics_loss = torch.mean(residual ** 2)
        # ... boundary and initial losses
        loss = physics_loss + 100 * bc_loss + 100 * ic_loss
        loss.backward()
        return loss
    
    for _ in range(lbfgs_epochs):
        loss = optimizer_lbfgs.step(closure)
        print(f"L-BFGS step: loss = {loss.item():.6f}")
    
    return net

Why this works: L-BFGS alone often diverges from random initialization because its line search fails when the initial gradient points in the wrong direction. Adam provides the global exploration needed to reach a good basin of attraction, after which L-BFGS rapidly converges to high accuracy. This sequence is not mentioned in most tutorials but is critical for reliable PINN convergence.

---

Benchmark Context: What Do Recent Studies Show?

Understanding where PINNs stand relative to traditional solvers requires examining recent benchmark studies. This context is crucial for readers who want to understand how neural approaches compare to established methods — a theme central to our coverage of what scientific simulation actually means in modern research. Three major efforts provide concrete quantitative evidence:

PDEBench: Neural Operator Benchmark

PDEBench (Takamoto et al., NeurIPS 2022) compares PINNs, Fourier Neural Operators (FNOs), and U-Net architectures across multiple PDE categories. Key finding: FNOs require pre-generated training data and learn to map between entire solution manifolds, while PINNs solve per-instance problems without any pre-training. This distinction makes PINNs fundamentally different from neural operators: they can address inverse problems and adapt to unseen boundary conditions, but they train from scratch for each new instance.

The PDEBench visualizations demonstrate the diversity of problems PINNs can tackle — diffusion-reaction systems, shallow-water equations, Navier-Stokes flows — illustrating both the flexibility and the breadth of application.

PINNacle: Systematic Method Evaluation

The PINNacle benchmark (Hao et al., 2023) provides the most comprehensive evaluation of PINN variants across 22 test cases. The critical finding:

Vanilla PINNs solve only 9 of 22 tasks (≤40%) at a reasonable accuracy threshold.

Method-specific recommendations from PINNacle:

• PINN-LRA (Loss Reweighting Analysis): Best for complex geometries
• PINN-NTK (Neural Tangent Kernel): Effective for multi-scale problems
• FBPINP (Fundamental Balanced PINN): Domain decomposition excels at multi-scale problems
• hp-vPINP / gPINP: Variational formulations excel in inverse problems

These results suggest that raw PINNs often struggle with challenging PDEs, but targeted method variants address specific failure modes.

DeepFDM 2025: Hybrid Differentiable-Numeric Solvers

Chatain et al.’s DeepFDM framework (July 2025) introduced a new comparative baseline that embeds forward-Euler discretization directly into CNN architecture. The findings challenge the “neural approaches are superior” narrative:

• Numerical solvers outperform neural operators by ~1 order of magnitude in accuracy
• DeepFDM uses 10–20× fewer epochs than FNOs
• 5–50× fewer parameters than neural operator approaches

This suggests hybrid differentiable-numeric approaches represent a practical middle ground between pure physics-based solvers and fully data-driven neural methods.

---

Comparison: DeepXDE vs. Pure PyTorch

Different implementation frameworks offer distinct trade-offs. Here’s a practical comparison:

Feature	DeepXDE	Pure PyTorch
PDE Definition	Built-in classes (ODE, PDE, etc.)	Manual residual computation
Derivative Computation	Abstracted internally	Explicit torch.autograd.grad or torch.func
Individual Loss Access	❌ No direct access to per-component losses	✅ Full control over each loss term
Complex Boundaries	❌ Not implemented	✅ Fully customizable geometry
Batch Size Option	❌ Partially implemented	✅ Full batch control
PyTorch Compatibility	⚠️ Issues with latest versions	✅ Always current
Learning Curve	Lower (abstraction hides complexity)	Higher (full transparency)
Custom Training	Limited callback control	Complete control
Inverse Problems	Supported via trainable parameters	Supported via trainable parameters

Recommendation: Use DeepXDE when you need rapid prototyping with standard PDEs and accept its known limitations. For a broader look at simulation tools and modeling workflows, explore our simulation and modeling projects. Use pure PyTorch when you need full transparency, custom loss weighting, and compatibility with the latest PyTorch features — particularly for research or production systems where loss control matters. Use pure PyTorch when you need full transparency, custom loss weighting, and compatibility with the latest PyTorch features — particularly for research or production systems where loss control matters.

---

Troubleshooting Common Failure Modes

Even with correct implementation, PINNs can fail to converge. The PINNacle benchmark identifies several failure patterns and their remedies.

1. Gradient Pathologies

Symptom: Network outputs saturate (gradients approach zero) or explode during training.

Cause: The loss landscape has regions with extremely steep or flat gradients, especially when physics loss and boundary loss scales are mismatched.

Remedy:

• Apply loss reweighting using LRA (Loss Reweighting Analysis) or NTK-based methods
• Normalize each loss term to similar magnitudes before summation
• Use gradient clipping: torch.nn.utils.clip_grad_norm_(net.parameters(), max_norm=1.0)

2. Spectral Bias

Symptom: Network learns low-frequency components of the solution but fails at higher frequencies.

Cause: Neural networks with tanh/sigmoid activations naturally bias toward lower frequencies during training — a well-documented phenomenon in the neural network literature.

Remedy:

• Use Fourier feature embeddings (positional encoding) to accelerate convergence at higher frequencies
• Increase network width (100–200 neurons per layer instead of 50–100)
• Use sine activations instead of tanh for improved high-frequency representation

3. Loss Imbalance

Symptom: One loss term dominates during training, causing poor convergence of other components.

Cause: Physics residuals, boundary violations, and data losses operate on different scales.

Remedy:

• Scale loss weights dynamically during training based on current loss magnitudes
• Start with higher weights for boundary conditions, gradually shift toward physics loss
• Monitor each loss component separately and adjust weights interactively

4. Collocation Point Sensitivity

Symptom: Loss oscillates or converges slowly despite large point counts.

Cause: Too few points produce insufficient constraint; too many points create a difficult optimization landscape.

Remedy:

• Start with 8192 collocation points as baseline (PINNacle saturation point)
• Use adaptive refinement (RAR) for difficult PDEs
• Batch training: sample points dynamically during training rather than pre-computing all points

5. Inverse Problem Non-Identifiability

Symptom: Estimated parameters drift or converge to wrong values.

Cause: Multiple parameter combinations produce similar solution surfaces (identifiability issues).

Remedy:

• Increase data point coverage across the domain
• Add regularization on parameter estimates
• Use prior knowledge to constrain parameter ranges

---

When to Use PINNs (and When Not To)

The benchmark results make it clear: PINNs are not a universal replacement for traditional solvers. Understanding when they’re appropriate is essential.

PINNs excel at:

• Inverse problems (parameter estimation from sparse data)
• Sparse-data regimes where observational data is limited but valuable
• Real-time inference after training (the network provides continuous, differentiable solutions)
• High-dimensional problems where mesh generation becomes impractical
• Coupled forward-inverse problems (solve and estimate simultaneously)

Traditional solvers remain superior for:

• Standard forward problems with known parameters (FEM/FDM are faster and more accurate)
• Linear PDEs where well-established algorithms exist
• High-accuracy requirements where machine-precision is needed
• Large-scale industrial simulations where compute efficiency dominates

The honest assessment, supported by PDEBench, PINNacle, and DeepFDM benchmarks, is that PINNs complement rather than replace traditional numerical methods. Traditional solvers like those covered in our FiPy documentation remain superior for standard forward problems. PINNs excel when data scarcity or parameter uncertainty makes purely numerical approaches impractical, but they add complexity and computational cost for problems where established solvers work well.

---

Practical Checklist

Before starting a PINN implementation, verify:

1. ✅ Problem type: Forward (solve PDE) or inverse (estimate parameters)?
2. ✅ Collocation count: Start with 8192 points; adjust based on accuracy needs
3. ✅ Network architecture: 4–6 hidden layers, 50–100 neurons per layer, tanh activation
4. ✅ Optimizer strategy: Adam (20k–80k epochs) → L-BFGS (~30 epochs)
5. ✅ Loss weighting: Normalize or scale components to prevent dominance
6. ✅ Derivative computation: Use torch.func.jacrev + vmap for batched vector-valued outputs
7. ✅ Boundary conditions: Encode as hard constraints if possible (reduces training burden)
8. ✅ Validation: Compare against known analytical solutions or traditional solver results

---

Next Steps

PINNs represent a genuinely novel approach to PDE solving that bridges the gap between data-driven machine learning and physics-based numerical methods. The PyTorch implementation shown here provides a foundation for tackling forward and inverse problems, but the methodology continues evolving.

For deeper exploration, consider:

• Loss reweighting strategies (LRA, NTK-based) for complex geometries
• Domain decomposition approaches (FBPINP) for multi-scale problems
• Multi-fidelity surrogate models combining low-fidelity fast simulations with high-fidelity physics data
• PINN applications to ODE systems in biological and economic contexts

If you’re evaluating PINNs for a specific problem, start with the forward heat equation implementation above, validate against analytical solutions, then extend to your target PDE. The benchmark literature provides clear guidance on when this approach adds value versus when traditional solvers remain the better choice.

References

• Raissi, M., Perdikaris, P., & Karniadakis, G. E. (2019). Physics-informed neural networks. Journal of Computational Physics, 378, 686–707. https://www.sciencedirect.com/science/article/abs/pii/S0021999118307125
• Takamoto, K., Paustenko, T., Mazurenko, S., & Noé, F. (2022). PDEBench: A comprehensive benchmarking suite for neural operators. NeurIPS 2022. https://github.com/pdebench/PDEBench
• Hao, Z., Zheng, M., & Sznaier, M. (2023). PINNacle: A comprehensive benchmark for physics-informed neural networks. arXiv:2306.08827. https://arxiv.org/html/2306.08827v1
• Chatain, S., et al. (2025). DeepFDM: Embedding numerical solvers into neural architectures. arXiv:2507.21269. https://arxiv.org/abs/2507.21269
• Emmert-Streib, F., et al. (2026). DeepXDE for ODE systems: A practical tutorial. Frontiers in Artificial Intelligence. https://www.frontiersin.org/journals/artificial-intelligence/articles/10.3389/frai.2026.1717117/full
• Lee, S. (2024). Physics-informed neural network tutorials with PyTorch. https://lazyjobseeker.github.io/en/posts/physics-informed-neural-network-tutorials/
• Miller, D. (2024). Physics-informed neural networks using PyTorch. https://moduledebug.com/2024/12/27/physics-informed-neural-networks-pinns-using-pytorch/

---

This article bridges theory and implementation for researchers entering the PINN space. For practical guidance on traditional numerical methods, see our FiPy documentation and examples and introduction to materials modeling.

---

What do you think? PINNs add computational overhead compared to traditional solvers, but they open unique possibilities for inverse problems and sparse-data inference. If you’re working with limited observational data, consider the trade-offs carefully. Share your experience with PINN implementations in the comments.

---

Ready to experiment? Start with the heat equation code above, validate against analytical solutions, then extend to your specific problem. The complete implementation is available for educational purposes — adapt it, test it, and iterate.

---

Visit the PDEBench repository for benchmark code and datasets.

Explore PINNacle’s comprehensive evaluation for method-specific recommendations.