import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.nn as nn
from scipy.integrate import odeint
import time

device = torch.device("cpu")

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

PINN for Non-Linear ODEs

Let’s evaluate an example of non-linear ODE based on a scenario from a damped pendulum of 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}\,\frac{d\theta}{dt} + \frac{g}{L}\sin\!\bigl(\theta(t)\bigr) = 0}

where b is the damping coefficient. We define c_p = g/L for convenience.

Derivation of the equation of motion

At any instant t, in a frame of reference attached to m, there are four forces acting on the mass (see Figure 1):

Tension T in the suspension wire, directed along the rod toward the pivot A.
Gravity W:

W = mg. \tag{1}

Inertia force F_I: the mass m travels along a circular arc of radius L, so its tangential acceleration is L\frac{d^2\theta}{dt^2}. By D’Alembert’s principle the inertia force acting in the tangent direction is:

F_I = mL\frac{d^2\theta}{dt^2}. \tag{2}

How is F_I = mL\frac{d^2\theta}{dt^2} derived?

The mass m is constrained to move along a circular arc of radius L. The position along the arc (measured from the lowest point) is the arc length:

s = L\,\theta

Differentiating once with respect to time gives the tangential velocity:

v_\text{tan} = \frac{ds}{dt} = L\,\frac{d\theta}{dt}

Differentiating again gives the tangential acceleration:

a_\text{tan} = \frac{d^2 s}{dt^2} = L\,\frac{d^2\theta}{dt^2}

Applying Newton’s second law (F = ma) in the tangential direction:

F_I = m \cdot a_\text{tan} = m \cdot L\,\frac{d^2\theta}{dt^2} = mL\frac{d^2\theta}{dt^2}

This is simply the conversion from angular acceleration (\frac{d^2\theta}{dt^2}, in rad/s^2) to tangential linear acceleration (a_\text{tan}, in m/s^2) via the arc-length relation, then multiplied by mass to get the force.

Drag force F_D from the surrounding air. In general this is a function of the velocity v = L\frac{d\theta}{dt}:

F_D = f_D(v) = f\!\left(L\frac{d\theta}{dt}\right). \tag{3}

In practice it is usually assumed that F_D is linearly proportional to v. The drag force can therefore be expressed in terms of the angular velocity as:

F_D = bL\frac{d\theta}{dt}, \tag{4}

where b is the damping coefficient of the air, which depends on the viscosity of the air and the size and shape of the mass m.

What is D’Alembert’s principle?

D’Alembert’s principle reformulates Newton’s second law so that a dynamics problem can be treated as a static equilibrium problem. The idea is simple: instead of writing \sum F = ma, we move the inertia term ma to the left-hand side and treat it as a fictitious “inertia force” F_I = -ma acting on the body. The equation becomes:

\sum F + F_I = 0 \qquad \text{(i.e. } \sum F - ma = 0\text{)}

With this rearrangement, the sum of all forces (real forces plus the inertia force) equals zero, exactly as in a statics problem. This lets us apply equilibrium methods, such as resolving forces along chosen directions and setting their sum to zero to systems that are actually accelerating.

For the pendulum, the real forces are gravity mg, tension T, and drag F_D. The inertia force is F_I = mL\frac{d^2\theta}{dt^2} (opposing the tangential acceleration). D’Alembert’s principle says we can sum all of these along the tangent direction and set the total to zero, which directly yields the equation of motion.

Applying D’Alembert’s principle

Projecting all forces onto the tangent direction u (see Figure 1) and setting the sum to zero gives:

mL\frac{d^2\theta}{dt^2} + bL\frac{d\theta}{dt} + mg\sin\theta = 0. \tag{5}

Dividing through by mL:

\frac{d^2\theta}{dt^2} + \frac{b}{m}\,\frac{d\theta}{dt} + \frac{g}{L}\sin\theta = 0,

which is exactly the governing ODE above.

Why is this ODE non-linear?

