import time

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from tqdm import tqdm

plt.rcParams["text.usetex"] = True

import warnings
import torch
import torch.nn as nn

# Suppress Intel Iris Xe "not officially supported" warning
warnings.filterwarnings(
    "ignore",
    message=".*detected GPU.*not officially supported by PyTorch XPU.*")

# Define the GPU type ("cuda" for NVIDIA GPUs, "xpu" for Intel GPUs)
gpu_type = "xpu"

if gpu_type == "cuda":
    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
elif gpu_type == "xpu":
    device = torch.device("xpu" if torch.xpu.is_available() else "cpu")
else:
    raise ValueError("Invalid GPU type. Use 'cuda' or 'xpu'.")

print(f"Using device: {device}")

Using device: cpu

Parameterised PINN for Non-Linear ODEs

In the previous PINN example, we have evaluated a non-linear ODE based on a scenario from a damped pendulum of fixed length L and mass m, which swings under gravity g. Air resistance introduces a damping force proportional to angular velocity. The angle from vertical, \theta(t), obeys the following equation of motion:

\boxed{\frac{d^2\theta}{dt^2} + \frac{b}{m}\,\underbrace{\frac{d\theta}{dt}}_{=\,\omega} + \underbrace{\frac{g}{L}}_{=c_p}\sin\!\bigl(\theta(t)\bigr) = 0}

where b is the damping coefficient. We define \omega = d\theta/dt and c_p = g/L for convenience reason. In this example, we will look at parameterised PINN, where L is treated as an unknown parameter and the network learns \hat{\theta}(t, L) for any L \in [0.5, 1.5] m. The equation remains the same, but what changes is that L is now a flexible parameter.

Why do we need a single model over a range of L values?

Training a separate PINN for each L value would require running a full training loop (thousands of epochs) for every design point, which is expensive and inefficient. A single parameterised PINN is trained once and can then be queried instantly for any L in the range. This is particularly valuable in design optimisation, sensitivity analysis, or real-time control applications where many L values must be evaluated quickly.

Formulate Parameterised PINN Loss Functions for Non-Linear ODE

Term 1: Initial Condition (IC)/Boundary Condition (BC) Loss

\mathcal{L}_\text{IC} penalises the network for predicting the wrong starting angle or starting angular velocity at t=0. It is purely a boundary/initial constraint.

Constant L parameter in the ODE-model

With a fixed pendulum length, the network only takes t as input, so there is a single IC to satisfy:

\mathcal{L}_\text{IC} = \left(\hat{\theta}(0) - \theta_0\right)^2 + \left(\frac{d\hat{\theta}}{dt}\bigg|_{t=0} - \omega_0\right)^2

Flexible L parameter in the ODE-model

Because the network now takes two inputs (t, L), the starting conditions must hold for all L values in the training range. A single BC check at one L would leave the network free to predict wrong initial angles at other lengths. The sum over N_\text{bc} sampled values of L_j ensures the network learns to start at \theta_0 with velocity \omega_0 regardless of which pendulum length is queried:

\mathcal{L}_\text{BC} = \frac{1}{N_\text{bc}}\sum_{j=1}^{N_\text{bc}}\left[\left(\hat{\theta}(0, L_j) - \theta_0\right)^2 + \left(\frac{d\hat{\theta}}{dt}\bigg|_{t=0,\, L_j} - \omega_0\right)^2\right]

Initial condition vs. boundary condition

With a fixed L, the network input is one-dimensional (t only), so the constraint at t = 0 is a standard initial condition (IC): the solution must start at the prescribed angle and velocity at a single point in time.

When L becomes a flexible input, the network operates over a two-dimensional domain (t, L). The constraint at t = 0 must now be enforced for every L in the training range [L_\text{min}, L_\text{max}), forming a line (or boundary) in the (t, L) plane rather than a single point. This makes it a boundary condition (BC) in the broader sense: we are prescribing values along the t = 0 boundary of the 2D input space for all sampled L_j values.

Term 2: Physics (Residual) Loss

How to formulate the physics loss from a given ODE

Step 1: Start from the governing ODE. For the damped pendulum the equation of motion is

\boxed{\frac{d^2\theta}{dt^2} + \frac{b}{m}\,\frac{d\theta}{dt} + \frac{g}{L}\sin\!\bigl(\theta(t)\bigr) = 0}

Step 2: Replace the exact solution with the network prediction. Substitute \theta \to \hat{\theta} (the neural-network output) and compute every derivative via automatic differentiation. This gives the residual:

r(t) = \frac{d^2\hat{\theta}}{dt^2} + \frac{b}{m}\,\frac{d\hat{\theta}}{dt} + \frac{g}{L}\sin\!\bigl(\hat{\theta}(t)\bigr)

If \hat{\theta} were the true solution, r(t) would be exactly zero everywhere. In practice it is not, so r(t) measures how badly the network violates the physics at each point.

Step 3: Define the physics loss as the mean squared residual. Evaluate r(t) at a set of collocation points \{t_i\}_{i=1}^{N_\text{phys}} and penalise non-zero residuals:

\mathcal{L}_\text{phys} = \frac{1}{N_\text{phys}}\sum_{i=1}^{N_\text{phys}} r(t_i)^2

No numerical ODE solver is involved: the ODE is enforced purely through this loss term, and the physics is fully satisfied when \mathcal{L}_\text{phys} \to 0.

Constant L parameter in the ODE-model

With a fixed L, the coefficient c_p = g/L is a constant computed once before training. The residual only depends on t:

\mathcal{L}_\text{phys} = \frac{1}{N_\text{phys}}\sum_{i=1}^{N_\text{phys}} \left[\frac{d^2\hat{\theta}}{dt^2}(t_i) + \frac{b}{m}\,\frac{d\hat{\theta}}{dt}(t_i) + c_p\sin\!\bigl(\hat\theta(t_i)\bigr)\right]^2

where \mathcal{L}_\text{phys} is the physics loss that penalises the network for violating the equation of motion at the collocation points.

Flexible L parameter in the ODE-model

When L varies, c_p can no longer be precomputed as a single constant. Instead, each collocation point carries its own L_i, and the ratio g/L_i must be evaluated per point inside the residual. The network predictions \hat{\theta}, d\hat{\theta}/{dt}, and d^2\hat{\theta}/{dt^2} now all depend on both t_i and L_i:

\mathcal{L}_\text{phys} = \frac{1}{N_\text{phys}}\sum_{i=1}^{N_\text{phys}} \left[\frac{d^2\hat{\theta}}{dt^2}(t_i, L_i) + \frac{b}{m}\,\frac{d\hat{\theta}}{dt}(t_i, L_i) + \frac{g}{L_i}\sin\!\bigl(\hat\theta(t_i, L_i)\bigr)\right]^2

This ensures the ODE is enforced across the full (t, L) domain, not just for one pendulum length.

Total Loss

\mathcal{L} = \lambda_\text{BC}\,\mathcal{L}_\text{BC} + \lambda_\text{phys}\,\mathcal{L}_\text{phys}

Understanding parameterised PINN loss functions

(1) Why does c_p = g/L_i vary at every collocation point? Each collocation point (t_i, L_i) represents a different physical pendulum (a different length). The ODE must be satisfied correctly for each one, and c_p = g/L_i is what makes the ODE specific to that length. Using a fixed c_p for all points would only enforce the physics for one particular L, defeating the purpose of parameterisation.

(2) Why do we need many L values in the IC loss, not just one? The network must satisfy \hat{\theta}(0, L) = \theta_0 for every L in the training range, not just one. If we only enforced the IC at a single L value, the network could satisfy the initial conditions perfectly at that L while failing to start correctly for other lengths. Sampling many L values at t=0 ensures the IC constraint is broadcast across the entire parameter space.

(3) Why are more collocation points needed here (2,000) compared to the previous example (500)? The training domain is now 2D - (t, L) \in [0,1] \times [0.5, 1.5] - which is a much larger space to cover than the 1D time domain in the previous example. To ensure the physics is well-enforced throughout this 2D space, we need proportionally more collocation points. With too few, large regions of the parameter space would have no physics supervision and the network could learn wrong solutions there.

(4) Why can the network predict at L = 1.23 m even though it was never explicitly trained on that exact value? The PINN learns a smooth, continuous mapping \hat{\theta}(t, L). During training, collocation points sample the (t, L) space randomly, covering the range L \in [0.5, 1.5] densely but not exhaustively. Because the network weights define a smooth function and the physics loss enforces the ODE everywhere it is sampled, the learned solution interpolates naturally to any L within the training range, including L = 1.23 m. This is interpolation within the training distribution.

Building and Configuring the PINN

Model architecture

Two-element input: at each forward pass, the network receives a pair (t_i, L_i), the time value and the pendulum length at that collocation point, concatenated into a 2-element input vector. During evaluation, all time points share the same fixed L value (the chosen test length).
Tanh activation is still required: we still need to compute d^2\hat{\theta}/{dt^2} via automatic differentiation with respect to t, which requires the network to be twice differentiable in t. Tanh is infinitely smooth, so this is satisfied. We do not differentiate with respect to L (it is a parameter, not a variable to differentiate through), but smooth derivatives in t remain essential.
Minimal capacity change: adding one input neuron to the first layer increases the parameter count by only 64 weights (one extra weight per hidden neuron in the first layer) out of roughly 8,000+ total parameters. This is less than a 1% increase. The network’s representational capacity is determined by the width (64) and depth (3 layers), not by whether the input is 1D or 2D.

Define the Model Parameters

t_start = 0.0
t_end   = 1.0

g  = 9.81          # gravitational acceleration (m/s^2)
L  = 1.0           # pendulum length (m)
m  = 1.0           # mass (kg)
b  = 0.01          # damping coefficient (kg/s)

# Initial angle (radians)
theta0 = 1.0
  
# Initial angular velocity
# starts from rest, so angular velocity is zero
omega0 = 0.0

# L is now a range, not a single value
# minimum pendulum length (metres)
L_min  = 0.5     

# maximum pendulum length (metres)
L_max  = 1.5

Network architecture

We use a fully connected (dense) feedforward neural network with tanh activations. The architecture has three parts:

Layer block	Abbreviation	Description
`fcs`	fully connected start	Input layer: maps inputs to `N_NEURON` neurons, followed by tanh.
`fch`	fully connected hidden	Hidden layers: `N_LAYERS - 1` identical blocks, each a linear layer + tanh.
`fce`	fully connected end	Output layer: a single linear layer that maps from `N_NEURON` to the single output \hat{\theta}(t).

Here we use N_NEURON = 64 and N_LAYERS = 3, giving two hidden blocks and one input block (three tanh layers total).

Each nn.Linear(in, out) layer computes y = xW^\top + b, where W is an out \times in weight matrix and b is a bias vector of length out. For example, nn.Linear(1, 64) takes the single scalar input t and produces a 64-dimensional vector, one value per hidden neuron.

Number of layers

N_LAYERS represents the total number of layers with activations. The fcs block already provides the first activated layer, so fch only needs to create the remaining N_LAYERS - 1. With N_LAYERS = 3:

Block	Layers created	Tanh activations
`fcs`	1 (input → 64)	1
`fch`	`N_LAYERS - 1` = 2 (64 → 64 each)	2
`fce`	1 (64 → output)	0 (no activation)
Total	-	`N_LAYERS` = 3

class FullyConnectedNN(nn.Module):
    def __init__(self, N_INPUT, N_OUTPUT, N_NEURON, N_LAYERS):
        super().__init__()

        # Using Tanh activation
        activation = nn.Tanh

        # First layer: input to hidden layers
        # Using nn.Sequential to stack layers and activations
        self.fcs = nn.Sequential(
            nn.Linear(N_INPUT, N_NEURON), activation())

        # Hidden layers: N_LAYERS - 1 hidden layers
        # Since fcs already provides one, fch only needs to add the 
        # remaining N_LAYERS - 1.
        self.fch = nn.Sequential(*[
            nn.Sequential(nn.Linear(N_NEURON, N_NEURON), activation())
            for _ in range(N_LAYERS - 1)])

        # Output layer: hidden to output
        self.fce = nn.Linear(N_NEURON, N_OUTPUT)

    # Forward pass through the network
    def forward(self, x):
        x = self.fcs(x)
        x = self.fch(x)
        x = self.fce(x)
        return x

Training Setup

The network is instantiated with N_INPUT = 2 (instead of 1) so the first linear layer accepts the concatenated vector (t, L). Everything else in the architecture stays the same.

Collocation points (physics loss)

We sample 2,000 random (t, L) pairs from the 2D domain [0, t_\text{end}] \times [L_\text{min}, L_\text{max}). Both t and L are drawn uniformly. Only t_phys has requires_grad=True because we differentiate \hat{\theta} with respect to t (not L) when computing the ODE residual.

Boundary condition points (IC loss)

We sample 500 random L values from [L_\text{min}, L_\text{max}) and pair each with t = 0. This ensures the initial conditions \hat{\theta}(0, L_j) = \theta_0 and (d\hat{\theta}/{dt})|_{t=0,\, L_j} = \omega_0 are enforced across the full range of pendulum lengths, not just at a single L.

Uniform sampling from [a, b]

The notation [0, 1) is called a half-open interval. The square bracket [ means the left endpoint is included, and the round parenthesis ) means the right endpoint is excluded. So [0, 1) contains every number x satisfying 0 \le x < 1: zero is a possible value, but one is not. This is exactly the range produced by torch.rand.

Both L_phys and L_bc are generated using the pattern:

L_phys = torch.rand(N, 1) * (L_max - L_min) + L_min
L_bc = torch.rand(N, 1) * (L_max - L_min) + L_min

This samples uniformly from [L_\text{min}, L_\text{max}) in three steps:

torch.rand(N, 1) produces N random values in [0, 1).
Multiplying by (L_\text{max} - L_\text{min}) scales the range to [0,\, L_\text{max} - L_\text{min}).
Adding L_\text{min} shifts it to [L_\text{min},\, L_\text{max}).

In general: x = \text{rand}() \times (b - a) + a to sample uniformly from [a, b).

t_phys uses the same idea but with a = 0:

t_phys = torch.rand(N_phys, 1, device=device) * t_end

When a = 0 the formula simplifies to x = \text{rand}() \times b. torch.rand gives values in [0, 1), and multiplying by t_end scales them to [0, t_\text{end}). No shift is needed because the interval already starts at zero.

# Set random seed for reproducibility
torch.manual_seed(42)

# Note: n_input is now 2 to account for (t, L)
n_input   = 2

# n_output is still 1 since we're predicting theta(t, L)
n_output  = 1
n_hidden  = 64
n_layers  = 3
epochs    = 10000
learning_rate = 1e-3

# ----------------------------------------------------------------------------
# Create the PINN model
pinn = FullyConnectedNN(
    n_input,
    n_output,
    n_hidden,
    n_layers).to(device)

# Create the optimiser using Adam algorithm
optimiser = torch.optim.Adam(
    pinn.parameters(),
    lr=learning_rate)

# ----------------------------------------------------------------------------
# Collocation points: randomly sample (t, L) space
N_phys = 2000
# t_phys is sampled uniformly from [0, t_end)
t_phys = torch.rand(N_phys, 1, device=device) * t_end
t_phys.requires_grad_(True)