An ODE is linear if the unknown function and all its derivatives appear only to the first power and are never composed with a nonlinear function. Inspecting each term:

Term	Explanation
\frac{d^2\theta}{dt^2}	linear in \theta
\frac{b}{m}\,\frac{d\theta}{dt}	linear in \theta
\frac{g}{L}\sin\!\bigl(\theta(t)\bigr)	nonlinear in \theta

The \sin(\cdot) function applied to the dependent variable \theta breaks linearity. If we replaced \sin(\theta) with \theta (the small-angle approximation), the ODE would become linear.

The pendulum starts at rest at an angle of 45^\circ:

\theta(0) = \frac{\pi}{4} \approx 0.785 \text{ rad}, \qquad \frac{d\theta}{dt}\bigg|_{t=0} = 0

Formulate PINN Loss Functions for Non-Linear ODE

Term 1: Initial Condition (IC) Loss

\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

\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.

Term 2: Physics (Residual) Loss

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

\mathcal{L}_\text{phys} penalises the network for violating the equation of motion, i.e. for predicting a \hat\theta(t) whose derivatives do not satisfy Newton’s second law at any of the collocation points.

A near-zero \mathcal{L}_\text{phys} means the network’s predicted trajectory satisfies the pendulum equation of motion at every collocation point. The network has learned a function that obeys Newton’s law for this system across the entire time window.

Total Loss

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

Building and Configuring the PINN

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)
cp = g / L         # c_p = g/L

# Initial angle (radians)
# 45 degrees converted to radians = pi/4
# 45 degrees = ~0.7854 rad
theta0 = np.pi / 4
  

# Initial angular velocity
# starts from rest, so angular velocity is zero
omega0 = 0.0

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 the single input t 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

Training points

t_ic: a single point at t = 0 used to enforce the initial conditions (\theta_0 and \omega_0).
t_phys: 500 evenly spaced collocation points over [0, 1] where the ODE residual is evaluated. Both tensors have requires_grad=True so that PyTorch can compute \frac{d\hat{\theta}}{dt} and \frac{d^2\hat{\theta}}{dt^2} via automatic differentiation.

Loss weights and optimiser

lambda_ic and lambda_physics are both set to 1.0, weighting the IC loss and physics loss equally.
The Adam optimiser is used with a learning rate of 10^{-3} for 5000 epochs.

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

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

pinn = FullyConnectedNN(
    N_INPUT=1, N_OUTPUT=1, N_NEURON=64, N_LAYERS=3
).to(device)

# Single point at t=0 (for initial conditions)
# Create a tensor of shape (1, 1) for the initial point
# Set requires_grad=True to enable gradient tracking for this tensor
t_ic = torch.tensor([[0.0]], device=device, requires_grad=True)
print(f"Initial point shape: {t_ic.shape}")

# 500 evenly spaced collocation points over [0, 1]
N_phys = 500

# Create a column vector of shape (N_phys, 1) for the collocation points
t_phys = torch.linspace(
    t_start, t_end, N_phys, device=device).view(-1, 1)

# Enable gradient tracking for the collocation points
t_phys.requires_grad_(True)
print(f"Collocation points shape: {t_phys.shape}")

# Set lambda_ic and lambda_physics to 1.0, 
# weighting the IC loss and physics loss equally.
lambda_ic = 1.0
lambda_physics  = 1.0

# Training parameters
epochs          = 5000
learning_rate   = 1e-3

# Using Adam optimiser for training the PINN
optimiser    = torch.optim.Adam(
    pinn.parameters(),
    lr=learning_rate)

print(
    f"Network created with {sum(p.numel() 
        for p in pinn.parameters())} parameters.")
print(f"Training for {epochs} epochs...")

Initial point shape: torch.Size([1, 1])
Collocation points shape: torch.Size([500, 1])
Network created with 8513 parameters.
Training for 5000 epochs...

Training Loop

Each epoch performs three steps:

Part 1: Initial Condition Loss

Pass t_ic (t = 0) through the network to get \hat{\theta}(0).
Use torch.autograd.grad to compute \frac{d\hat{\theta}}{dt}\big|_{t=0}, which is the predicted angular velocity at t = 0.
Compute the IC loss as the sum of squared errors against the known initial values:

\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

Part 2: Physics (ODE Residual) Loss

Pass all 500 collocation points t_phys through the network to get \hat{\theta}(t_i).
Differentiate twice using torch.autograd.grad to obtain \frac{d\hat{\theta}}{dt}\big|_{t_i} and \frac{d^2\hat{\theta}}{dt^2}\big|_{t_i}.
Evaluate the ODE residual at every collocation point:

r_i = \frac{d^2\hat{\theta}}{dt^2}\bigg|_{t_i} + \frac{b}{m}\,\frac{d\hat{\theta}}{dt}\bigg|_{t_i} + c_p\sin\!\bigl(\hat{\theta}(t_i)\bigr)

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

Part 3: Backpropagation

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

Note that create_graph=True is passed to torch.autograd.grad so that the derivative computations themselves are part of the computational graph, which allows gradients to flow back through the physics loss during loss.backward().

# List to store the history of loss values during training
total_loss_history = []
ic_loss_history = []
physics_loss_history = []

start_time = time.time()

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

    # -----------------------------------------------------------------------
    # Part 1: Initial Condition Loss
    theta_ic     = pinn(t_ic)
    dtheta_dt_ic = torch.autograd.grad(
        outputs=theta_ic,
        inputs=t_ic,
        grad_outputs=torch.ones_like(theta_ic),
        create_graph=True)[0]

    # Squared error between predicted angle at t=0 and theta0 
    loss_theta_ic  = torch.mean((theta_ic - theta0)**2)

    # Squared error between predicted velocity at t=0 and omega0
    loss_omega_ic  = torch.mean((dtheta_dt_ic - omega0)**2)

    # Total boundary loss is the sum of angle and velocity losses
    total_loss_ic  = loss_theta_ic + loss_omega_ic

    # -----------------------------------------------------------------------
    # Part 2: Physics (ODE Residual) Loss
    theta_phys       = pinn(t_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]

    # ODE residual = theta'' + (b/m)*theta' + cp*sin(theta)
    # This equals zero when the network perfectly satisfies the pendulum ODE
    r_physics  = (d2theta_dt2_phys 
        + (b/m)*dtheta_dt_phys 
        + cp*torch.sin(theta_phys))

    loss_physics = torch.mean(r_physics ** 2)

    # -----------------------------------------------------------------------
    # Part 3: Total Loss
    total_loss = lambda_ic * total_loss_ic + lambda_physics * loss_physics

    total_loss.backward()
    optimiser.step()
    total_loss_history.append(total_loss.item())
    ic_loss_history.append(total_loss_ic.item())
    physics_loss_history.append(loss_physics.item())

elapsed_time = time.time() - start_time
print(f"\nTraining complete in {elapsed_time:.2f}s." 
    + f" Final loss: {total_loss_history[-1]:.6f}")


Training complete in 28.70s. Final loss: 0.000014

Saving the Model Weights

After training, we save the learned parameters using torch.save(pinn.state_dict(), ...). The state_dict() method returns a dictionary mapping each layer name to its weight and bias tensors. Only these numerical values are stored and not the network architecture or training code, so the resulting .pth file is very small (typically a few tens of KB for a network of this size).

To reload the model later, you would recreate the same FullyConnectedNN architecture and call:

pinn.load_state_dict(torch.load("pinn_damped_pendulum_model_weights.pth"))

This restores the trained network exactly, ready for inference without re-training.

# Save the trained PINN model weights for later reuse
import os
model_path = "pinn_damped_pendulum_model_weights.pth"
torch.save(pinn.state_dict(), model_path)