# L_phys is sampled uniformly from [L_min, L_max)
L_phys = torch.rand(
    N_phys, 1, device=device) * (L_max - L_min) + L_min

# ----------------------------------------------------------------------------
# Boundary points: t=0 for 500 random L values
N_boundary = 500

# L_bc is sampled uniformly from [L_min, L_max)
L_bc = torch.rand(
    N_boundary, 1, device=device) * (L_max - L_min) + L_min

# t_bc is set to 0 for all boundary points
t_bc = torch.zeros_like(
    L_bc, device=device, requires_grad=True)


print(f"Collocation points: {N_phys}   Boundary points: {N_boundary}")
print(f"Training for {epochs} epochs...")

Collocation points: 2000   Boundary points: 500
Training for 10000 epochs...

Inspecting PINN parameters in the network

pinn.parameters() is a generator that yields every learnable tensor (weight matrix or bias vector) in the network. For each parameter tensor p:

p.shape shows its dimensions (e.g. torch.Size([64, 2]) for the first weight matrix, which maps 2 inputs to 64 neurons).
p.numel() returns the total number of scalar values in that tensor (i.e. the product of all dimensions).

sum(p.numel() for p in pinn.parameters()) adds up the element counts across all tensors to give the total number of trainable parameters in the network.

With N_INPUT = 2, N_NEURON = 64, N_LAYERS = 3, and N_OUTPUT = 1, the loop prints eight tensors (one weight matrix and one bias vector per linear layer):

Block	Layer	Tensor	Shape	Parameters
`fcs`	`Linear(2, 64)`	weight	(64 \times 2)	128
-	-	bias	(64,)	64
`fch[0]`	`Linear(64, 64)`	weight	(64 \times 64)	4,096
-	-	bias	(64,)	64
`fch[1]`	`Linear(64, 64)`	weight	(64 \times 64)	4,096
-	-	bias	(64,)	64
`fce`	`Linear(64, 1)`	weight	(1 \times 64)	64
-	-	bias	(1,)	1
-	-	-	Total	8,577

Note that the fcs weight matrix has shape (64, 2) rather than (64, 1) because N_INPUT = 2. This is the only layer affected by the change from 1D to 2D input, adding just 64 extra weights compared to the single-input PINN (which has 64 \times 1 = 64 weights in the first layer, totalling 8,513 parameters).

for index, p in enumerate(pinn.parameters(), 1):
    print(f"Layer {index}: {p.numel()} parameters with shape: {p.shape}")

print("--" * 50)
print("Total parameters in the network:")
print(sum(p.numel() for p in pinn.parameters()))

Layer 1: 128 parameters with shape: torch.Size([64, 2])
Layer 2: 64 parameters with shape: torch.Size([64])
Layer 3: 4096 parameters with shape: torch.Size([64, 64])
Layer 4: 64 parameters with shape: torch.Size([64])
Layer 5: 4096 parameters with shape: torch.Size([64, 64])
Layer 6: 64 parameters with shape: torch.Size([64])
Layer 7: 64 parameters with shape: torch.Size([1, 64])
Layer 8: 1 parameters with shape: torch.Size([1])
----------------------------------------------------------------------------------------------------
Total parameters in the network:
8577

Training Loop

Each epoch performs three steps. The structure is the same as the fixed-L PINN, but every forward pass now receives a 2D input (t, L) instead of just t, and c_p = g/L is computed per collocation point rather than once before training.

Part 1: Boundary Condition Loss

Concatenate inputs: torch.cat([t_bc, L_bc], dim=1) joins the boundary time values (t = 0) and boundary lengths into an (N_\text{bc}, 2) tensor. Each row is one (0, L_j) pair. In the fixed-L version, the network received only t_bc (a scalar column of zeros).
Forward pass: pinn(inputs_bc) produces \hat{\theta}(0, L_j) for all N_\text{bc} boundary points simultaneously.
Differentiate w.r.t. t only: torch.autograd.grad(theta_bc, t_bc, ...) computes (d\hat{\theta}/{dt})|_{t=0,\, L_j}. The gradient is taken with respect to t_bc, not L_bc, because the IC constrains the time derivative at t = 0.
IC loss is the mean squared error over all boundary points:

\mathcal{L}_\text{bc} = \frac{1}{N_\text{bc}}\sum_{j=1}^{N_\text{bc}}\left[\left(\hat{\theta}(0, L_j) - \theta_0\right)^2 + \left(\frac{d\hat{\theta}}{dt}\bigg|_{t=0,\, L_j} - \omega_0\right)^2\right]

In the fixed-L version there was only one (t=0) point, so the loss was a single squared-error pair rather than a mean over many L values.

Part 2: Physics (ODE Residual) Loss

Concatenate inputs: torch.cat([t_phys, L_phys], dim=1) builds the (N_\text{phys}, 2) collocation input. Previously the network received only t_physics.
Forward pass: pinn(inputs_phys) produces \hat{\theta}(t_i, L_i) at all 2,000 collocation points.
First and second derivatives: two calls to torch.autograd.grad with respect to t_phys yield d\hat{\theta}/{dt} and d^2\hat{\theta}/{dt^2}. This is identical to the fixed-L version.
Per-point c_p: cp = g / L_phys produces a (N_\text{phys}, 1) tensor, giving each collocation point its own coefficient. In the fixed-L version, cp was a single Python float computed once before training.
ODE residual (element-wise, since cp is now a tensor):

r_i = \frac{d^2\hat{\theta}}{dt^2}(t_i, L_i) + \frac{b}{m}\,\frac{d\hat{\theta}}{dt}(t_i, L_i) + \frac{g}{L_i}\sin\!\bigl(\hat{\theta}(t_i, L_i)\bigr)

The physics loss is the mean squared residual: \mathcal{L}_\text{phys} = \frac{1}{N_\text{phys}}\sum_i r_i^2.

Part 3: Backpropagation

Combine the two losses: \mathcal{L} = \mathcal{L}_\text{BC} + \mathcal{L}_\text{phys} (both weights are 1 here).
Call loss.backward() to compute gradients of \mathcal{L} with respect to all network weights.
Call optimiser.step() to update the weights via Adam.

create_graph=True is passed to every torch.autograd.grad call so that the derivative computations remain part of the computational graph, allowing gradients to flow back through the physics loss during loss.backward().

total_loss_history = []
bc_loss_history = []
phys_loss_history = []

start_time = time.time()

for i in range(epochs + 1):
    optimiser.zero_grad()

    # ------------------------------------------------------------------------
    # Part 1: Boundary Condition (BC) Loss
    # concatenate t_bc and L_bc → shape (N_boundary, 2)
    # Each row is one (t=0, L_j) pair
    # the network sees both inputs simultaneously
    inputs_bc    = torch.cat([t_bc, L_bc], dim=1)
    theta_bc     = pinn(inputs_bc)
    dtheta_dt_bc = torch.autograd.grad(
        outputs=theta_bc,
        inputs=t_bc, 
        grad_outputs=torch.ones_like(theta_bc),
        create_graph=True)[0]

    loss_bc = torch.mean((theta_bc - theta0)**2) + torch.mean(
        (dtheta_dt_bc - omega0)**2)

    # ------------------------------------------------------------------------
    # Part 2: Physics (ODE Residual) Loss
    inputs_phys      = torch.cat([t_phys, L_phys], dim=1)
    theta_phys       = pinn(inputs_phys)
    dtheta_dt_phys   = torch.autograd.grad(
        outputs=theta_phys,
        inputs=t_phys,
        grad_outputs=torch.ones_like(theta_phys),
        create_graph=True)[0]
    d2theta_dt2_phys = torch.autograd.grad(
        outputs=dtheta_dt_phys,
        inputs=t_phys,
        grad_outputs=torch.ones_like(dtheta_dt_phys),
        create_graph=True)[0]

    # Using a different value at each collocation point
    # cp = g / L_phys 
    # This correctly encodes the physics for each individual pendulum length
    cp  = g / L_phys

    # ODE residual = d²θ/dt² + (b/m)·dθ/dt + cp·sin(θ)
    # cp is now a tensor (N, 1), so this is an element-wise operation
    res = d2theta_dt2_phys + (b/m)*dtheta_dt_phys + cp*torch.sin(theta_phys)

    loss_phys = torch.mean(res**2)

    # ------------------------------------------------------------------------
    # Part 3: Total Loss
    total_loss = loss_bc + loss_phys
    total_loss.backward()
    optimiser.step()

    # Store the loss values for plotting later
    total_loss_history.append(total_loss.item())
    bc_loss_history.append(loss_bc.item())
    phys_loss_history.append(loss_phys.item())

elapsed_time_pinn_training = time.time() - start_time
print(f"\nTraining complete in {elapsed_time_pinn_training:.4f} seconds.")
print(f"Final loss: {total_loss_history[-1]:.6f}")


Training complete in 91.9363 seconds.
Final loss: 0.000513

PINN Prediction

After training, we query the parameterised PINN at a single chosen pendulum length L_\text{test} (here 0.86 m, but any value in [L_\text{min}, L_\text{max}) works).

Choose a test length: L_test_val = 0.86 selects the pendulum length to evaluate.
Build the time grid: torch.linspace creates 100 evenly spaced t values over [t_\text{start}, t_\text{end}], reshaped to a (100, 1) column.
Create a constant L column: torch.full_like(t_pred, L_test_val) produces a (100, 1) tensor where every entry is 0.86. This pairs each time point with the same L.
Concatenate: torch.cat([t_pred, L_pred], dim=1) builds the (100, 2) input tensor the network expects, with each row being (t_i, 0.86).
Forward pass: inside torch.no_grad() (no gradient tracking needed at inference), the network returns \hat{\theta}(t_i, 0.86) for all 100 time points.

In the fixed-L PINN the network received only the time column t_pred. Here, because the network was trained on 2D inputs (t, L), we must always supply both dimensions, even when evaluating at a single L.

start_time_pred = time.time()

# Choose a test value of L within [L_min, L_max)
# Any value in [0.5, 1.5] is correct — 0.86 is the default
L_test_val = 0.86