# Get the file size of the saved model weights and print it in KB
file_size = os.path.getsize(model_path)
print(f"Model weights saved to '{model_path}' ({file_size / 1024:.1f} KB)")

Model weights saved to 'pinn_damped_pendulum_model_weights.pth' (37.3 KB)

PINN Prediction

After training, we evaluate the network on a dense grid of 200 evenly spaced points over [t_\text{start},\, t_\text{end}].

The torch.no_grad() context manager disables gradient tracking during inference as no derivatives are needed here, so this saves memory and speeds up the forward pass. The result theta_pred contains the PINN’s predicted \hat{\theta}(t) at each test point.

#  PINN prediction
start_time = time.time()
t_eval = torch.linspace(
    t_start, t_end, 200,
    device=device).view(-1, 1)

with torch.no_grad():
    theta_pred = pinn(t_eval)

elapsed_pred_time = time.time() - start_time
print(f"\nPrediction complete in {elapsed_pred_time:.4f} s.")


Prediction complete in 0.0012 s.

Numerical Ground Truth

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, cp):
    theta, omega = states

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

start_time = time.time()
t_num   = np.linspace(t_start, t_end, 100)
sol_num = odeint(pendulum_system, [theta0, omega0],
                 t_num, args=(b, m, cp))
theta_num = sol_num[:, 0]
elapsed_num_time = time.time() - start_time
print(f"\nNumerical solution complete in {elapsed_num_time:.4f} s.")


Numerical solution complete in 0.0004 s.

Visualising the Results

We produce two plots:

Top panel: overlays the PINN prediction (black solid line) on top of the numerical ground truth (red dashed line), showing \theta(t) over the full time window.
- Does the PINN prediction match the numerical solution well? The PINN prediction should closely follow the numerical solution across the full [0, 1] s window. There may be small discrepancies, typically the two curves overlap almost perfectly near t=0 (where the IC loss forces agreement) and remain close throughout. Any deviation is larger later in time, as small errors can compound. After 5,000 epochs the agreement is generally very good for this problem.
- Why is the numerical solver not truly exact? The odeint solver (based on LSODA, an adaptive Runge-Kutta / Adams method) uses a finite step size and makes local truncation errors at each step. It is not analytically exact, it produces a discrete approximation. However, for well-conditioned ODEs like this one with tight tolerance settings, the error is extremely small (typically below 10^{-8}), making it a reliable reference.
Bottom panel: plots three loss curves on a logarithmic y-axis against epoch number, revealing the convergence behaviour.
- Total Loss (black, thick): the weighted sum \mathcal{L} = \lambda_\text{IC}\,\mathcal{L}_\text{IC} + \lambda_\text{phys}\,\mathcal{L}_\text{phys}. This is the single scalar that the optimiser minimises at each step. It captures the overall training progress.
- IC Loss (orange): \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. This term measures how well the network satisfies the initial conditions. It typically drops quickly in early epochs because matching two scalar values (\theta_0 and \omega_0) is a relatively easy constraint.
- Physics Loss (red): \mathcal{L}_\text{phys} = \frac{1}{N_\text{phys}}\sum_i r_i^2, the mean squared ODE residual over all 500 collocation points. This term measures how well the predicted \hat{\theta}(t) satisfies the pendulum equation across the entire time domain. It usually starts higher and takes longer to converge than the IC loss, since it requires the network to learn the correct dynamics everywhere, not just at a single point.
- Why use a logarithmic scale? The loss drops by many orders of magnitude during training (e.g. from \sim 1 down to \sim 0.001). On a linear scale, the entire interesting early-training behaviour would be compressed into the left edge of the plot and invisible. A log scale spreads the detail across the full plot, making both rapid early convergence and slow late-stage refinement visible simultaneously.
- Approximately what is the loss at epoch 0? At epoch 5,000? The initial loss at epoch 0 is typically in the range of 10^{0}–10^{1} (roughly 1–10), since the randomly initialised network has no knowledge of the physics. By epoch 5,000, the loss should have fallen to approximately 10^{-3}–10^{-4}, a reduction of several orders of magnitude. (Exact values depend on the random seed and hardware.)
- What would happen if we increase the epochs from 5000 to 20,000 epochs? The loss would continue to decrease, potentially reaching 10^{-5} or lower, though with diminishing returns. The curve typically flattens (the network approaches a local or global minimum). The solution quality would improve marginally as the PINN curve would be closer to the numerical reference. For this relatively simple problem, 5,000 epochs already gives a good solution. There is a risk of overfitting to the collocation points at the expense of smoothness, though this is rare with physics losses.

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="Numerical Solver (Ground Truth)", linewidth=2.5)
ax1.plot(t_eval.cpu().numpy(), theta_pred.cpu().numpy(), 'k-',
         label="PINN Prediction", alpha=0.9)