# ----------------------------------------------------------------------------
# PINN prediction at the chosen L
t_pred  = torch.linspace(t_start, t_end, 100, device=device).view(-1, 1)
L_pred  = torch.full_like(t_pred, L_test_val)
inputs_pred = torch.cat([t_pred, L_pred], dim=1)

print("\nPredicted inputs (t, L):")
print(inputs_pred)
print(inputs_pred.shape)
print("--" * 50)

with torch.no_grad():
    theta_pred = pinn(inputs_pred).cpu().numpy()

elapsed_time_pinn_single_L_pred = time.time() - start_time_pred
print("Prediction complete in " 
      + f"{elapsed_time_pinn_single_L_pred:.4f} seconds.")


Predicted inputs (t, L):
tensor([[0.0000, 0.8600],
        [0.0101, 0.8600],
        [0.0202, 0.8600],
        [0.0303, 0.8600],
        [0.0404, 0.8600],
        [0.0505, 0.8600],
        [0.0606, 0.8600],
        [0.0707, 0.8600],
        [0.0808, 0.8600],
        [0.0909, 0.8600],
        [0.1010, 0.8600],
        [0.1111, 0.8600],
        [0.1212, 0.8600],
        [0.1313, 0.8600],
        [0.1414, 0.8600],
        [0.1515, 0.8600],
        [0.1616, 0.8600],
        [0.1717, 0.8600],
        [0.1818, 0.8600],
        [0.1919, 0.8600],
        [0.2020, 0.8600],
        [0.2121, 0.8600],
        [0.2222, 0.8600],
        [0.2323, 0.8600],
        [0.2424, 0.8600],
        [0.2525, 0.8600],
        [0.2626, 0.8600],
        [0.2727, 0.8600],
        [0.2828, 0.8600],
        [0.2929, 0.8600],
        [0.3030, 0.8600],
        [0.3131, 0.8600],
        [0.3232, 0.8600],
        [0.3333, 0.8600],
        [0.3434, 0.8600],
        [0.3535, 0.8600],
        [0.3636, 0.8600],
        [0.3737, 0.8600],
        [0.3838, 0.8600],
        [0.3939, 0.8600],
        [0.4040, 0.8600],
        [0.4141, 0.8600],
        [0.4242, 0.8600],
        [0.4343, 0.8600],
        [0.4444, 0.8600],
        [0.4545, 0.8600],
        [0.4646, 0.8600],
        [0.4747, 0.8600],
        [0.4848, 0.8600],
        [0.4949, 0.8600],
        [0.5051, 0.8600],
        [0.5152, 0.8600],
        [0.5253, 0.8600],
        [0.5354, 0.8600],
        [0.5455, 0.8600],
        [0.5556, 0.8600],
        [0.5657, 0.8600],
        [0.5758, 0.8600],
        [0.5859, 0.8600],
        [0.5960, 0.8600],
        [0.6061, 0.8600],
        [0.6162, 0.8600],
        [0.6263, 0.8600],
        [0.6364, 0.8600],
        [0.6465, 0.8600],
        [0.6566, 0.8600],
        [0.6667, 0.8600],
        [0.6768, 0.8600],
        [0.6869, 0.8600],
        [0.6970, 0.8600],
        [0.7071, 0.8600],
        [0.7172, 0.8600],
        [0.7273, 0.8600],
        [0.7374, 0.8600],
        [0.7475, 0.8600],
        [0.7576, 0.8600],
        [0.7677, 0.8600],
        [0.7778, 0.8600],
        [0.7879, 0.8600],
        [0.7980, 0.8600],
        [0.8081, 0.8600],
        [0.8182, 0.8600],
        [0.8283, 0.8600],
        [0.8384, 0.8600],
        [0.8485, 0.8600],
        [0.8586, 0.8600],
        [0.8687, 0.8600],
        [0.8788, 0.8600],
        [0.8889, 0.8600],
        [0.8990, 0.8600],
        [0.9091, 0.8600],
        [0.9192, 0.8600],
        [0.9293, 0.8600],
        [0.9394, 0.8600],
        [0.9495, 0.8600],
        [0.9596, 0.8600],
        [0.9697, 0.8600],
        [0.9798, 0.8600],
        [0.9899, 0.8600],
        [1.0000, 0.8600]])
torch.Size([100, 2])
----------------------------------------------------------------------------------------------------
Prediction complete in 0.0047 seconds.

Numerical Ground Truth at the same L

To validate the PINN’s prediction we compute a reference solution using scipy.integrate.odeint, a classical numerical ODE solver.

The second-order pendulum ODE is first rewritten as a system of two first-order ODEs by introducing the state vector [\theta,\;\omega] where \omega = \frac{d\theta}{dt}:

\frac{d\theta}{dt} = \omega, \qquad \frac{d\omega}{dt} = -c_p\sin\theta - \frac{b}{m}\,\omega

pendulum_system returns these two derivatives for a given state, and odeint integrates the system forward from the initial conditions [\theta_0,\;\omega_0] over 100 time points. The result theta_num is the numerically solved \theta(t), which we can plot alongside the PINN prediction to assess accuracy.

#  Numerical ground truth  #
def pendulum_system(states, t, b, m):
    theta, omega = states

    # d(theta)/dt = omega
    # d(omega)/dt = -cp*sin(theta) - (b/m)*omega   [rearranged ODE]
    return [omega, -(g / L_test_val)*np.sin(theta) - (b/m)*omega]