ax1.set_title(f"Damped Pendulum  (b={b},  cp={cp:.2f})", 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(
    ic_loss_history,
    color='orange',
    linewidth=0.8,
    label="IC Loss")
ax2.plot(
    physics_loss_history,
    color='red',
    linewidth=1.2,
    alpha=0.8,
    label="Physics Loss")
ax2.set_yscale("log")
ax2.set_title("PINN Training Loss", fontsize=16)
ax2.set_xlabel("Epoch", fontsize=12)
ax2.set_ylabel("PINN Losses (log scale)", fontsize=12)
ax2.grid(True, which="both", linestyle="-", alpha=0.2)
ax2.legend(fontsize=13)

plt.tight_layout()
plt.show()

Why Use a PINN Instead of a Numerical ODE Solver?

For a simple, well-posed ODE like this damped pendulum, a classical numerical solver (odeint) is faster, more accurate, and perfectly sufficient. So why bother with a PINN at all?

The value of PINNs emerges when the problem moves beyond what classical solvers handle easily:

Scenario	Numerical Solver	PINN
Inverse problems: some parameters (e.g. b, c_p) are unknown and must be inferred from sparse observations	Requires a separate optimisation loop wrapped around the solver	The unknown parameters become additional trainable variables inside the same loss function; the physics loss and data loss are optimised jointly
Sparse or noisy data: only a few measured data points are available	Cannot incorporate scattered observations directly	Naturally blends a data-fitting loss with the physics loss, using the ODE as a regulariser to fill gaps between measurements
High-dimensional PDEs: e.g. Navier–Stokes in 3-D	Requires mesh generation and scales poorly with dimension (curse of dimensionality)	Mesh-free: collocation points are sampled randomly, and the neural network approximates the solution as a continuous function over the entire domain
Complex or coupled multi-physics systems	Each new coupling often needs a new specialised solver or re-meshing	The same framework (network + physics loss) extends to coupled systems by adding residual terms for each equation
Transfer learning and real-time prediction: need fast predictions for varying parameters.	Must re-solve from scratch for every new parameter set.	A trained PINN (or a meta-learned PINN) can generalise across parameter ranges, enabling near-instant inference once trained.
Data and model versioning: tracking how models and datasets evolve over time	Data and model versioning with numerical models are difficult due to large model size from meshing (especially with commercial numerical solver).	The trained network weights, architecture, and training data (collocation points, observations) form a self-contained artefact that can be versioned, checkpointed, and shared using standard ML tooling (e.g. MLflow, DVC).
Model storage: persisting a trained/solved model for later reuse	Saving numerical models requires large storage space (often gigabytes for fine meshes in 3-D), since the full mesh, field data, and solver state must all be stored.	Only the network weights need to be saved (typically kilobytes to a few megabytes). The same model can be fully reproduced later by reloading the weights into the architecture, without re-running the solver.

In short, simple forward ODE problems like this one serve as a pedagogical stepping stone. They let us validate that the PINN framework works correctly (by comparing against a trusted solver) before applying it to problems where classical methods struggle or fail entirely.

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.