start_time_ode_single_L = time.time()

t_num   = np.linspace(t_start, t_end, 100)
sol_num = odeint(pendulum_system, [theta0, omega0],
                 t_num, args=(b, m))
theta_num = sol_num[:, 0]

elapsed_time_ode_single_L = time.time() - start_time_ode_single_L
print(f"ODE solution complete in {elapsed_time_ode_single_L:.4f} seconds.")

ODE solution complete in 0.0006 seconds.

Visualising the Results

fig, axes = plt.subplots(2, 1, figsize=(9, 10))

# ----------------------------------------------------------------------------
# Plot 1: PINN prediction vs numerical solution
ax1 = axes[0]
ax1.plot(
    t_num,
    theta_num,
    'r--',
    label=f"Numerical Solver (Ground Truth, L={L_test_val})",
    linewidth=2.5)
ax1.plot(
    t_pred.cpu().numpy(),
    theta_pred,
    'k-',
    label=f"PINN Prediction (L={L_test_val})",
    alpha=0.9)
ax1.set_title(
    f"Parameterised PINN with flexible length (L={L_test_val} m, b={b})",
    fontsize=16)
ax1.set_xlabel("Time (s)", fontsize=12)
ax1.set_ylabel(r"Angle $\theta$ (rad)", fontsize=12)
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# ----------------------------------------------------------------------------
# Plot 2: Training loss convergence
ax2 = axes[1]
ax2.plot(
    total_loss_history,
    marker='o',
    markersize=2,
    color='black',
    linewidth=3,
    label="Total Loss"
    )
ax2.plot(
    bc_loss_history,
    color='orange',
    linewidth=0.8,
    label="BC Loss")
ax2.plot(
    phys_loss_history,
    color='red',
    linewidth=1.2,
    alpha=0.8,
    label="Physics Loss")
ax2.set_yscale("log")
ax2.set_title("Parameterised PINN Training Loss", fontsize=16)
ax2.set_xlabel("Epoch", fontsize=12)
ax2.set_ylabel("Total Loss (log scale)", fontsize=12)
ax2.grid(True, which="both", linestyle="-", alpha=0.2)
ax2.legend(fontsize=13)

plt.tight_layout()
plt.show()

Why Are There Spikes in the Loss?

The training loss plot above shows dramatic, regularly-spaced spikes that persist throughout training. There are several contributing factors:

1. Gradient competition between loss terms (primary cause)

The two losses (i.e., BC loss and physics loss) pull the network in opposing directions. When Adam makes good progress minimising the physics residual, it can push the weights into a region where the boundary conditions are violated (and vice versa). With equal weighting (\lambda = 1 for both) and no adaptive balancing, these competing gradients periodically destabilise each other, producing the sawtooth pattern.

2. Fixed collocation points (no stochastic smoothing)

The 2,000 collocation points and 500 boundary points are sampled once before training and reused every epoch. There is no resampling. This means the loss landscape is deterministic and the optimiser can get trapped oscillating in narrow valleys. Periodic resampling of collocation points would introduce stochastic averaging that dampens these oscillations.

3. Constant learning rate with Adam’s momentum

Adam at a fixed lr = 1e-3 accumulates momentum. When progress on one loss term builds up momentum, the resulting step can overshoot and spike the other loss. Without learning rate decay or gradient clipping, there is nothing to dampen these overshoots as training progresses.

Possible mitigations:

Resample collocation points every N epochs to smooth the loss landscape.
Learning rate scheduling (e.g. cosine decay or ReduceLROnPlateau) to reduce step sizes as training progresses.
Adaptive loss weighting (e.g. the method of Wang et al., 2021) so the two loss terms remain balanced.
Gradient clipping (torch.nn.utils.clip_grad_norm_) to prevent large update steps.

Evaluating across multiple pendulum lengths

The key advantage of the parameterised PINN is that we can evaluate it at any L \in [L_\text{min}, L_\text{max}) without retraining. This cell demonstrates that by looping over five test lengths: L = 0.5, 0.75, 1.0, 1.25, 1.5 m.

For each L value the code:

Builds the PINN input: creates a constant L column with torch.full_like and concatenates it with the shared time grid to form a (300, 2) input tensor.
Queries the trained network: a single pinn(test_input) forward pass (inside torch.no_grad()) returns \hat{\theta}(t, L) for all 300 time points at that length.
Computes the numerical reference: odeint solves the pendulum ODE with c_p = g/L_\text{val} for the same time grid, providing the ground truth \theta(t) to compare against.
Plots both curves in the same colour (dashed for numerical, solid for PINN), so agreement between the two is immediately visible for each L.

Shorter pendulums (L = 0.5 m) oscillate faster because c_p = g/L is larger, while longer ones (L = 1.5 m) oscillate more slowly. The plot shows that the PINN captures this frequency variation across the full range from a single training run.

# Plot PINN predictions for multiple L values on the same graph
# No retraining needed, just change L_pred in the evaluation inputs

# all within training range
L_param_range   = np.linspace(L_min, L_max, 10).tolist()
t_range_pred     = torch.linspace(
    t_start, t_end, 300, device=device).view(-1, 1)

colors = plt.cm.viridis(np.linspace(0.15, 0.85, len(L_param_range)))

start_time_multi_pred = time.time()

L_param_pinn_pred_dict = {}

for L_val, col in zip(L_param_range, colors):

    # PINN prediction
    L_col   = torch.full_like(t_range_pred, L_val)

    # Concatenate t_range_pred and L_col to create input pairs (t, L_val)
    test_input     = torch.cat([t_range_pred, L_col], dim=1)

    # Evaluate the PINN at the new L value
    with torch.no_grad():
        pred = pinn(test_input).cpu().numpy()
        L_param_pinn_pred_dict[L_val] = pred

elapsed_time_multi_pred = time.time() - start_time_multi_pred
print("Multi-L PINN prediction complete in " 
      + f"{elapsed_time_multi_pred:.4f} seconds.")

Multi-L PINN prediction complete in 0.0042 seconds.

t_num_ext  = np.linspace(t_start, t_end, 300)

start_time_multi_num = time.time()

L_param_ode_num_dict = {}

for L_val, col in zip(L_param_range, colors):

    # Numerical ground truth at the same L value
    ode_sol = odeint(
        lambda s, t: [s[1], -(g / L_val)*np.sin(s[0]) - (b/m)*s[1]],
        [theta0, omega0], t_num_ext
    )

    L_param_ode_num_dict[L_val] = ode_sol[:,0]

elapsed_time_multi_ode = time.time() - start_time_multi_pred
print("Multi-L ODE computation complete in " 
      + f"{elapsed_time_multi_ode:.4f} seconds.")

Multi-L ODE computation complete in 0.0095 seconds.

fig, ax = plt.subplots(figsize=(9, 6))

colors = plt.cm.viridis(np.linspace(0.15, 0.85, len(L_param_range)))

for L_val, col in zip(L_param_range, colors):

    ax.plot(
        t_range_pred.cpu().numpy(),
        L_param_pinn_pred_dict[L_val],
        '-',
        color=col,
        linewidth=1.2,
        alpha=0.7,
        label=f'PINN L={L_val:.2f}')
    
    ax.plot(
        t_num_ext,
        L_param_ode_num_dict[L_val],
        '--', color=col, linewidth=2,
        label=f'Numerical L={L_val:.2f}')

ax.set_xlabel('Time (s)', fontsize=16)
ax.set_ylabel(r'$\theta$ (rad)', fontsize=16)
ax.set_title(
    'Parameterised PINN with multiple L Values', 
    fontsize=16)
ax.legend(
    fontsize=11, ncol=5,
    loc='upper center', bbox_to_anchor=(0.5, -0.12))
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Computational Time Comparison

The bar chart below compares the wall-clock time for each stage of the workflow. PINN training (10,000 epochs) is by far the most expensive step, but it is a one-time cost. Once trained, the parameterised PINN can predict \hat{\theta}(t, L) at any new L value with a single forward pass through the network without requiring any no ODE integration.

The key insight is in the prediction columns: querying the trained PINN at 10 different L values takes roughly the same time as querying it at one, because each evaluation is just a cheap matrix multiplication through the network. By contrast, the ODE solver must re-integrate the full system from scratch for every new L value, so its cost scales linearly with the number of test lengths.

This is precisely the advantage of the parameterised approach: the upfront training cost is amortised over all future evaluations. In settings where many L values must be explored, such as design optimisation, sensitivity analysis, or real-time control, the PINN offers orders-of-magnitude speedup over repeatedly calling a numerical solver.

labels = [
    "PINN Training\n(10k epochs)",
    "PINN Predict\n(single L)",
    "ODE Solve\n(single L)",
    "PINN Predict\n(10 L values)",
    "ODE Solve\n(10 L values)",
]

times = [
    elapsed_time_pinn_training,
    elapsed_time_pinn_single_L_pred,
    elapsed_time_ode_single_L,
    elapsed_time_multi_pred,
    elapsed_time_multi_ode,
]

bar_colors = ["black", "grey", "salmon", "grey", "salmon"]

fig, ax = plt.subplots(figsize=(9, 5))
bars = ax.bar(
    labels,
    times,
    color=bar_colors,
    edgecolor="black",
    linewidth=0.5)

for bar, t in zip(bars, times):
    ax.text(
        bar.get_x() + bar.get_width() / 2,
        bar.get_height() * 1.15,
        f"{t:.4f} s",
        ha="center",
        va="bottom",
        fontsize=12,
    )

ax.set_yscale("log")
ax.set_ylabel("Elapsed Time (s, log scale)", fontsize=14)
ax.set_title("Computation Time Comparison", fontsize=16)
ax.set_ylim(bottom=1e-4, top=max(times)*10)
ax.tick_params(axis="x", labelsize=14)
ax.grid(axis="y", alpha=0.3, which="both")
plt.tight_layout()
plt.show()

References

Francis Fernandes (2026). Mastering Dynamic PINNs.
Dao, Duc Long. “Experimental evaluation of damping models for a nonlinear pendulum system.” Physics Education 58, no. 5 (2023): 055003.