Scientific Machine Learning:
Physics-informed Neural Networks (PINN)
A physics-informed neural network is a neural network trained to satisfy a differential equation. Instead of, or in addition to, fitting data, its loss measures how badly the network output violates the governing equation, the boundary conditions and the initial conditions, with every derivative supplied by automatic differentiation. This notebook explains why a purely data-driven model needs physics, how the physics enters the loss, and how the resulting network solves forward and inverse problems with little or no data.
Learning Objectives¶
After studying this notebook, you should be able to
- explain what a purely data-driven model gives up when data is scarce, place physics-informed learning on the two axes of how knowledge is represented and where it is injected, and state what the universal approximation theorem does and does not promise;
- write the physics-informed loss as a weighted sum of residuals of the governing equation at collocation points and of the initial and boundary conditions, and explain why the network must be differentiated with respect to its input;
- explain how observed data enter the same loss to turn an under-determined inverse problem into an over-determined one, and how unknown parameters become trainable tensors;
- implement a PINN by hand in PyTorch, with
requires_gradon the input coordinate andtorch.autograd.grad(..., create_graph = True)for the derivative, and predict the loss at initialisation from the equation before running anything; - solve a fourth-order beam equation and a two-dimensional plane-stress problem with a PINN, read the residual and boundary terms in the training trace, and explain which fields the unit weighting neglects and how measured data repair them;
- set up an inverse problem as the same optimization with a longer decision vector, and explain the roles of data and prior knowledge in making it well posed.
Overview
Consider the position of an engineer who wants to predict the temperature field inside a turbine blade. One option is to instrument blades, run them, record temperatures, and fit a model to the measurements. The other is to write down the heat equation, supply the material properties and the boundary conditions, and solve. The first approach needs data and no physics; the second needs physics and no data. Both have been available for decades, and the interesting question is what happens when the two are combined rather than chosen between.
This section motivates that combination. It begins with an honest account of what data-driven modelling does well, because the case for physics-informed learning is not that machine learning fails, but that it fails in specific, identifiable ways when data is thin. It then organizes the space of possible combinations along two axes, namely how domain knowledge is represented and where in the learning pipeline it is injected, and locates PINNs at a particular point in that space. Finally it examines the theoretical result that makes the enterprise plausible at all, the universal approximation theorem, and is careful about what the theorem does and does not promise.
Advantages of data-driven approaches
When sufficient data is available, a data-driven model can reach a high level of predictive accuracy without requiring extensive domain knowledge from its builder. This is a genuine and important strength. A convolutional network that classifies surface defects does not need a theory of how the defects form; it needs labelled photographs. The method scales with data rather than with understanding, and in domains where data accumulates automatically, that is an enormous practical advantage. Nothing in this notebook argues against using data when you have it. Section 6.3 in fact argues the opposite, and shows measurable gains from adding observations to a model that was already solving the correct equation.
The difficulty with black boxes
The trouble begins when the model is asked to do something the data did not cover. Four related failures recur often enough in engineering practice to be worth naming individually.
The first is physical inconsistency. A purely data-driven predictor has no notion of conservation, so nothing prevents its predicted velocity field from creating mass out of nothing, or its predicted stress field from violating equilibrium. The predictions may be numerically close to the training data and still be physically impossible, which in a design setting is worse than being merely inaccurate, because the error is invisible to any check the model itself can perform.
The second is spurious correlation. A model minimizes error on the data it is given, and any regularity present in that data is fair game, whether or not it reflects a mechanism. If every high-temperature run in the training set happened also to be a high-pressure run, the model has no way to tell which variable is doing the work, and it will confidently attribute the effect to the wrong cause. Explainable-AI methods are largely about detecting this situation after the fact; the physics-informed approach tries to prevent it in advance by constraining what the model is allowed to believe.
The third is poor generalization outside the training distribution. This is the same point stated in terms of inputs rather than variables. A model interpolates well and extrapolates badly, and the boundary between the two is precisely the edge of the region the training data covered. Engineering design very often asks exactly the extrapolation question: what happens if we make it thinner, run it hotter, spin it faster than anything we have measured.
The fourth is lack of interpretability. A prediction that comes with no mechanism attached cannot be argued with, cannot be checked against known limiting cases, and cannot contribute to scientific understanding. For a recommendation engine this may be acceptable. For a component whose failure is expensive, it usually is not, and understanding the mechanism that governs a system is a large part of what engineering is for.
What physics-informed learning proposes
Physics-informed neural networks respond to these limitations by integrating domain knowledge directly into the model's training. The claim is threefold. Combining data-driven flexibility with established physical law improves accuracy and robustness, because the physics constrains the model in exactly the regions where data is absent. Enforcing known constraints mitigates data scarcity and class imbalance, because the governing equation holds everywhere in the domain, including at points where nothing was measured, and can therefore supply supervision there. And requiring outputs to be consistent with fundamental principles improves interpretability, because a solution that satisfies the governing equation can be reasoned about with the whole apparatus of the underlying science.
By embedding prior knowledge from physics and engineering, PINNs aim at a framework that is simultaneously more generalizable, more interpretable, and more physically consistent than a purely data-driven one. The rest of this section makes that aim concrete, and Section 2 makes it computational.
Before building anything it is worth seeing the shape of the design space, because "add physics to the network" is not one idea but a family of them, and the members of the family differ in ways that matter. Two questions organize the family. First, in what form is the knowledge available? Second, at which point in the modelling pipeline is it inserted? The figure below arranges the possibilities along these two axes, with the network architectures listed on one side.
Deep neural networks
The left column of the figure lists architectures, most of which have appeared earlier in this course. Artificial neural networks are the fully connected building block of Artificial Neural Networks, and they are the architecture used in every laboratory here. Convolutional neural networks specialize to data laid out on a grid, recurrent neural networks handle sequences, graph neural networks operate on data whose natural structure is a graph, such as a molecular structure or a mesh, and generative models, from autoencoders to adversarial and diffusion models, produce new samples from a learned distribution.
It is worth noticing that the physics-informed idea is orthogonal to this column. Nothing in the construction of Section 2 depends on the network being fully connected. A physics-informed convolutional network and a physics-informed graph network are both sensible objects, and both exist. Every laboratory in this notebook uses a small fully connected network with $\tanh$ activations, and it will help to understand why: the loss involves second, third, and even fourth derivatives of the network output, and a $\tanh$ network is smooth to all orders while a ReLU network has a second derivative that is zero almost everywhere. A ReLU network cannot represent the curvature of a beam, because its curvature is identically zero wherever it is defined. This is the first of several places where the physics dictates a modelling choice that would be arbitrary in ordinary supervised learning.
Knowledge representation
The middle column asks what form the domain knowledge takes, because that determines whether it can be written as a differentiable expression. Differential equations describe continuous dynamical systems and are the representation used throughout this notebook. Algebraic equations express relationships and constraints that hold pointwise without derivatives, such as an equation of state. Knowledge graphs capture relational and semantic structure between entities. Simulation results are the structured output of a computational model, and they are knowledge in the sense that they encode the behaviour of a solver we trust. Human feedback covers expert judgement and manual annotation.
The representations differ in how directly they can be turned into a loss. A differential equation is nearly ideal for this purpose: it is an equation asserting that a certain expression equals zero, and the deviation from zero is immediately a number that can be squared and minimized. That is exactly why PINNs took the form they did.
Knowledge integration
The right column asks where the knowledge is inserted, and the three answers correspond to three quite different engineering practices. Feature engineering inserts knowledge into the inputs, by constructing quantities that the modeller knows to be relevant, such as a Reynolds number rather than separate velocity, length, and viscosity. Designing inserts knowledge into the architecture, by building the constraint into the structure of the model so that it cannot be violated. Regularizing inserts knowledge into the objective, by adding a penalty term that grows when the constraint is broken.
The standard PINN of Sections 2 through 6 is a regularizing method: the governing equation appears as an additive penalty in the loss, and the network is free to violate it, merely at a cost. The hard-constraint method, which Regularization, MAP Estimation, and Loss Weighting takes up, is a designing method: the boundary condition is built into the functional form of the solution, so it holds exactly and no penalty is required. Keeping these two categories distinct will make the hard-constraint construction read as a natural alternative rather than a trick. Embedding Physical Knowledge into Neural Networks takes up the same taxonomy from a broader angle.
Before we ask a neural network to represent the solution of a differential equation, we should ask whether it can represent that solution at all. The answer is supplied by a family of results collectively known as universal approximation theorems, the best known of which dates from 1989 and gives this subsection its title.
What the theorem says
Informally, a feedforward network with a single hidden layer, a finite number of neurons, and a suitable non-linear activation function can approximate any continuous function to arbitrary accuracy. Stated with the quantifiers in the right order, which is where the content lies:
Let $K \subset \mathbb{R}^n$ be compact, let $f : K \to \mathbb{R}$ be continuous, and let $\varepsilon > 0$ be given. Then there exist an integer $m$, weights $w_i \in \mathbb{R}^n$, biases $b_i \in \mathbb{R}$, and output coefficients $c_i \in \mathbb{R}$ such that the single-hidden-layer network
$$ N(x) = \sum_{i=1}^{m} c_i \, \sigma\!\left(w_i^{\top} x + b_i\right) $$
satisfies
$$ \sup_{x \in K} \left| f(x) - N(x) \right| < \varepsilon. $$
The 1989 result is stated for Borel measurable functions with respect to an appropriate notion of closeness, which is more general than the continuous case written above; the continuous version is the one worth memorizing. The activation $\sigma$ must be non-linear in an essential way. A polynomial activation fails, because a sum of polynomials of bounded degree is again a polynomial of bounded degree and cannot approximate an arbitrary continuous function. Every activation used in practice, including the $\tanh$ used throughout this notebook, satisfies the requirement.
Reading the quantifiers
The order is: given $f$ and given $\varepsilon$, there exists $m$. The width is chosen after the target function and the tolerance are known, and nothing in the statement bounds it. That is the first thing the theorem does not say.
Three things the theorem does not say
Students routinely overread this result, and the overreadings are worth naming, because each of them is a false expectation that will be disappointed in Section 4 and later.
First, the theorem is an existence statement, not a construction. It asserts that suitable weights exist; it gives no procedure for finding them. Training a network is a non-convex optimization problem, and gradient descent is guaranteed to find a stationary point of a smooth objective and nothing more. The gap between "a good network exists" and "training will find it" is not a technicality. It is where most of the practical difficulty in this notebook lives, and Regularization, MAP Estimation, and Loss Weighting is entirely devoted to it. When a training loop in a later section stops descending well short of zero, as the one in Section 4.1 does, no theorem has been violated. The approximating network exists; the optimizer simply did not reach it.
Second, the theorem says nothing useful about $m$. The number of hidden units required may grow very rapidly with the input dimension and with $1/\varepsilon$, and for a badly behaved $f$ it may be astronomically large. A guarantee that a network of unspecified size exists is not a guarantee that a network of practical size exists. In practice the depth and width are chosen empirically, which is why the laboratories in this notebook use networks of thirty-two, sixty-four, and two hundred and fifty-six units without deriving those numbers from anything.
Third, and most important for us, the theorem controls the value of the function, not its derivatives. The quantity bounded by $\varepsilon$ is $|f(x) - N(x)|$. It is a fact that two functions can be uniformly close while their derivatives are arbitrarily far apart. Take
$$ f_n(x) = \frac{\sin\!\left(n^2 x\right)}{n}. $$
Then $\sup_x |f_n(x)| = 1/n$, which tends to zero, so $f_n$ converges uniformly to the zero function. But $f_n'(x) = n \cos(n^2 x)$, so $\sup_x |f_n'(x)| = n$, which tends to infinity. A tiny, very rapidly oscillating wiggle is invisible in the function value and enormous in the derivative.
This matters here more than in any other notebook of the course, because a PINN's loss is built almost entirely from derivatives. In Section 5 we will differentiate the network four times. A network that matches the true deflection of a beam to within $10^{-4}$ everywhere may still have a completely wrong fourth derivative, so approximation in value alone would be no guarantee at all. Stronger versions of the approximation theorem do control derivatives, requiring correspondingly stronger smoothness of the activation function, and it is those versions rather than the elementary one that make the PINN construction theoretically sensible. This is a second reason, independent of the ReLU argument above, to prefer smooth activations here.
The practical lesson is worth stating once and remembering: representability and findability are different questions, and only the second one is decided by your training loop.
Neural networks have emerged as a flexible alternative to traditional numerical methods for differential equations. The appeal has three parts, and it is worth separating them from the beginning so that later claims can be judged against the right standard.
The first is that the representation is mesh-free. A finite element or finite difference method must first partition the domain, and for complicated geometry that partition is often the most laborious and least automatable step in the whole analysis. A neural network takes coordinates as input and returns a value; it needs points, but it does not need a mesh connecting them, and points can be scattered anywhere, including just inside an irregular boundary.
The second is that the solution obtained is a continuous, differentiable function rather than a table of nodal values. Once trained, the network can be evaluated at any point of the domain, and its derivatives at that point are available exactly by automatic differentiation rather than by finite differencing between neighbouring nodes.
The third is that the cost scales with the number of sampled points rather than with the number of grid cells, which for a $d$-dimensional domain grows like $N^d$. This is the argument that makes the approach attractive for high-dimensional and strongly non-linear problems, where conventional discretization becomes impractical.
Against these advantages stands one substantial disadvantage, which the laboratories will make plain: for a low-dimensional, well-posed forward problem with a regular geometry, a classical solver is faster, more accurate, and more reliable than a PINN. Section 3 states this openly, and it is the reason the second half of this notebook concentrates on inverse problems, where classical methods are far weaker and the comparison runs the other way.
The idea is older than the current wave of interest. Two foundational works illustrate the line of development, the first proposing a neural algorithm for solving differential equations in 1990.
The second, from 1998, developed artificial neural networks for both ordinary and partial differential equations, and already contained the essential ingredients: a network taking the independent variables as input, derivatives obtained by differentiating the network, and a cost function assembled from the residual of the equation. What the modern formulation adds is not a new principle but a practical substrate, namely automatic differentiation as a standard, efficient, and reliable service provided by every deep learning framework, together with the optimizers and hardware that make training a large network routine.
The following paper has played a significant role in advancing the field of physics-informed neural networks, and it is the paper that gave them their name. It is highly recommended for those interested in this emerging area of research:
Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations in Journal of Computational Physics (2019)
This groundbreaking work explores how deep learning can be used for both forward and inverse problems in non-linear partial differential equations. Notice that forward and inverse problems appear together in its title, and that the two are handled by the same machinery. In the forward problem the equation is known completely and the solution is wanted; in the inverse problem part of the equation is unknown and observations of the solution are available. It is a genuine feature of the PINN formulation, and the subject of Section 3 and of the inverse problems in Physics-Informed Neural Networks: Labs, that moving from one to the other requires no new algorithm, only a longer list of trainable variables.
The following review surveys the prevailing trends in embedding physics into machine learning, and describes current capabilities, limitations, and applications of physics-informed learning, including forward and inverse problems, the discovery of hidden physical laws, and high-dimensional problems.
- Physics-informed machine learning in Nature Reviews Physics (2021)
If you read only one thing beyond this notebook, the review is the better choice, because it places PINNs within the wider family of methods sketched in Section 1.1 rather than treating them as the whole subject.
Overview
We now build the object. The construction is short enough to state in one paragraph, and the rest of the section is spent making each piece of it precise.
Let a neural network $u_\theta(x,t)$ take the independent variables of the problem as its inputs and return the unknown field as its output. Because the network is a composition of affine maps and smooth activations, it is a smooth function of $x$ and $t$, and its partial derivatives of any order can be computed exactly by automatic differentiation. Substitute those derivatives into the governing equation. The equation will not in general hold, and the amount by which it fails is a number, the residual, computable at any point of the domain we like. Square it, average it over a set of sampled points, add analogous terms measuring the violation of the initial and boundary conditions, and minimize the total over $\theta$. If the total reaches zero, the network satisfies the equation and its side conditions, and is therefore the solution.
Two features of this construction deserve emphasis before the details. The first is that no solved examples appear anywhere in it. Ordinary supervised learning requires pairs $(x_i, u_i)$ in which $u_i$ is the correct answer, obtained somehow. Here the correct answer is never known and never needed: the governing equation itself tells the network how wrong it is. The second is that the physics enters as an additive term in the objective function, not as an architectural restriction and not as a preprocessing of data. Regularization, MAP Estimation, and Loss Weighting will show the alternative, and the contrast is instructive, but the standard construction is the one described here.
Physics-informed neural networks integrate the laws of physics, represented by partial differential equations, along with initial and boundary conditions, into the training process of a neural network. This is what makes them powerful for solving forward and inverse problems even with limited or noisy data.
Key formulation
Consider a generic time-dependent partial differential equation defined on a domain $\Omega \times [0,T]$:
$$ \mathcal{N}[u(x, t)] = f(x, t), \quad (x, t) \in \Omega \times [0, T], $$
in which $\mathcal{N}$ is a differential operator involving derivatives of $u(x,t)$, the function $u(x,t)$ is the unknown solution to be approximated, $f(x,t)$ is a source or forcing term, $\Omega$ is the spatial domain, and $[0,T]$ is the temporal interval. The operator notation is deliberately abstract because it has to cover every equation in this notebook. For the beam of Section 5, $\mathcal{N}[y] = \partial^4 y / \partial x^4$ and there is no time variable at all. For steady heat conduction, $\mathcal{N}[T] = \kappa(\partial^2 T/\partial x^2 + \partial^2 T/\partial y^2)$ and $f = 0$. For a flow past a cylinder, $\mathcal{N}$ is the non-linear Navier-Stokes operator, and it acts on three coupled fields at once. The construction below does not care which.
Initial and boundary conditions
A differential equation on its own does not determine a solution. The equation constrains how the field changes from point to point; it says nothing about where the field starts. Two families of side conditions supply what is missing.
(1) The initial condition specifies the solution at the start of the time domain:
$$ u(x, 0) = g(x), \quad x \in \Omega. $$
(2) The boundary conditions specify the behaviour of the solution at the edges of the spatial domain. Two kinds appear repeatedly. A Dirichlet condition prescribes the value of the field on the boundary,
$$ u(x_b, t) = h_1(x_b, t), \quad x_b \in \partial \Omega, \, t \in [0, T], $$
and a Neumann condition prescribes its normal derivative there,
$$ \frac{\partial u}{\partial n}(x_b, t) = h_2(x_b, t), \quad x_b \in \partial \Omega, \, t \in [0, T], $$
where $n$ denotes the outward normal to the boundary. The distinction is physical, not merely mathematical. A Dirichlet condition on temperature is a wall held at a fixed temperature by contact with a reservoir; a Neumann condition on temperature is a prescribed heat flux, and the homogeneous case $\partial T/\partial n = 0$ is a perfectly insulated wall. For the beam of Section 5 the two types appear together at a clamped end, where both the deflection and the slope are prescribed. For the plate of Section 6 the free edges are Neumann-type conditions on stress, which is itself built from derivatives of the displacement.
Together, the equation, the initial condition, and the boundary conditions constitute a well-posed problem: there is exactly one $u$ satisfying all of them. That phrase exactly one is what makes the forward problem solvable with no data. When one of these ingredients goes missing, as it does when a boundary condition is unknown, uniqueness goes with it, and data must be brought in to restore it.
From the equation to a loss
We approximate the solution $u(x,t)$ by a neural network $u_\theta(x,t)$ parameterized by $\theta$, the collection of all weights and biases. Before writing the loss, it is worth being completely explicit about the computation that makes the whole scheme possible, because it is the one step that is genuinely new relative to every earlier notebook.
Two different derivatives
Training any network requires $\partial L / \partial \theta$, the gradient of the loss with respect to the weights, and that is what backpropagation computes. A PINN additionally requires $\partial u_\theta / \partial x$, $\partial^2 u_\theta / \partial x^2$, and so on: derivatives of the network's output with respect to its input. These are different objects, but they are computed by the same machinery, because automatic differentiation does not care which node of the computational graph you differentiate with respect to. If the framework can propagate a derivative back to a weight, it can propagate it back to an input coordinate just as easily.
Three properties of these input-derivatives matter.
They are exact. This is not finite differencing; there is no step size $h$, no truncation error of order $h^2$, and no catastrophic cancellation when $h$ is made small. Automatic differentiation applies the chain rule to the elementary operations that make up the network, and the result is the derivative to machine precision. For a fourth derivative, as in Section 5, this distinction is decisive: a fourth-order finite difference formula divides by $h^4$, and rounding error is amplified accordingly.
They are available at any point, not only at grid nodes. There is no mesh, so there is no interpolation and no question of where a derivative is defined. A point in the domain is simply a valid input to the network.
They cost extra. Each additional order of differentiation builds a further layer onto the computational graph, so a fourth-order equation is materially more expensive per training step than a first-order one, and the memory required to hold the graph grows correspondingly. The differentiation must also itself be recorded, since the loss built from a derivative has in turn to be differentiated with respect to the weights; this second-order use of automatic differentiation is what every laboratory in this notebook depends on.
The residual
Given the network, define the residual of the equation at a point $(x,t)$ as
$$ r_\theta(x,t) = \mathcal{N}[u_\theta(x,t)] - f(x,t). $$
The residual is a computable number for any $\theta$ and any point. It is zero exactly where the network happens to satisfy the equation, and its magnitude measures the local failure of the physics. This is the object that replaces the label in ordinary supervised learning: instead of comparing $u_\theta(x_i)$ against a known $u_i$, we compare $\mathcal{N}[u_\theta]$ against the known $f$. The known quantity is the forcing term, which we were given, rather than the solution, which we were not.
The four loss terms
The total loss is assembled from the residual and the side conditions, each averaged over its own set of points.
(1) Physics equation:
$$ \mathcal{L}_{\text{physics}} = \frac{1}{N_f} \sum_{i=1}^{N_f} \left| \mathcal{N}[u_\theta(x_i, t_i)] - f(x_i, t_i) \right|^2, $$
where $\{(x_i,t_i)\}_{i=1}^{N_f}$ are collocation points sampled from the interior $\Omega \times [0,T]$.
(2) Initial condition:
$$ \mathcal{L}_{\text{initial}} = \frac{1}{N_I} \sum_{i=1}^{N_I} \left| u_\theta(x_i, 0) - g(x_i) \right|^2, $$
where $\{x_i\}_{i=1}^{N_I}$ are points sampled from $\Omega$ at $t = 0$.
(3) Boundary conditions. For Dirichlet conditions,
$$ \mathcal{L}_{\text{Dirichlet}} = \frac{1}{N_B} \sum_{i=1}^{N_B} \left| u_\theta(x_b, t_i) - h_1(x_b, t_i) \right|^2, $$
with $\{(x_b, t_i)\}_{i=1}^{N_B}$ sampled on the boundary $\partial \Omega$, and for Neumann conditions,
$$ \mathcal{L}_{\text{Neumann}} = \frac{1}{N_B} \sum_{i=1}^{N_B} \left| \frac{\partial u_\theta}{\partial n}(x_b, t_i) - h_2(x_b, t_i) \right|^2. $$
Notice that the Neumann term is structurally identical to the physics term: both involve derivatives of the network, both are obtained by automatic differentiation, and neither requires knowing the solution. Only the set of points differs. This uniformity is why the implementations that follow look so repetitive; every term in every laboratory is the mean square of some expression built from $u_\theta$ and its derivatives, evaluated on some set of points.
(4) The total loss combines the components:
$$ \mathcal{L}_{\text{total}} = \mathcal{L}_{\text{physics}} + \lambda_I \mathcal{L}_{\text{initial}} + \lambda_B \mathcal{L}_{\text{boundary}}, $$
where $\lambda_I$ and $\lambda_B$ weight the contributions of the initial and boundary terms.
The weights are a modelling decision
It is tempting to read the $\lambda$'s as tuning knobs of no conceptual importance and set them all to one, which is what most laboratories in this notebook in fact do. That is a defensible default and it is not a neutral one, and the reasons are worth spelling out.
The first reason is that the terms are not measured in the same units. For the beam of Section 5 the physics residual has the units of $\partial^4 y / \partial x^4$, that is, length divided by length to the fourth power, while the Dirichlet residual $y(0)$ has the units of length. Adding them is adding a curvature-like quantity to a displacement, and the sum has no dimensional meaning unless $\lambda_B$ carries units to fix it. Setting $\lambda_B = 1$ is therefore a statement about the numerical scale of the problem rather than about its physics, and it is only reasonable because the problem has been nondimensionalized so that all quantities are of order one. This is why the laboratories go to the trouble of nondimensionalizing, and it is a good habit to import into your own problems: choose units so that every residual is of comparable size, and the default weights become defensible.
The second reason is that the sum is a scalarization of a multi-objective problem, in the usual sense of multi-objective optimization. There are three things we want to be small, they compete, and by adding them with weights we have chosen a single point on the trade-off surface. A different $\lambda_B$ selects a different point. Nothing in the mathematics privileges $\lambda_B = 1$; what privileges it is that the residuals happen to have comparable magnitudes.
The third reason is that the failure modes at the two extremes are different and both are recognizable in practice. If $\lambda_B$ is far too small, the network will satisfy the equation beautifully in the interior and drift away from the correct boundary values, and since a differential equation has infinitely many solutions before the boundary conditions are imposed, the result is a perfectly valid solution of the wrong problem. This is the more dangerous failure, because the physics loss will look excellent. If $\lambda_B$ is far too large, the network will nail the boundary and ignore the interior, producing a field that is correct on a set of measure zero and meaningless elsewhere. Regularization, MAP Estimation, and Loss Weighting returns to the practical management of this balance and shows what the weights actually decide.
By encoding physical laws and constraints directly into the learning process, PINNs offer a robust and data-efficient framework for problems in which measurements are expensive and equations are not.
Collocation points are not data
As you may have noticed, no data is required to find the solution of a partial differential equation with a PINN. Since the equation, the initial condition, and the boundary condition are all available, the solution is uniquely determined, and measurements would be redundant. What the training process does require is a mechanism for evaluating the loss and hence for updating the weights, and that mechanism needs points at which to evaluate it. To avoid confusion we call these collocation points rather than data. They are usually sampled at random from the input space, or laid out on a regular lattice, and the crucial fact about them is that a collocation point carries no label. It is an input coordinate and nothing more; the loss at that point is computed from the network's own output and its own derivatives, compared against the equation, and no knowledge of the true solution enters anywhere.
The word collocation is borrowed from classical numerical analysis, where a collocation method is one that enforces an equation exactly at a chosen finite set of points rather than in an integrated or weighted sense. That is precisely what the physics loss does, up to the fact that the enforcement is approximate and by penalty.
If observational data does happen to be available, it can be added as a further term of the same shape, and doing so improves accuracy. That is the subject of Section 6.3 and, in a more essential role, of the inverse problems in Physics-Informed Neural Networks: Labs. But it is an option, not a requirement, and the distinction between "cannot train without labels" and "can use labels if offered" is the practical difference between conventional deep learning and this method.
A common misconception
Students sometimes describe a PINN as a network trained on solutions of the equation, perhaps generated by a conventional solver, with the equation used as a regularizer. That describes a different and also useful method, and it is closer to what Deep Neural Operators and Fourier Neural Operator will do, but it is not what is happening here. In a PINN there is no solver, there are no reference solutions, and the equation is not a regularizer on a data-fitting objective; it is the entire objective. The exact solutions that appear in the laboratories of Sections 4 and 5 are used only after training, to draw a comparison curve, and are never seen by the optimizer.
Overview
We have now seen PINNs as a neural network approach to forward problems, and it is worth stating plainly what the last two sections did not claim. For a well-posed forward problem on a simple geometry, the method offers no decisive advantage over conventional computational techniques. A finite element solver for the beam of Section 5 would return the answer in milliseconds, to more digits, with an error estimate and a convergence theory behind it, where the PINN needs a thousand gradient steps and gives none of those guarantees. Section 5 is worth doing because it teaches the mechanics on a problem whose exact solution we can write down, not because it is the way one should analyse a beam.
The application of PINNs to inverse problems is a different matter, and it is the direction in which the method genuinely advances scientific computing. Inverse problems involve determining unknown parameters, unknown functions, or unknown physical laws from limited information. They are typically ill-posed, meaning that the available information does not determine a unique and stable answer, and this is not a defect of the solution method but a property of the problem. To obtain a usable answer, a combination of observed data and physical laws is essential, and the PINN framework happens to be an unusually convenient place to combine them.
One fact from linear algebra is worth having in mind for the whole of this section, and it can be stated in three sentences. A square non-singular system $Ax=b$ has exactly one solution, and it corresponds to the well-posed forward problem. An under-determined system, with fewer equations than unknowns, has infinitely many solutions, and an algorithm handed one will return a particular member of the family, typically the one of smallest norm, with no indication that others exist; that is the ill-posed inverse problem with insufficient data. An over-determined system, with more equations than unknowns, has in general no exact solution and is handled by least squares, returning the point that comes closest to satisfying everything at once. Adding observations to an ill-posed inverse problem moves it from the second situation towards the third, and the PINN loss, being a sum of squared residuals, is performing least squares in a non-linear setting. Section 5 in Physics-Informed Neural Networks: Labs stages that progression deliberately.
PINNs with data therefore provide a framework for inverse problems that integrates data-driven evidence with physics-based constraints derived from the governing equations, and the hybrid is more stable and more accurate than either ingredient alone.
Let us examine how PINNs estimate an unknown parameter, as the simplest example of solving an inverse problem.
Formulation of an inverse problem with PINNs
Consider a system governed by a partial differential equation in which the operator depends on a parameter $\lambda$:
$$ \mathcal{N}[u(x, t; \lambda)] = f(x, t), \quad (x, t) \in \Omega \times [0, T], $$
where $\mathcal{N}$ acts on the unknown solution $u(x,t)$, the quantity $\lambda$ is an unknown parameter or function to be inferred, $f(x,t)$ is a known forcing term, $\Omega$ is the spatial domain and $[0,T]$ the temporal interval. In the flow laboratory of Physics-Informed Neural Networks: Labs the role of $\lambda$ is played by the density and the viscosity of a fluid, both unknown; in its heat laboratory the unknown is not a number at all but an entire missing boundary condition.
The goal is to estimate $\lambda$ given a set of observed data points
$$ \{ (x_i, t_i, u_{\text{obs}}(x_i, t_i) ) \}_{i=1}^{N_d}. $$
These are data in the full sense of the word: input coordinates paired with measured values of the solution. They are what a sensor produces, and unlike collocation points they are expensive, few, and noisy.
PINNs solve this problem by learning the solution $u_\theta(x,t)$ and the unknown parameter $\lambda$ simultaneously. That word simultaneously is where the elegance lies, and it is worth pausing on before the loss is written. In a classical treatment, parameter identification is an outer optimization wrapped around an inner solve: guess $\lambda$, solve the equation for $u$, compare $u$ against the measurements, adjust $\lambda$, and repeat, with each iteration paying for a complete solution of the equation. In the PINN treatment there is no inner solve, because $u_\theta$ and $\lambda$ are found by the same gradient descent at the same time.
The same optimization with a larger decision vector
Here is the point stated as precisely as possible. In the forward problem, the decision variables of the optimization are the network parameters, and we minimize over $\theta \in \mathbb{R}^{p}$. In the inverse problem, the unknown $\lambda$ is declared to be a trainable variable as well, and we minimize over the concatenated vector $(\theta, \lambda) \in \mathbb{R}^{p+1}$. The objective is a differentiable function of $\lambda$ because $\lambda$ appears in the residual, so $\partial \mathcal{L} / \partial \lambda$ is available from the same automatic differentiation pass that produces $\partial \mathcal{L}/\partial \theta$, and the same optimizer updates both.
Nothing else changes. There is no separate estimation algorithm, no outer loop, no sensitivity analysis to be derived by hand. The entire modification is to create the parameter as a trainable tensor and to hand it to the optimizer alongside the weights, which is literally what that flow laboratory does in two lines. An inverse problem is the same optimization problem with a longer decision vector. This is the observation that makes PINNs attractive for identification, and it generalizes without effort: ten unknown parameters means ten extra components, and an unknown field means an extra network or an extra structural constraint.
Loss function for the inverse problem
To estimate $\lambda$ we define a composite loss with four parts.
(1) Physics residual loss, ensuring the equation holds at the collocation points:
$$ \mathcal{L}_{\text{physics}} = \frac{1}{N_f} \sum_{i=1}^{N_f} \left| \mathcal{N}[u_\theta(x_i, t_i; \lambda)] - f(x_i, t_i) \right|^2. $$
(2) Initial condition loss, enforcing the initial state of the system:
$$ \mathcal{L}_{\text{initial}} = \frac{1}{N_I} \sum_{i=1}^{N_I} \left| u_\theta(x_i, 0) - g(x_i) \right|^2. $$
(3) Boundary condition loss:
$$ \mathcal{L}_{\text{boundary}} = \frac{1}{N_B} \sum_{i=1}^{N_B} \left| u_\theta(x_b, t_i) - h_1(x_b, t_i) \right|^2. $$
(4) Data loss, ensuring the predicted solution matches the observations:
$$ \mathcal{L}_{\text{data}} = \frac{1}{N_d} \sum_{i=1}^{N_d} \left| u_\theta(x_i, t_i) - u_{\text{obs}}(x_i, t_i) \right|^2. $$
(5) Total loss:
$$ \mathcal{L}_{\text{total}} = \mathcal{L}_{\text{physics}} + \lambda_I \mathcal{L}_{\text{initial}} + \lambda_B \mathcal{L}_{\text{boundary}} + \lambda_D \mathcal{L}_{\text{data}}. $$
where $\lambda_I$, $\lambda_B$ and $\lambda_D$ balance the contributions. Only the fourth term is new, and it is the one that would be the entire loss in a conventional supervised model. A PINN for an inverse problem is thus literally a supervised regression on the sensor readings, plus penalties that force the fitted field to obey known physics everywhere else.
Notice how the two kinds of information cooperate. The data alone are far too few to determine a field over the whole domain; seventeen scattered velocity readings would let a purely data-driven model interpolate almost anything between them. The physics alone cannot determine the field either, because the parameter it depends on is unknown. Together they can: the equation propagates the information contained in the sparse measurements across the entire domain, and the measurements pin down the parameter that the equation leaves free. Neither ingredient is sufficient and each supplies exactly what the other lacks.
Why observations restore well-posedness
Return to the picture drawn above. Without data, an inverse problem is under-determined: many pairs $(u, \lambda)$ satisfy the equation and its side conditions, so a PINN trained on physics alone will drive its loss to a small value and report one of those pairs, with no indication that others exist. This is not a training failure and no amount of extra epochs will fix it. Adding $N_d$ observations adds $N_d$ further conditions on the same unknowns; when there are enough of them the system becomes over-determined, no exact solution exists, and minimizing the sum of squares returns the best compromise, which is stable and, if the measurements are informative, close to the truth.
Summary of the idea
A PINN alone cannot approximate the solution of an ill-posed problem. If observational data is available, however, it can be used to construct an over-determined system, and that is what makes the inverse problem tractable. Physics-Informed Neural Networks: Labs carries this out twice, once for unknown scalar parameters and once for an unknown boundary condition, and additionally shows what happens when prior structural knowledge, namely that the missing boundary temperature is constant, is added on top of the data.
Overview
We have discussed the fundamental theory behind physics-informed neural networks; now we build one. Section 4.1 implements a PINN by hand, with every derivative, every residual, and every gradient step written out explicitly.
Doing it by hand first is not ceremony. A PINN implementation contains exactly one unfamiliar idea, namely differentiating the network output with respect to the network input, and once you have written those three lines yourself the rest of the notebook reads as variations on a theme. Every laboratory from Section 5 onward has the same skeleton: a helper that computes a derivative, a function that assembles the residual of the governing equation, one function per boundary condition, and a loop that squares everything, adds it up, and calls the optimizer. Learn the skeleton here and the later laboratories become exercises in translating physics rather than exercises in programming.
There is a second reason to write the loop out by hand. In a framework that hides training behind a single fit call, the three operations that actually constitute gradient descent are invisible; PyTorch never hides them, and the loop below is the same zero_grad, backward, step triple as in every earlier lab.
Here they appear explicitly: clear the accumulated gradient, evaluate it, move against it. In an ordinary supervised problem this transparency is a convenience. In a PINN it is close to a necessity, because the loss itself contains a derivative, so the backward pass is being asked to differentiate through a differentiation, and it helps a great deal to be able to see where each of the two differentiations happens.
Note. The first two laboratories train in seconds and are run in place. The thin plate of Section 6 trains for ten thousand epochs on 2,500 points, which takes minutes on a GPU and much longer without one, so its two training cells may be skipped: the models they produce, together with the finite element reference data, are downloaded from Google Drive in the cells that follow, and every result of Section 6 is drawn from those files.
We begin with about the simplest differential equation that is not trivial: a first-order ordinary differential equation whose exact solution can be written down in one line, so that every claim about the network can be checked against the truth. This first laboratory is implemented from scratch in PyTorch.
- Let's look at the ODE:
$$\frac{du}{dt} = \cos2\pi t$$
- Initial condition:
$$u(0) = 1$$
- The exact solution:
$$u(t) = \frac{1}{2\pi}\sin2\pi t + 1$$
Where the exact solution comes from
The problem is simple enough to solve by inspection, but writing out the two steps is worth the space because it fixes what the network is being asked to reproduce. Integrating both sides of the equation with respect to $t$,
$$ u(t) = \int \cos 2\pi t \; dt = \frac{1}{2\pi}\sin 2\pi t + C, $$
where the factor $1/2\pi$ appears because differentiating $\sin 2\pi t$ produces a factor $2\pi$ by the chain rule, which must be cancelled. The constant $C$ is fixed by the initial condition: setting $t=0$ gives $u(0) = \frac{1}{2\pi}\sin 0 + C = C$, and since $u(0)=1$ we have $C = 1$, giving the exact solution stated above.
Two features of this solution are worth noticing in advance, because they explain design decisions that follow. The solution oscillates with period $1$, and the training domain will be $t \in [0,2]$, so the network must represent two full periods. And the solution takes values between $1 - 1/2\pi \approx 0.841$ and $1 + 1/2\pi \approx 1.159$; it is centred on $1$, not on $0$, and its total variation is small compared with its mean. Predicting a nearly constant function with a small ripple on top is harder than it looks, because the loss is built from the derivative, where the ripple is everything and the constant contributes nothing.
Make a neural network and loss functions like below:
Import Library
The implementation uses PyTorch, whose automatic differentiation records operations as they are executed, so that gradients can be taken inside an ordinary Python loop with no separate graph-construction step, together with NumPy and Matplotlib for the arithmetic and the plotting.
import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
%matplotlib inline
plt.rcParams['figure.dpi'] = 72
np.random.seed(0)
torch.manual_seed(0)
Define Network and Hyper-parameters
The network takes one input, $t$, and returns one output, the predicted $u$. Between them are two hidden layers of thirty-two units with $\tanh$ activations. Before running the cell, work out the parameter count for yourself: the first layer holds $1 \times 32$ weights and $32$ biases, giving $64$; the second holds $32\times32 + 32 = 1056$; the output layer holds $32 \times 1 + 1 = 33$. The total is $1153$, which is what the printed count reports.
Two choices here are physics-driven rather than conventional, and both are easy to pass over.
The activation is $\tanh$ rather than ReLU. The loss will contain $du_\theta/dt$, so the network must be differentiable, and Section 5 will need four derivatives, for which a ReLU network is useless: it is piecewise linear, so its second derivative vanishes wherever it exists. A $\tanh$ network is smooth to all orders, which is exactly the property the method requires.
The output layer has no activation at all. This is deliberate and not merely conventional: the target values lie near $1$, and a $\tanh$ output would be confined to $(-1,1)$, making the correct answer unreachable at every point where $u > 1$, which is half the domain. A linear output layer places no restriction on the range. Keep this observation in mind when you reach Section 5, where the beam network does end with a $\tanh$ and gets away with it only because the true deflection never leaves $[-0.125, 0]$.
The optimizer is Adam at a learning rate of $0.001$, the default throughout this notebook. From the optimizer's point of view this is an ordinary differentiable objective in $1153$ variables; nothing about the physics-informed setting changes how it works. One detail of its construction matters and matters again in every inverse problem: the optimizer holds a fixed list of the tensors it is allowed to update, here NN.parameters(), and anything not in that list is not trained however many gradients flow through it. Appending an unknown physical constant to that list is the whole difference between a forward and an inverse problem in code.
NN = nn.Sequential(
nn.Linear(1, 32),
nn.Tanh(),
nn.Linear(32, 32),
nn.Tanh(),
nn.Linear(32, 1)
)
print(NN)
print('parameters:', sum(p.numel() for p in NN.parameters()))
optm = torch.optim.Adam(NN.parameters(), lr = 0.001)
Define ODE System
Now the loss. Following Section 2, it has two terms.
- ODE loss, the mean squared residual over the collocation points:
$$L_{\text{ODE}} = \frac{1}{N_f}\sum_{i=1}^{N_f} \left(\frac{d \text{u}(t_i)}{dt}- \cos2\pi t_i\right)^2$$
- Initial condition loss:
$$L_{IC} = \frac{1}{N_I}\sum_{i=1}^{N_I} \left({\text{u}(0)}- 1\right)^2$$
- Total loss:
$$L_{\text{Total}} = L_{\text{ODE}} + L_{\text{IC}}$$
The two terms are added with equal weight, which in the notation of Section 2 means $\lambda_I = 1$. Since $u$ is of order one and $du/dt$ is of order one, the two residuals have comparable magnitudes and the choice is reasonable.
The cell that follows implements this, and one idea in it carries the whole method: the input coordinate $t$ is created with requires_grad = True, so that a graph connecting $t$ to the network output exists and can be walked backwards, and torch.autograd.grad(u, t, ...) walks it. Two arguments of that call matter. grad_outputs = torch.ones_like(u) asks for the derivative of each output with respect to its own input, which is what a pointwise network gives, and create_graph = True records the differentiation itself so that the loss built from $du_\theta/dt$ can be differentiated again with respect to the weights. Ordinarily one tracks the weights, because the weights are what will be updated. Here the input is tracked as well, not because $t$ will ever be updated, but because we want $du_\theta/dt$. This is the single most common omission in a first hand-written PINN, and it produces not a wrong answer but no answer at all.
The rest of the function is arithmetic. The equation residual $du_\theta/dt - \cos 2\pi t$ is evaluated at all twenty collocation points at once and its mean square is $L_{\text{ODE}}$; the initial-condition error $u_\theta(0) - 1$ is a single number, since that condition is enforced at the one point $t=0$, and its square is $L_{IC}$. The two are added with no explicit $\lambda_I$, which is the equal weighting stated above.
def ode_system(t, net):
t = torch.tensor(t, dtype = torch.float32, requires_grad = True)
u = net(t)
u_t = torch.autograd.grad(u, t, grad_outputs = torch.ones_like(u), create_graph = True)[0]
ode_loss = torch.mean((u_t - torch.cos(2*np.pi*t))**2) # ODE residual at the collocation points
IC_loss = (net(torch.zeros(1, 1)) - 1)**2 # initial condition at t = 0
return ode_loss + IC_loss
Train & Prediction
The training loop is ordinary gradient descent with one unusual feature: the collocation points are drawn afresh at every iteration, twenty of them uniformly distributed on $[0,2]$. The network therefore never sees the same set twice, and over three thousand iterations the domain is covered by sixty thousand distinct points. This is a cheap way of ensuring that the residual is small everywhere rather than at a fixed set of locations where the network could overfit, and it has a price that will become visible in the plotted loss curve. PINN: Limitations and Extensions discusses this and related sampling strategies at length.
Two nested differentiations happen at every step, and separating them explicitly is worth the effort, since this is the one structural feature that distinguishes a PINN's training step from an ordinary one. The inner differentiation produces $du_\theta/dt$, the derivative of the output with respect to the input. The outer differentiation produces $\partial L/\partial\theta$, the derivative of the loss with respect to the weights, and because the loss contains $du_\theta/dt$, the second differentiation must pass through the first. That is second-order automatic differentiation, and it is what makes a derivative-based loss trainable at all; in PyTorch it is exactly what create_graph = True buys, and without it the inner derivative would be a constant to the outer backward().
Predict the first printed number
Before running the cell, work out what the loss should be at iteration zero. The weights are freshly initialized, so the network output is close to zero everywhere and so is its derivative. With $u_\theta \approx 0$ we get $\text{IC\_loss} \approx 0 - 1 = -1$, contributing $1$ to the total. With $du_\theta/dt \approx 0$ the ODE residual is approximately $-\cos 2\pi t$, and the average of $\cos^2$ over a whole number of periods is $1/2$, contributing about $0.5$. The predicted initial loss is therefore about $1.5$.
train_loss_record = []
for itr in range(3000):
train_t = np.random.rand(20, 1)*2
optm.zero_grad()
train_loss = ode_system(train_t, NN)
train_loss.backward()
optm.step()
train_loss_record.append(train_loss.item())
if itr % 1000 == 0:
print('Epoch: {} Loss: {:.4f}'.format(itr, train_loss.item()))
plt.figure(figsize = (6, 4))
plt.plot(train_loss_record)
plt.show()
The printed values are $1.3936$ at iteration $0$, $0.3519$ at iteration $1000$, and $0.0004$ at iteration $2000$. The opening value confirms the prediction to two figures, which is the useful part: the arithmetic above was a genuine test of whether we understood what the loss measures, and it passed.
The rest of the trace is the interesting part, and the plotted curve carries information the three printed numbers do not. For roughly the first fifteen hundred iterations the loss wanders noisily between about $0.2$ and $0.6$, declining but nowhere near zero; the value of $0.3519$ at iteration $1000$ means that the mean squared equation residual is still of the order it started from, and that the network has learned the level of the solution far better than its oscillation. Then, over a few hundred iterations, the loss drops by three orders of magnitude and settles near $4\times10^{-4}$, where it stays; the value of $0.0004$ at iteration $2000$ is the settled level. That shape, a long plateau followed by a sudden descent, is typical of PINNs on oscillatory solutions. The plateau is the network fitting the mean and the overall bend of $u$; the descent is the moment at which its derivative locks onto $\cos 2\pi t$. Nothing in the loss told us in advance when, or whether, that moment would come.
The noise on the curve has an identifiable source: the collocation points are resampled at every iteration, so the quantity being plotted is not one fixed function evaluated at successive iterates but a fresh random estimate at each step. Twenty points is a small sample from which to estimate a mean square, and an unlucky draw that happens to concentrate near the steep parts of $\cos 2\pi t$ will report a larger loss than the previous iterate even though the weights improved. Under a fixed learning rate, a stochastic gradient method does not converge to a minimizer at all; it converges to a neighbourhood whose size is set by the product of the step size and the gradient noise, and then wanders inside it, which is what the flat tail of the curve shows. When you read a PINN loss curve, ask first whether the objective is fixed or resampled, because the two produce very different-looking plots for the same underlying progress.
One more fact belongs here. The same six cells, run in the original version of this lab under a different framework whose default initialisation differs from PyTorch's, printed $1.4421$, $0.4021$ and $0.3037$ at the three checkpoints and never left the plateau in three thousand iterations. Same equation, same architecture, same optimizer, same step size, same iteration count; a different random starting point; a factor of a thousand in the final loss. Nothing here contradicts Section 1.2, and it is worth stating the relationship in the direction that actually holds. The universal approximation theorem guarantees that a network of this shape exists which makes this loss essentially zero. Whether a particular run finds it is a fact about the optimization, not about the representation, and it is exactly the gap between representability and findability that Section 1.2 warned about. Regularization, MAP Estimation, and Loss Weighting and PINN: Limitations and Extensions are devoted to closing it.
test_t = np.linspace(0, 2, 100).reshape(-1, 1)
true_u = np.sin(2*np.pi*test_t)/(2*np.pi) + 1
with torch.no_grad():
pred_u = NN(torch.tensor(test_t, dtype = torch.float32)).numpy()
plt.figure(figsize = (6, 4))
plt.plot(test_t, true_u, 'k', label = 'True', linewidth = 3, alpha = 0.3)
plt.plot(test_t, pred_u, 'r--', label = 'Prediction')
plt.legend()
plt.xlabel('t')
plt.ylabel('u')
plt.show()
The comparison plot places the exact solution and the network prediction on the same axes, and it should be read against the loss just discussed rather than in isolation.
The relationship between a loss value and the corresponding picture is quantitative and runs in a direction students often get backwards. An error of amplitude $\varepsilon$ in a solution oscillating at angular frequency $2\pi$ produces an error of amplitude $2\pi\varepsilon$ in the derivative, so a derivative-based loss magnifies displacement errors by a factor of about six here. A PINN loss therefore looks worse than the corresponding picture, systematically and by a knowable factor, and a loss that has not reached zero can still correspond to a plot in which the two curves are close. The converse is the more important half: a plot in which the curves lie on top of one another is weak evidence, since it is blind to exactly the errors the loss is most sensitive to. Judge a PINN by its residual, and use the plot to see where the residual lives.
That is the entire method, implemented in six cells. Everything that follows adds physics, not machinery.
Overview
Having implemented a PINN by hand in PyTorch, we now extend our study to a more practical and familiar problem: the Euler beam bending problem, a fundamental topic in solid mechanics. This implementation uses PyTorch, and every laboratory from here to the end of the notebook is written in it directly, without a library in between.
The Euler beam is a good next step for three reasons. It is a genuine engineering problem with a hundred and fifty years of history and a closed-form answer, so every number the network produces can be checked. It is fourth order, which stretches the automatic differentiation machinery in a way the first-order problem of Section 4 did not, and which makes concrete the remark in Section 2 that each additional derivative costs something. And its boundary conditions are of two different kinds at two different places, so the composite loss of Section 2 finally has all its parts.
There is also a pedagogical contrast waiting at the end. The problem of Section 4.1 was first order and spent its first fifteen hundred iterations on a plateau near $0.35$; this problem is fourth order and its loss is down to $10^{-4}$ within one thousand epochs. Order of differentiation, it turns out, is not what makes a PINN hard.
This is a one-dimensional problem: the beam occupies the interval $x\in[0,1]$ and the unknown is its transverse deflection $y(x)$.
Problem properties (nondimensionalized)
$$E = 1\,\text{Pa}, \quad I = 1\,\text{m}^4, \quad q = 1\,\text{N/m}, \quad l = 1\,\text{m}$$
Here $E$ is Young's modulus, $I$ the second moment of area of the cross-section, $q$ the intensity of a uniformly distributed load, and $l$ the length of the beam. All four have been set to unity. This is not laziness; it is the nondimensionalization urged in Section 2. With every quantity of order one, the equation residual and the boundary residuals have comparable magnitudes, and the equal weighting used in the loss below becomes defensible rather than arbitrary. If you repeat this laboratory with a steel beam in SI units, where $E \approx 2\times 10^{11}$, the loss terms will differ by twenty orders of magnitude and the training will fail for reasons that have nothing to do with the method.
Governing equation and boundary conditions
$${\partial^4 y \over \partial x^4} + q = 0, \qquad \text{where} \quad x \in [0,1]$$
This is the Euler-Bernoulli beam equation. In its dimensional form it reads $EI\, d^4y/dx^4 = -q$, and it is derived by combining three ingredients: the kinematic assumption that plane sections remain plane and normal to the deformed axis, so that curvature is $d^2y/dx^2$; the constitutive relation $M = EI\, d^2y/dx^2$ between bending moment and curvature; and the equilibrium relations $dM/dx = V$ and $dV/dx = -q$ between moment, shear force, and distributed load. Differentiating the constitutive relation twice and substituting the equilibrium relations gives the fourth-order equation. With $EI = 1$ and $q = 1$ it reduces to $y'''' = -1$.
The chain of derivatives is worth keeping in view, because it is what gives the boundary conditions their meaning:
$$ y \; \longrightarrow \; y' = \text{slope} \; \longrightarrow \; EIy'' = M \; \longrightarrow \; EIy''' = V. $$
Why four boundary conditions
A fourth-order ordinary differential equation must be integrated four times to recover $y$, and each integration introduces an arbitrary constant. Four constants require four independent conditions to pin them down. This is the general rule: the number of side conditions equals the order of the equation. In this problem they are supplied two at each end.
- One Dirichlet boundary condition on the left boundary:
$$y(0) = 0$$
- One Neumann boundary condition on the left boundary:
$$y'(0) = 0$$
- Two boundary conditions on the right boundary:
$$y''(1) = 0, \quad y'''(1) = 0$$
Read physically, the left end at $x=0$ is clamped, built into a wall: the wall permits neither displacement, hence $y(0)=0$, nor rotation, hence $y'(0)=0$. The right end at $x=1$ is free, hanging in the air: nothing there applies a bending moment, hence $EIy''(1) = M(1) = 0$, and nothing applies a transverse force, hence $EIy'''(1) = V(1) = 0$. The beam is therefore a cantilever carrying a uniformly distributed load, which is the first structure in every strength-of-materials course.
Notice how the two ends divide the four conditions. A clamped end constrains the kinematic quantities, displacement and slope; a free end constrains the static quantities, moment and shear. Every combination of end conditions in beam theory, simply supported, guided, pinned, is another way of choosing two from these four, and the choice is always made so that each end contributes exactly two.
The exact solution
$$y(x) = -{1 \over 24}x^4 + {1 \over 6}x^3 - {1 \over 4}x^2$$
Verify it rather than accepting it, since every later claim is measured against it. Differentiating successively,
$$ \begin{align*} y'(x) &= -\tfrac{1}{6}x^3 + \tfrac{1}{2}x^2 - \tfrac{1}{2}x,\\[2pt] y''(x) &= -\tfrac{1}{2}x^2 + x - \tfrac{1}{2} = -\tfrac{1}{2}(x-1)^2,\\[2pt] y'''(x) &= -x + 1 = -(x-1),\\[2pt] y''''(x) &= -1. \end{align*} $$
The last line gives $y'''' + q = -1 + 1 = 0$, so the equation holds. At $x=0$ we get $y(0)=0$ and $y'(0)=0$, the clamped conditions. At $x=1$ the factored forms make the free conditions immediate: $y''(1) = -\tfrac12 (0)^2 = 0$ and $y'''(1) = -(0) = 0$. All four conditions are satisfied, so this quartic is the solution, uniqueness being guaranteed by the well-posedness of the problem.
Two by-products of the algebra are worth recording, because they connect the exercise to results the reader already knows. The tip deflection is
$$ y(1) = -\frac{1}{24} + \frac{1}{6} - \frac{1}{4} = \frac{-1 + 4 - 6}{24} = -\frac{3}{24} = -\frac{1}{8}, $$
which is the classical cantilever result $\delta = ql^4/8EI$ with $q = l = EI = 1$, the negative sign recording that the deflection is downward. And at the clamped root, $M(0) = EIy''(0) = -\tfrac12$ and $V(0) = EIy'''(0) = 1$, which are the familiar fixed-end values $ql^2/2$ and $ql$. The network we are about to train has no access to any of this; it will be given only the differential equation and the four boundary conditions.
Since the deflection never leaves the interval $[-0.125,\, 0]$, keep that range in mind when the network architecture appears.
Make a neural network and loss functions like below:
Nothing new is imported; the single differentiation of Section 4.1 is wrapped in a reusable helper and then applied four times in succession.
Define Collocation Points
Twenty-five points, evenly spaced on $[0,1]$ with spacing $1/24 \approx 0.0417$, serve as the collocation set at which the equation residual will be enforced. The two ends, $x = 0$ and $x = 1$, are the boundary points.
Three remarks. First, the points are laid out on a regular lattice rather than sampled at random, and unlike Section 4.1 they are computed once and reused for every epoch. This makes the objective a fixed function of the weights and the gradient exact rather than noisy, and it is the first of the differences that will explain this laboratory's much better convergence. Second, twenty-five points for a fourth-order equation is a small budget; it works here because the solution is a quartic polynomial, about as gentle a function as exists. Third, the two boundary points remain members of the collocation set, so the equation residual is enforced at $x=0$ and $x=1$ as well as the boundary conditions. Nothing is wrong with that: the equation does hold at the endpoints.
The plot marks the interior points in blue and the two boundaries in red, which is a habit worth keeping in higher dimensions, where an error in constructing the boundary sets is otherwise invisible until training fails.
x_pts = np.linspace(0, 1, 25)
plt.figure(figsize = (6, 3))
plt.scatter(x_pts, np.zeros_like(x_pts), s = 10)
plt.scatter(0, 0, label = 'x = 0', s = 10, c = 'red')
plt.scatter(1, 0, label = 'x = 1', s = 10, c = 'red')
plt.xlabel('x')
plt.legend()
plt.title('Collocation Points with Boundaries')
plt.show()
Define Network and Hyper-parameters
The network maps one input to one output through four hidden layers of thirty-two units with $\tanh$ activations, a little larger than the network of Section 4 and still very small by the standards of the rest of this course.
The architecture contains one feature that should stop you: there is a $\tanh$ after the final linear layer. The output of this network is therefore confined to the interval $(-1,1)$.
In Section 4.1 the same choice would have been fatal, because the solution there took values near $1.16$. Here it is survivable, because the true deflection lies in $[-0.125, 0]$, comfortably inside the reachable range. It is nevertheless a fragile choice, and worth understanding as a hazard rather than a technique. If the load $q$ were increased by a factor of ten, the true tip deflection would be $-1.25$ and no setting of the weights could produce it; the training would converge to something pressed against the saturation of the output $\tanh$, and the failure would look like a modelling error rather than an architectural one. A second, subtler cost is that $\tanh$ has small derivative when its argument is large, so gradients reaching the network through a saturated output layer are attenuated.
The optimizer is Adam at $10^{-3}$, and every term of the composite loss will be a mean squared residual.
torch.manual_seed(0)
model = nn.Sequential(
nn.Linear(1, 32),
nn.Tanh(),
nn.Linear(32, 32),
nn.Tanh(),
nn.Linear(32, 32),
nn.Tanh(),
nn.Linear(32, 32),
nn.Tanh(),
nn.Linear(32, 1),
nn.Tanh()
)
optimizer = torch.optim.Adam(model.parameters(), lr = 1e-3)
Define a Function to Compute Derivatives
One short helper does all the work of this and every remaining laboratory. It returns the derivative of one tensor with respect to another, and it is the operation of Section 4.1 given a name. Two details of it carry the method.
The differentiation is asked for at every collocation point at once, and the vector of ones passed to it extracts exactly the per-point derivatives. The reason this works is worth understanding once, properly: the underlying routine computes a vector-Jacobian product rather than a Jacobian, returning $g^{\top}\,\partial y/\partial t$ for a supplied vector $g$, and choosing $g = \mathbf{1}$ returns the column sums of the Jacobian. That looks like the wrong object until one notices that the network is applied pointwise to a batch, so the $i$-th output depends only on the $i$-th input, the Jacobian is diagonal, and each column sum collapses to its single surviving diagonal entry $\partial y_i/\partial t_i$. This same trick recurs in every PINN implementation in every framework.
The differentiation is also asked to build a differentiable graph for itself. Without that, the returned derivative would carry no history, and asking for its derivative in turn, which is precisely what the fourth-order chain does next, would fail. It is also what allows the final loss, built from these derivatives, to be differentiated with respect to the network weights. Every higher-order derivative in this notebook depends on it.
Any tensor with respect to which we intend to differentiate must be created with requires_grad = True, for the reason given in Section 4.1; the training cell below does this for the three point sets.
def derivative(y, t):
return torch.autograd.grad(y, t, grad_outputs = torch.ones_like(y), create_graph = True)[0]
Define PDE System
$$\frac{\partial^4 y}{\partial x^4} + q = 0, \qquad \text{where} \quad x \in [0,1]$$
The function below evaluates the residual by applying the derivative helper four times in succession, each call differentiating the output of the previous one. The returned quantity is $y'''' + 1$, that is, $y'''' + q$ with $q = 1$, and it should be zero when the equation holds.
This is where the cost of a high-order equation becomes visible. Each call extends the computational graph, so the graph differentiated at the end of the step is several times as deep as the network itself, and the memory and time per epoch scale accordingly. It is also where the smoothness requirement of Section 1.2 becomes non-negotiable: a ReLU network differentiated four times returns zero everywhere, so the residual would be identically $1$ and no amount of training would reduce it.
Note finally that the residual is returned unsquared and unaveraged. Squaring and averaging happen in the training loop.
def PDE(model, x):
u = model(x)
u_x = derivative(u, x)
u_xx = derivative(u_x, x)
u_xxx = derivative(u_xx, x)
u_xxxx = derivative(u_xxx, x)
pde = u_xxxx + 1
return pde
Define Boundary Condition
- One Dirichlet boundary condition on the left boundary:
$$y(0) = 0$$
- One Neumann boundary condition on the left boundary:
$$y'(0) = 0$$
- Two boundary conditions on the right boundary:
$$y''(1) = 0, \quad y'''(1) = 0$$
The two functions below return the four residuals: at the left boundary, the value of the network and its first derivative; at the right boundary, its second and third derivatives. Because all four prescribed values are zero, each residual is simply the corresponding derivative of the network. Had the beam been clamped at a prescribed non-zero settlement $y(0) = \delta$, the residual would have been $y_\theta(0) - \delta$, following the pattern of Section 2 exactly; the homogeneous case merely makes the subtraction invisible.
The homogeneity has a consequence for the training output that is worth predicting now. An untrained network produces values near zero, so all four boundary residuals begin small, and essentially all of the initial loss will come from the equation term.
def BC_left(model, x_0):
# x = 0
u_0 = model(x_0)
bc1 = u_0
u_x_0 = derivative(u_0, x_0)
bc2 = u_x_0
return bc1, bc2
def BC_right(model, x_f):
# x = 1
u_f = model(x_f)
u_x_f = derivative(u_f, x_f)
u_xx_f = derivative(u_x_f, x_f)
bc3 = u_xx_f
u_xxx_f = derivative(u_xx_f, x_f)
bc4 = u_xxx_f
return bc3, bc4
Train
The three point sets are created as tracked tensors and the loop runs for a thousand and one epochs. Each epoch computes the equation residual on the twenty-five collocation points and the four boundary residuals, takes the mean square of each, and adds them:
$$ \mathcal{L} = \mathcal{L}_{\text{PDE}} + \mathcal{L}_{\text{BC1}} + \mathcal{L}_{\text{BC2}} + \mathcal{L}_{\text{BC3}} + \mathcal{L}_{\text{BC4}}. $$
All five weights are $1$, which is the nondimensionalization argument from Section 5.1 being cashed in.
Predict the first printed line
At epoch $0$ the network output and all its derivatives are near zero, so $y'''' \approx 0$ and the residual is $0 + 1 = 1$, giving $\mathcal{L}_{\text{PDE}} \approx 1$. The four boundary residuals are near zero, so $\mathcal{L}_{\text{BC}} \approx 0$ and the total is about $1$.
x = torch.tensor(x_pts.reshape(-1, 1), dtype = torch.float32, requires_grad = True)
x0 = torch.tensor([[0.0]], requires_grad = True)
x1 = torch.tensor([[1.0]], requires_grad = True)
for epoch in range(1001):
pde = PDE(model, x)
bc1, bc2 = BC_left(model, x0)
bc3, bc4 = BC_right(model, x1)
loss_pde = torch.mean(pde**2)
loss_bc = torch.mean(bc1**2) + torch.mean(bc2**2) + torch.mean(bc3**2) + torch.mean(bc4**2)
loss = loss_pde + loss_bc
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 500 == 0:
print('Epoch: {} Loss: {:.4f} PDELoss: {:.4f} BCLoss: {:.4f}'.format(epoch,
loss.item(),
loss_pde.item(),
loss_bc.item()))
The cell prints a total of $1.0029$ at epoch $0$, of which $0.9173$ is the equation term and $0.0855$ the boundary terms, which is the predicted split: the equation term near $1$ and dominant, the boundary terms small. (The equation term is somewhat below $1$ because the freshly initialised network's fourth derivative is not exactly zero.) By epoch $500$ the total is $0.0005$, and by epoch $1000$ it is $0.0001$, with both the equation and the boundary terms printing as $0.0000$ to four decimal places: the network satisfies the fourth-order equation and all four boundary conditions to about that accuracy.
Compare this with Section 4.1, where a first-order equation, after three times as many iterations, was still reporting a loss around $0.35$ at its thousandth iteration, and did not leave that level until about iteration $1500$. The order of differentiation was clearly not the obstacle there. Three differences account for the contrast. The collocation points here are fixed rather than resampled, so the gradient is exact and Adam is minimizing a genuine function rather than chasing a moving target, and the loss printed here is correspondingly an exact value rather than a noisy estimate. The solution here is a quartic polynomial, smooth and monotone over the whole domain, whereas the solution in Section 4.1 oscillated through two full periods, and a network must work considerably harder to produce an oscillation than a bend. And this network is deeper, with four hidden layers against two. The general lesson, which PINN: Limitations and Extensions develops, is that PINN difficulty tracks the character of the solution and the conditioning of the optimization far more than it tracks the order of the equation.
Define Exact Solution
The quartic verified in Section 5.1 is coded up for comparison. It is used only for plotting, after training is complete, and never enters the loss; the network was trained on the equation and the boundary conditions alone.
def exact_solution(x):
return -(x ** 4) / 24 + x ** 3 / 6 - x ** 2 / 4
Result
The prediction is evaluated on a hundred points, four times as many as the network was trained on and at positions that mostly fall between the training points. This is the meaningful test. A model that had merely driven its residual to zero at twenty-five isolated locations, by whatever contortion, would betray itself on the finer grid; a model that has actually learned the solution will interpolate smoothly, because the true solution is smooth and so is the network.
x_ = np.linspace(0, 1, 100)
x_ = torch.tensor(x_.reshape(-1, 1), dtype = torch.float32)
with torch.no_grad():
u_ = model(x_)
plt.figure(figsize = (6, 4))
plt.plot(x_, exact_solution(x_), 'k', label = 'Exact Solution', linewidth = 5, alpha = 0.3)
plt.plot(x_, u_.numpy(), 'r--', label = 'Predicted u')
plt.legend(fontsize = 10)
plt.xticks(fontsize = 10)
plt.yticks(fontsize = 10)
plt.xlabel('x', fontsize = 10)
plt.ylabel('Displacement, y', fontsize = 10)
plt.show()
The exact solution and the prediction are drawn on the same axes. Given a training loss of about $10^{-4}$, the two curves should be visually indistinguishable over the whole span, falling from zero at the clamped root to about $-0.125$ at the free tip.
It is worth appreciating what has happened, because it is easy to lose sight of it among the code. A small network was shown no solution, no data, and no finite element mesh. It was shown a differential equation, four boundary conditions, and twenty-five bare coordinates. From those it recovered the deflection curve of a loaded cantilever, including the correct tip deflection $ql^4/8EI$, which is a result that took the nineteenth century some effort to establish. The supervision came entirely from the physics.
Overview
Every example so far has been one-dimensional: a single independent variable, a single unknown field. The thin plate is the first two-dimensional problem, and the step up costs more than an extra input. There are now two coupled unknown fields, the displacements $u$ and $v$; there are mixed partial derivatives; the boundary is a curve rather than two points, so each edge needs its own point set; and the boundary conditions are expressed in terms of stresses, which are themselves combinations of derivatives of the unknowns. Every one of these complications is handled by the machinery already built, which is the point of doing the problem.
Section 6.2 solves the plate as a pure equation solver, with no data of any kind, and compares against a finite element solution. Section 6.3 then adds sensor data to the same model and measures the improvement. That comparison is the most valuable part of this section, and it is the cleanest demonstration in the notebook of what the data term actually contributes.
We will solve the thin plate equations to find the displacement and stress distribution of a thin plate in two dimensions.
The formulation rests on Kirchhoff-Love plate theory, which makes three kinematic hypotheses:
- straight lines normal to the mid-surface remain straight after deformation;
- straight lines normal to the mid-surface remain normal to the mid-surface after deformation;
- the thickness of the plate does not change during a deformation.
These three assumptions are the plate analogue of the assumption that plane sections remain plane in beam theory, and they do the same work. They say that the displacement everywhere through the thickness is determined by the displacement of the mid-surface, so a three-dimensional elasticity problem collapses into a two-dimensional one posed on that surface. Without them the unknowns would be functions of three variables and the network would need a third input. They are accurate when the plate is thin, which is the case here: the thickness is $h = 1$ mm against a width of $20$ mm.
A non-uniform stretching force is applied to a square elastic plate. Because both the geometry and the in-plane loading are symmetric about the two centre-lines, only one quarter of the plate need be considered, the yellow region in the figure below. This is standard practice in finite element analysis, and it is worth being precise about why it is legitimate and what replaces the material that has been removed. Symmetry of geometry and loading implies symmetry of the solution, and a symmetric displacement field has no component normal to a plane of symmetry: material on the symmetry plane cannot move across it, since the mirror-image material on the other side would have to move the opposite way. The removed three quarters are therefore replaced by exactly two conditions, $u = 0$ on the plane $x = 0$ and $v = 0$ on the plane $y = 0$. Those are the first two boundary conditions listed below, and reading them as symmetry conditions rather than as clamps is essential to understanding the problem.
Problem properties
$$ E = 50\,\text{MPa}, \quad \nu = 0.3, \quad \omega = 20\,\text{mm}, \quad h = 1\,\text{mm}, \quad f = 1\,\text{MPa} $$
Here $E$ is Young's modulus, $\nu$ Poisson's ratio, $\omega$ the width of the full plate, $h$ its thickness, and $f$ the magnitude of the applied in-plane traction. Since only a quarter is modelled, the computational domain is the square $[0, \omega/2] \times [0, \omega/2] = [0,10]\times[0,10]$ in millimetres. Working in millimetres and megapascals keeps every number within an order of magnitude or two of unity, which is again the nondimensionalization argument of Section 2 applied by a judicious choice of units rather than by explicit scaling.
Governing equations
The governing equations (Föppl-von Kármán equations) for the isotropic elastic plate are
$$ \begin{align*} &{E \over 1 - \nu^2}\left({\partial^2 u \over \partial x^2} + {1 - \nu \over 2}{\partial^2 u \over \partial y^2} + {1 + \nu \over 2}{\partial^2 v \over \partial x \partial y} \right) = 0\\\\ &{E \over 1 - \nu^2}\left({\partial^2 v \over \partial y^2} + {1 - \nu \over 2}{\partial^2 v \over \partial x^2} + {1 + \nu \over 2}{\partial^2 u \over \partial x \partial y} \right) = 0 \end{align*} $$
These are the equations of in-plane equilibrium written in terms of the displacement components, and they are worth deriving rather than accepting, both because the derivation is short and because it explains where the awkward coefficients $\tfrac{1-\nu}{2}$ and $\tfrac{1+\nu}{2}$ come from.
Start from equilibrium in stress form. With no body force, a material element balances if
$$ \frac{\partial \sigma_{xx}}{\partial x} + \frac{\partial \sigma_{xy}}{\partial y} = 0, \qquad \frac{\partial \sigma_{xy}}{\partial x} + \frac{\partial \sigma_{yy}}{\partial y} = 0. $$
Next substitute the plane-stress constitutive relations, which express stress in terms of the displacement gradients:
$$ \begin{align*} \sigma_{xx} &= \frac{E}{1 - \nu^2} \left( \frac{\partial u}{\partial x} + \nu \frac{\partial v}{\partial y} \right), \\ \sigma_{yy} &= \frac{E}{1 - \nu^2} \left( \frac{\partial v}{\partial y} + \nu \frac{\partial u}{\partial x} \right), \\ \sigma_{xy} &= \frac{E}{2(1 + \nu)} \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right). \end{align*} $$
Differentiating the first with respect to $x$ and the third with respect to $y$ and adding gives
$$ \frac{E}{1-\nu^2}\left(\frac{\partial^2 u}{\partial x^2} + \nu \frac{\partial^2 v}{\partial x \partial y}\right) + \frac{E}{2(1+\nu)}\left(\frac{\partial^2 v}{\partial x \partial y} + \frac{\partial^2 u}{\partial y^2}\right) = 0. $$
To combine the two brackets, note that the shear modulus prefactor can be rewritten with the same denominator as the first, since $1-\nu^2 = (1-\nu)(1+\nu)$ and therefore
$$ \frac{E}{2(1+\nu)} = \frac{E}{1-\nu^2}\cdot\frac{1-\nu}{2}. $$
Factoring $E/(1-\nu^2)$ out of both terms leaves
$$ \frac{E}{1-\nu^2}\left(\frac{\partial^2 u}{\partial x^2} + \frac{1-\nu}{2}\frac{\partial^2 u}{\partial y^2} + \left[\nu + \frac{1-\nu}{2}\right]\frac{\partial^2 v}{\partial x \partial y}\right) = 0, $$
and the coefficient of the mixed derivative simplifies:
$$ \nu + \frac{1-\nu}{2} = \frac{2\nu + 1 - \nu}{2} = \frac{1+\nu}{2}. $$
This is exactly the first governing equation. The second follows in the same way from the second equilibrium equation, with the roles of $x$ and $y$, and of $u$ and $v$, interchanged. The mixed derivative $\partial^2 u / \partial x \partial y$ in that equation is what couples the two: neither displacement component can be solved for independently, which is why the network has two outputs rather than two networks with one each.
Boundary conditions
- Two Dirichlet boundary conditions at $x = 0$ and $y = 0$ (B.C. ① and B.C. ②):
$$ \begin{align*} v(x,y) &= 0 \qquad \text{at} \quad y = 0\\\\ u(x,y) &= 0 \qquad \text{at} \quad x = 0 \end{align*} $$
As explained above, these are the symmetry conditions that stand in for the three quarters of the plate we discarded, not physical restraints.
- Two free boundary conditions at $y = \omega / 2$ (B.C. ③):
$$ \sigma_{yy} = 0,\quad \sigma_{yx} = 0 $$
The top edge is a free surface: nothing pushes on it and nothing drags along it, so both the normal and the shear traction on a plane with normal $y$ vanish.
- Free boundary condition and in-plane force boundary condition at $x = \omega / 2$ (B.C. ④):
$$\sigma_{xx} = P \cdot h,\quad \sigma_{xy} = 0$$
The right edge is where the plate is pulled. The normal stress there is prescribed by the applied load and the shear stress is zero. The load is the non-uniform stretching force mentioned above, and the code makes its distribution explicit: the applied traction is $f h \cos\!\left(\pi y / 2b\right)$ with $b = \omega/2 = 10$, taking its maximum value $fh = 1$ at $y = 0$, the centre-line of the full plate, and falling to zero at $y = 10$, where it must vanish to be compatible with the free top edge. A discontinuity between a non-zero pull on one edge and a traction-free adjacent edge would produce a stress singularity at the corner, and choosing a cosine distribution avoids it.
Note that these last four conditions constrain stresses, not displacements, and stresses are built from first derivatives of $u$ and $v$. In the language of Section 2 they are Neumann-type conditions, and implementing them requires nothing new: the network is differentiated at the boundary points, the derivatives are assembled into the stress components using the constitutive relations above, and the result is compared against the prescribed value.
The network and the loss terms are defined below.
Numerical Solution (treated as the exact solution)
A finite element solution of the same problem is shown below, giving the displacement components $u$ and $v$ and the stresses $\sigma_{xx}$ and $\sigma_{yy}$ over the quarter domain. We will solve the problem with a PINN and compare against it.
Two habits of mind are worth adopting here. First, the finite element result is called the exact solution throughout this section, and it is not; it is a numerical approximation with its own discretization error, and it is being used as a reference because it is the best available answer, not because it is the truth. Second, look at the four fields before reading on and predict which will be hardest for the network. The plate is being pulled in $x$, so $u$ is the large, smooth, directly driven response, while $v$ is the smaller transverse contraction produced by the Poisson effect, an order of magnitude smaller and therefore contributing correspondingly less to a loss that weights all residuals equally. That prediction is confirmed in Section 6.3, where it turns out to be exactly the field the data term is needed for.
Select Device & Define Properties
The first cell selects a device, since this laboratory is the one that benefits from a GPU; the second defines the material and geometric constants, which the residual and boundary functions will read. The constants are $E = 50$, $\nu = 0.3$, $a = b = 10$, $h = 1$ and $f = 1$, matching Section 6.1 in millimetre-megapascal units. Only $b$ is used in what follows, since the quarter domain is square.
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
print(device)
# properties
E = 50.0
nu = 0.3
a = 10.0
b = 10.0
h = 1.0
f = 1.0
Define Collocation Points
A $50 \times 50$ lattice on $[0,10]^2$ gives $2500$ interior collocation points, and the four edges are then extracted by testing each coordinate against $0$ or $b$, giving fifty points per edge.
Note that the four corner points satisfy two edge tests each, so they appear in two of the boundary sets and are constrained twice. This is harmless, and it is in fact what one wants, since both conditions do hold at a corner.
The scatter plot colours the interior grey and each edge differently. In two dimensions this check is no longer optional. An error in a boundary mask does not raise an exception; it silently produces a model that solves a different problem, and the plot is the cheapest way to catch it.
Nx = 50 # Number of samples
Ny = 50 # Number of samples
x = torch.linspace(0, b, Nx) # Input data for x (Nx x 1)
y = torch.linspace(0, b, Ny) # Input data for y (Ny x 1)
xy = torch.meshgrid(x, y, indexing = 'ij')
Y_domain = torch.cat([xy[0].reshape(-1, 1), xy[1].reshape(-1, 1)], dim = 1)
Y_left = Y_domain[Y_domain[:, 0] == 0]
Y_right = Y_domain[Y_domain[:, 0] == b]
Y_top = Y_domain[Y_domain[:, 1] == b]
Y_bottom = Y_domain[Y_domain[:, 1] == 0]
plt.figure(figsize = (5, 5))
plt.scatter(Y_domain[:, 0], Y_domain[:, 1], s = 4, c = 'gray', alpha = 0.3)
plt.scatter(Y_bottom[:, 0], Y_bottom[:, 1], s = 4, c = 'red')
plt.scatter(Y_top[:, 0], Y_top[:, 1], s = 4, c = 'blue')
plt.scatter(Y_left[:, 0], Y_left[:, 1], s = 4, c = 'green')
plt.scatter(Y_right[:, 0], Y_right[:, 1], s = 4, c = 'purple')
plt.title('Collocation Points')
plt.xlabel('X')
plt.ylabel('Y')
plt.grid(True)
plt.show()
Define Network and Hyper-parameters
The network takes the two coordinates as input and returns the two displacement components as output, through six hidden layers of sixty-four $\tanh$ units. Three points of construction deserve comment.
The two outputs are produced by one network with a shared trunk, not by two independent networks. This matters because the governing equations are coupled through the mixed derivative $\partial^2 v/\partial x\partial y$, so the two fields must be learned together in any case, and a shared representation is the natural choice.
The weights are initialized by the Xavier, or Glorot, scheme, which scales the weight variance by the fan-in and fan-out of each layer so that the variance of the activations neither grows nor decays with depth. With six hidden layers this matters more than it did with two. The network is written as a class, because the pre-trained models loaded later in this section were saved as whole objects of this class and are unpickled against it.
There is no activation after the final linear layer, so the outputs are unbounded. Contrast Section 5, where a final $\tanh$ was tolerable only because the answer happened to be small. Here the displacements are of order $0.1$ and the unbounded choice is simply correct.
class PINN(nn.Module):
def __init__(self):
super(PINN, self).__init__()
self.net = nn.Sequential(
nn.Linear(2, 64),
nn.Tanh(),
nn.Linear(64, 64),
nn.Tanh(),
nn.Linear(64, 64),
nn.Tanh(),
nn.Linear(64, 64),
nn.Tanh(),
nn.Linear(64, 64),
nn.Tanh(),
nn.Linear(64, 64),
nn.Tanh(),
nn.Linear(64, 2)
)
for m in self.net:
if isinstance(m, nn.Linear):
nn.init.xavier_normal_(m.weight)
nn.init.zeros_(m.bias)
def forward(self, x):
output = self.net(x)
u, v = output[:, 0:1], output[:, 1:2]
return u, v
model = PINN().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr = 1e-3)
Define Functions to Compute Physics Loss
The derivative helper is the two-dimensional version of the one from Section 5. The gradient with respect to the full input has two columns, so the helper slices out the two partial derivatives and returns them as a pair, $\partial/\partial x$ first and $\partial/\partial y$ second. Everything said in Section 5 about the vector of ones and about recording the differentiation applies unchanged; the tensor of ones is created on the same device as the model.
Returning both partials at once has a pleasant consequence for second derivatives. Applying the helper to $\partial u/\partial x$ returns $\partial^2 u/\partial x^2$ and $\partial^2 u/\partial x \partial y$ in one call, so the mixed derivative is obtained for free alongside the pure one, and the residual function below exploits this.
def derivative(y, t):
df = torch.autograd.grad(y, t, grad_outputs = torch.ones_like(y), create_graph = True)[0]
df_x = df[:, 0:1]
df_y = df[:, 1:2]
return df_x, df_y
Define PDE with Boundary Conditions
Governing equations (Föppl-von Kármán equations) for the isotropic elastic plate:
$$ \begin{align*} &{E \over 1 - \nu^2}\left({\partial^2 u \over \partial x^2} + {1 - \nu \over 2}{\partial^2 u \over \partial y^2} + {1 + \nu \over 2}{\partial^2 v \over \partial x \partial y} \right) = 0\\\\ &{E \over 1 - \nu^2}\left({\partial^2 v \over \partial y^2} + {1 - \nu \over 2}{\partial^2 v \over \partial x^2} + {1 + \nu \over 2}{\partial^2 u \over \partial x \partial y} \right) = 0 \end{align*} $$
The function below computes six second derivatives and assembles the two residuals. Read it with the derivation of Section 6.1 alongside and the correspondence is exact. One naming convention will otherwise confuse you: the $x$-displacement is called $w$ throughout the code, so a second derivative named dw_xx is $\partial^2 u/\partial x^2$ and not a derivative of any out-of-plane deflection. The letter $w$ conventionally denotes transverse deflection in plate theory, and it is used here for the in-plane $u$ instead.
The residuals are returned unsquared, as before.
def PDE(model, Y):
w, v = model(Y)
dw_x, dw_y = derivative(w, Y)
dw_xx, dw_xy = derivative(dw_x, Y)
_, dw_yy = derivative(dw_y, Y)
dv_x, dv_y = derivative(v, Y)
dv_xx, dv_xy = derivative(dv_x, Y)
_, dv_yy = derivative(dv_y, Y)
force_eq_x = (dw_xx + 0.5*(1 - nu)*dw_yy + 0.5*(1 + nu)*dv_xy)*E/(1 - nu**2)
force_eq_y = (dv_yy + 0.5*(1 - nu)*dv_xx + 0.5*(1 + nu)*dw_xy)*E/(1 - nu**2)
return force_eq_x, force_eq_y
Two Dirichlet boundary conditions at $x = 0$ and $y = 0$ (B.C. ①, left, and B.C. ②, bottom):
$$ \begin{align*} v(x,y) &= 0 \qquad \text{at} \quad y = 0 \text{: (bottom)} \\ \\ u(x,y) &= 0 \qquad \text{at} \quad x = 0 \text{: (left)} \end{align*} $$
These are the simplest residuals in the notebook. Each evaluates the network on the relevant edge and returns the component that should vanish, discarding the other. Note that only one of the two displacement components is constrained on each edge: the left edge forbids motion across the symmetry plane but permits sliding along it, and a condition on $v$ there would be wrong, since material on the symmetry plane certainly does move in $y$ when the plate is stretched.
def BC_left(model, Y_left):
# x = 0
w_left, _ = model(Y_left)
return w_left
def BC_bottom(model, Y_bottom):
# y = 0
_, v_bottom = model(Y_bottom)
return v_bottom
Two free boundary conditions at $y = \omega / 2$ (B.C. ③, top):
$$ \sigma_{yy} = 0,\quad \sigma_{yx} = 0 $$
Free boundary condition and in-plane force boundary condition at $x = \omega / 2$ (B.C. ④, right):
$$ \begin{align*} \sigma_{xx} &= P \cdot h \\ \sigma_{xy} &= 0 \end{align*} $$
Stress
$$ \begin{align*} \sigma_{xx} &= \frac{E}{1 - \nu^2} \left( \frac{\partial u}{\partial x} + \nu \frac{\partial v}{\partial y} \right) \\ \sigma_{yy} &= \frac{E}{1 - \nu^2} \left( \frac{\partial v}{\partial y} + \nu \frac{\partial u}{\partial x} \right)\\ \sigma_{xy} &= \sigma_{yx} = \frac{E}{2(1 + \nu)} \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right) \end{align*} $$
The next cell implements the two edges. On each, the network is evaluated, its first derivatives are taken, the stress components are assembled from them, and the residual is the stress minus its prescribed value. On the top edge both prescribed values are zero, so the residuals are the stresses themselves. On the right edge the shear residual is the stress itself, while the normal residual is the difference between $\sigma_{xx}$ and
$$ \sigma_{\text{ex}}(y) = f\,h\,\cos\!\left(\frac{\pi y}{2b}\right), $$
which is the non-uniform applied traction described in Section 6.1, equal to $1$ at $y=0$ and to $0$ at $y=b=10$.
This cell is the clearest illustration in the notebook of a point worth stating explicitly. A conventional numerical method must treat a Neumann condition specially, because its discrete unknowns are nodal values and a derivative at the boundary has to be reconstructed from them by a one-sided formula. A PINN needs nothing of the kind. A stress boundary condition is written by differentiating the network at the boundary points and combining the derivatives with the constitutive law, which is the same operation used in the interior, applied to a different point set. In the PINN framework, a condition on a derivative is no harder than a condition on a value.
def BC_top(model, Y_top):
# y = b
w, v = model(Y_top)
dw_x, dw_y = derivative(w, Y_top)
dv_x, dv_y = derivative(v, Y_top)
sig_xy = (dv_x + dw_y)*E/(1 + nu)/2
sig_yy = (dv_y + nu*dw_x)*E/(1 - nu**2)
return sig_yy, sig_xy
def BC_right(model, Y_right):
# x = b
w, v = model(Y_right)
dw_x, dw_y = derivative(w, Y_right)
dv_x, dv_y = derivative(v, Y_right)
sig_xx = (dw_x + nu*dv_y)*E/(1 - nu**2)
sig_xy = (dv_x + dw_y)*E/(1 + nu)/2
sig_ex = f*h*torch.cos(np.pi/(2*b)*Y_right[:, 1:2])
return sig_xx - sig_ex, sig_xy
Train
The training loop assembles eight residual terms: two equation residuals over the $2500$ interior points, one on each of the two symmetry edges, two on the top edge, and two on the right edge. Each becomes a mean squared error against zero and all are added with unit weight, giving
$$ \mathcal{L} = \underbrace{\mathcal{L}_{\text{PDE},x} + \mathcal{L}_{\text{PDE},y}}_{\text{equilibrium}} \; + \; \underbrace{\mathcal{L}_{\text{bottom}} + \mathcal{L}_{\text{left}} + \mathcal{L}_{\text{top},1} + \mathcal{L}_{\text{top},2} + \mathcal{L}_{\text{right},1} + \mathcal{L}_{\text{right},2}}_{\text{boundary}}. $$
Training a PINN can be computationally expensive and time-consuming. To facilitate experimentation and reproducibility, a trained model is provided below, so that the results can be evaluated and analysed without waiting for a long run. You may still run the loop yourself to watch the behaviour; the printed trace is from a run on a GPU.
Before running anything, notice one consequence of the unit weighting that the discussion of Section 6.1 anticipated. The equation residuals are multiplied by $E/(1-\nu^2) \approx 54.9$, so a given displacement error contributes far more through the equilibrium terms than through the displacement boundary terms. The stress boundary conditions carry the same large factor. The transverse displacement $v$, being the small Poisson response, influences all of these only weakly, and it is therefore the field the optimizer has the least incentive to get right. Keep this in mind when the results appear.
Y_domain = Y_domain.to(device).requires_grad_(True)
Y_bottom = Y_bottom.to(device).requires_grad_(True)
Y_top = Y_top.to(device).requires_grad_(True)
Y_left = Y_left.to(device).requires_grad_(True)
Y_right = Y_right.to(device).requires_grad_(True)
epochs = 10001
for epoch in range(epochs):
# PDE loss
force_eq_x, force_eq_y = PDE(model, Y_domain)
loss_PDE = torch.mean(force_eq_x**2) + torch.mean(force_eq_y**2)
# BC loss
bottom_v = BC_bottom(model, Y_bottom)
left_w = BC_left(model, Y_left)
top_sigyy, top_sigxy = BC_top(model, Y_top)
right_sigxx, right_sigxy = BC_right(model, Y_right)
loss_BC = (torch.mean(bottom_v**2) + torch.mean(left_w**2)
+ torch.mean(top_sigyy**2) + torch.mean(top_sigxy**2)
+ torch.mean(right_sigxx**2) + torch.mean(right_sigxy**2))
loss = loss_PDE + loss_BC
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 1000 == 0:
print('Epoch: {} Loss: {:.4f} PDELoss: {:.4f} BCLoss: {:.4f}'.\
format(epoch, loss.item(), loss_PDE.item(), loss_BC.item()))
Load the Model Parameters
The trained model will be used to avoid the long training time.
The first cell downloads the file; the second loads the model onto the selected device. The file holds the whole PINN object, so the class above must have been defined before it is loaded.
!gdown 1u64cvvA4GNSl3un0L8UAz40Pa4hSMSQx -O /content/lab2_adam_model.pt
model = torch.load('/content/lab2_adam_model.pt', map_location = device, weights_only = False)
Plot Displacement and Stress
Two helpers follow. The first evaluates the network at a set of points, takes the two first derivatives it needs, and returns $\sigma_{xx}$ and $\sigma_{yy}$ from the constitutive relations; this is the same computation performed inside the right-edge boundary function, now applied over the whole domain rather than on an edge. The second draws two scatter maps side by side, each with its own fixed colour range.
Observe that the stresses are derived quantities here. The network outputs displacements only; the stress field is obtained afterwards by differentiating it. A finite element code proceeds the same way, computing nodal displacements first and recovering stresses from them, with the difference that the finite element derivative is taken of a piecewise-polynomial interpolant and is discontinuous between elements, while the PINN derivative is taken of a smooth function and is exact.
The test grid is then refined to $100 \times 100$, giving ten thousand evaluation points, four times as many as were used for training. As in Section 5, most of them lie between training points, so this is a genuine test of whether the network has learned a field rather than memorized a lattice. Note that the test points carry requires_grad = True even though nothing is trained here: a stress is a derivative, and a derivative needs a graph.
The colour limits are fixed by hand so that the PINN maps can be compared directly against the finite element maps shown earlier: $u$ on $[0, 0.182]$, $v$ on $[-0.06, 0.011]$, $\sigma_{xx}$ on $[-0.0022, 1.0]$ and $\sigma_{yy}$ on $[-0.15, 0.45]$. Fixing the scale by hand is essential for a comparison of this kind, since a colour map with automatically chosen limits makes every field look equally well resolved.
def check_stress(model, Y):
w, v = model(Y)
w_x, _ = derivative(w, Y)
_, v_y = derivative(v, Y)
sig_xx = E/(1 - nu**2)*(w_x + nu*v_y)
sig_yy = E/(1 - nu**2)*(v_y + nu*w_x)
return sig_xx, sig_yy
def plot_results(points, data, color_legend, titles):
plt.figure(figsize = (8, 3))
for idx in range(2):
plt.subplot(1, 2, idx + 1)
plt.scatter(points[:, 0], points[:, 1],
c = data[idx].detach().cpu().numpy(),
cmap = 'rainbow',
s = 5)
plt.clim(color_legend[idx])
plt.title(titles[idx])
plt.axis('square')
plt.xlim((0, b))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.show()
Nx = 100
Ny = 100
x_test = torch.linspace(0, b, Nx)
y_test = torch.linspace(0, b, Ny)
Y_test = torch.meshgrid(x_test, y_test, indexing = 'ij')
test = torch.cat([Y_test[0].reshape(-1, 1), Y_test[1].reshape(-1, 1)], dim = 1)
test_tensor = test.clone().to(device).requires_grad_(True)
u_test, v_test = model(test_tensor)
sigma_xx, sigma_yy = check_stress(model, test_tensor)
titles = ['x-displacement ($u$)', 'y-displacement ($v$)',
'$\sigma_{xx}$', '$\sigma_{yy}$']
color_legend = [[0, 0.182], [-0.06, 0.011], [-0.0022, 1.0], [-0.15, 0.45]]
plot_results(test, [u_test, v_test], color_legend[:2], titles[:2])
plot_results(test, [sigma_xx, sigma_yy], color_legend[2:], titles[2:])
Read the four maps against the physics. The plate is pulled in $x$, so $u$ should grow from zero on the left symmetry edge to its largest value on the loaded right edge, and $\sigma_{xx}$ should be largest near $y=0$ where the applied traction is largest and fall towards zero at the top edge, following the imposed cosine. The transverse displacement $v$ is negative, the plate contracting laterally as it is stretched, and its magnitude is a few tenths of the magnitude of $u$, which is the Poisson effect with $\nu = 0.3$.
We now generate two sets of plots, one for the Finite Element Method (FEM) and another for the PINN, to compare their performance quantitatively. Additionally, to enable a more precise comparison, we plot the solution values along a diagonal line. A curve along a line is a far more sensitive instrument than two colour maps side by side: differences of a few percent that are invisible in a rainbow colormap are obvious as a gap between two lines.
Comparison with FEM data along a diagonal line
The finite element data file contains, for each node, the two coordinates, the two displacement components, and the two stress components, and the cell below splits it accordingly. The array named w again holds the $x$-displacement $u$, following the convention adopted in the network code.
The line used for the comparison is the anti-diagonal $x + y = 10$, running from the corner $(0,10)$ to the corner $(10,0)$. It is a good line to choose, because it crosses the domain from the free top edge to the loaded right edge and therefore passes through the regions where each boundary condition is active. The next cells evaluate the trained model along a hundred and one evenly spaced points on that line, select the finite element nodes lying on it, and draw four panels, one per field, each carrying the PINN prediction as a solid line and the finite element reference as a dashed black line.
Before looking, predict which of the four panels will show the largest discrepancy. The argument was made above: the loss weights all residuals equally, the equation and stress residuals carry a factor of $E/(1-\nu^2)\approx 55$, and the transverse displacement $v$ is the smallest of the fields, so it exerts the least influence on the objective and has the least incentive to be right.
!gdown 1NVLCxoF9HX-K2wa1Ofk7TaRCdYzktRCY -O /content/lab2_Plate_data.npy
Plate_data = np.load('/content/lab2_Plate_data.npy')
loc = Plate_data[:, 0:2]
w = Plate_data[:, 2:3]
v = Plate_data[:, 3:4]
stress = Plate_data[:, 4:6]
def RESULT(test_model, diag):
diag = torch.tensor(diag, dtype = torch.float32, requires_grad = True).to(device)
pde_w, pde_v = test_model(diag)
pde_disp = np.hstack([pde_w.cpu().detach().numpy(), pde_v.cpu().detach().numpy()])
sig_x, sig_y = check_stress(test_model, diag)
pde_sig = np.hstack([sig_x.cpu().detach().numpy(), sig_y.cpu().detach().numpy()])
return pde_disp, pde_sig
diag_ = np.where(np.isclose(loc[:, 0] + loc[:, 1], 10))[0]
diag_x = np.linspace(0, 10, 101).reshape(-1, 1)
diag_y = -diag_x + 10
diag = np.concatenate((diag_x, diag_y), 1)
FEM = [w[diag_], v[diag_], stress[diag_, 0], stress[diag_, 1]]
disp, sig = RESULT(model, diag)
PINN_pred = [disp[:, 0], disp[:, 1], sig[:, 0], sig[:, 1]]
plt.figure(figsize = (7, 6))
for idx, title in enumerate(titles):
plt.subplot(2, 2, idx + 1)
plt.plot(diag[:, 0], PINN_pred[idx], linewidth = 2, label = 'PINN')
plt.plot(diag[:, 0], FEM[idx], linestyle = '--', c = 'k', linewidth = 3, alpha = 0.5, label = 'FEM')
plt.xlabel('x')
plt.ylabel(title)
plt.legend()
plt.tight_layout()
plt.show()
The results for the $y$-displacement and the $y$-direction stress deviate visibly from the ground truth, while the $x$-direction quantities agree well. This is the outcome predicted above, and it is a genuine limitation rather than a bug: the network is solving the correct equations, and it is solving them in a way that concentrates its accuracy where the loss is most sensitive.
The conventional remedies are available. We could tune the loss weights to raise the influence of the terms involving $v$, increase the number of collocation points, enlarge the network, or train longer. Instead we take a different route, one that is only available to a physics-informed model and that illustrates the central claim of this notebook from the opposite direction.
Assume that observed data is available, such as displacement values obtained from physical sensors placed on the thin plate or extracted from an earlier analysis of a similar component. Section 3 introduced the data term for inverse problems, where it was indispensable. Here the problem is a perfectly well-posed forward problem that needs no data at all, and we add data anyway, to see what it buys.
Adding data constraints for regularization
$$ \begin{align*} \mathcal{L} &= \omega_{\text{data}} \mathcal{L}_{\text{data}} + \omega_{\text{PDE}} \mathcal{L}_{\text{PDE}} + \omega_{\text{BC}} \mathcal{L}_{\text{BC}}\\\\\\ \mathcal{L}_{\text{data}} &= \frac{1}{N_{\text{data}}} \sum \left(u - \hat u \right)^2 \\ \end{align*} $$
The interpretation is worth stating carefully, because it is the reverse of the usual one. In ordinary machine learning the data term is the objective and everything else is regularization; here the physics is the objective and the data is the regularizer. What the data supplies is information about which of the many nearly-equation-satisfying fields is the right one, and it supplies it most usefully in exactly the places where the equation residual is least sensitive, which is where the pure solver was weakest.
Define Observed Data
The observation set is selected by a condition on both coordinates at once: a finite element node is retained when its $x$-coordinate is one of $0, 2, 5, 10$ and its $y$-coordinate is also one of those four values. The retained nodes therefore sit at the sixteen intersections of four vertical lines with four horizontal ones, a coarse rectangular array of sensors rather than a dense sampling of anything. For each of them the two displacement components and the two stress components are retained.
Read the logic of the selection rather than the numbers, because a single-coordinate test would have given something quite different. Testing $x$ alone would retain four complete vertical lines, several hundred nodes; conjoining the two tests reduces that to the handful of points where the two families of lines cross. Sixteen scattered readings against the $2500$ collocation points of Section 6.2 is a very thin data budget, which is what makes the improvement worth measuring.
One honest caveat about the experiment. The observations are drawn from the same finite element solution that will serve as the reference, so the comparison is favourable by construction; real sensor data would be noisy, would not lie on a tidy lattice, and would not include stresses, which are not directly measurable and would have to be inferred from strain gauges.
x_values = [0, 2, 5, 10]
tr_sample = np.where(np.logical_and(np.isin(loc[:, 0], x_values),
np.isin(loc[:, 1], x_values)))[0]
x_sample = loc[tr_sample, :]
w_sample = w[tr_sample, :]
v_sample = v[tr_sample, :]
stress_sample = stress[tr_sample, :]
data = np.concatenate([x_sample, w_sample, v_sample, stress_sample], axis = 1)
plt.figure(figsize = (5, 5))
plt.scatter(x_sample[:, 0], x_sample[:, 1], s = 10)
plt.xlim([-0.5, 10.5])
plt.xlabel('x')
plt.ylabel('y')
plt.show()
Define Data-Driven Loss Term
The data function evaluates the network at the observation points and returns the two displacements and the two stresses, the latter assembled from derivatives exactly as before.
That the stresses are included is the most important detail of Section 6.3 and it is easy to read past. A data term on displacement alone would supervise the network's values. Because the observed stresses are functions of the network's derivatives, the stress part of the data term supervises the network's derivatives directly, at the observation points. Recall the discussion of Section 1.2: closeness in value does not imply closeness in derivative, and it is precisely the derivative-valued field $\sigma_{yy}$ that the pure solver got wrong. Supervising it directly attacks the weakness at its source.
def DATA(model, loc):
w, v = model(loc)
w_x, _ = derivative(w, loc)
_, v_y = derivative(v, loc)
sigma_xx = E/(1 - nu**2)*(w_x + nu*v_y)
sigma_yy = E/(1 - nu**2)*(v_y + nu*w_x)
sigma = torch.cat((sigma_xx, sigma_yy), dim = 1)
return w, v, sigma
Train
A new model of the same class is created and trained for the same ten thousand epochs, with the loss now carrying three groups of terms:
$$ \mathcal{L} = \mathcal{L}_{\text{PDE}} + \mathcal{L}_{\text{BC}} + \mathcal{L}_{\text{data}}, $$
where the data group is itself the sum of three mean squared errors, on $u$, on $v$, and on the pair of stresses. All weights are again $1$, so in the notation of the formula above $\omega_{\text{data}} = \omega_{\text{PDE}} = \omega_{\text{BC}} = 1$. As in Section 6.2 the training cell may be skipped; the trained model is loaded below.
Note that the comparison being set up is a fair one: the same architecture, the same optimizer, the same learning rate, and the same number of epochs, with the data term as the only difference. Whatever improvement appears is attributable to the data and to nothing else.
model = PINN().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr = 1e-3)
epochs = 10001
loc_data = torch.tensor(data[:, 0:2], dtype = torch.float32, requires_grad = True).to(device)
gt_w = torch.tensor(data[:, 2:3], dtype = torch.float32).to(device)
gt_v = torch.tensor(data[:, 3:4], dtype = torch.float32).to(device)
gt_sigma = torch.tensor(data[:, 4:6], dtype = torch.float32).to(device)
for epoch in range(epochs):
# PDE
force_eq_x, force_eq_y = PDE(model, Y_domain)
loss_PDE = torch.mean(force_eq_x**2) + torch.mean(force_eq_y**2)
# BC
bottom_v = BC_bottom(model, Y_bottom)
left_w = BC_left(model, Y_left)
top_sigyy, top_sigxy = BC_top(model, Y_top)
right_sigxx, right_sigxy = BC_right(model, Y_right)
loss_BC = (torch.mean(bottom_v**2) + torch.mean(left_w**2)
+ torch.mean(top_sigyy**2) + torch.mean(top_sigxy**2)
+ torch.mean(right_sigxx**2) + torch.mean(right_sigxy**2))
# DATA
data_w, data_v, data_sigma = DATA(model, loc_data)
loss_Data = torch.mean((data_w - gt_w)**2) + torch.mean((data_v - gt_v)**2) + torch.mean((data_sigma - gt_sigma)**2)
# Total
loss = loss_PDE + loss_BC + loss_Data
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 1000 == 0:
print('Epoch: {} Loss: {:.4f} PDELoss: {:.4f} BCLoss: {:.4f} DATALoss: {:.4f}'.format(
epoch, loss.item(), loss_PDE.item(), loss_BC.item(), loss_Data.item()
))
The printed trace repays a moment's reading. At epoch $0$ the total is $22.6604$, split as $3.2873$ from the equilibrium equations, $6.9588$ from the boundary conditions and $12.4143$ from the data. The data term dominates at initialization for a reason that should now be predictable: an untrained network outputs values near zero, and the measured displacements and stresses are not near zero, so the data residual is essentially the size of the measurements themselves.
By epoch $1000$ the total has fallen to $0.0031$ and by the end of the run to about $10^{-4}$, with the three groups roughly balanced and the equation term the largest of them. The trace is not monotone: at epoch $2000$ the total rises to $0.0178$, almost all of it in the data term, before falling again. Adam is not a descent method in the strict sense; its momentum and its per-parameter scaling can carry it into a worse region, and on a non-convex objective an occasional excursion of this kind is unremarkable so long as the run recovers. It is a useful reminder that a single printed loss is a sample of a trajectory rather than a measurement of quality, and it is one more argument for reading the whole trace rather than the last line.
Results
The trained model is provided so that the outcome can be inspected without waiting. It is loaded and evaluated on the same $100 \times 100$ grid as in Section 6.2.
The four maps are drawn with the same fixed colour limits as before, so they can be laid alongside the earlier figures and the two sets compared panel by panel.
!gdown 1oKqRzA3GvJuz1H98C6Kj_x_4xeLMKALL -O /content/lab2_data_model.pt
model = torch.load('/content/lab2_data_model.pt', map_location = device, weights_only = False)
u_test, v_test = model(test_tensor)
sigma_xx, sigma_yy = check_stress(model, test_tensor)
plot_results(test, [u_test, v_test], color_legend[:2], titles[:2])
plot_results(test, [sigma_xx, sigma_yy], color_legend[2:], titles[2:])
Comparison with FEM data along the diagonal
The decisive plot is the next one. The cell below evaluates the data-trained model along the same anti-diagonal as in Section 6.2 and draws the four panels, this time with the data-trained PINN in green against the same dashed black reference.
Because the two figures are drawn on the same axes and along the same line, the comparison the section exists for is made by looking at this figure and the corresponding figure of Section 6.2 in turn, one panel at a time. In the two $x$-direction panels the pure solver was already close to the reference, and there is little room for the data to help. In the $y$-displacement and $\sigma_{yy}$ panels the pure solver deviated, and the data-trained curve should sit closer to the black line. That is the contribution of the data term, stated as a picture: it does not improve the model uniformly, it improves it where the physics-only loss was insensitive.
disp, sig = RESULT(model, diag)
PINN_data = [disp[:, 0], disp[:, 1], sig[:, 0], sig[:, 1]]
plt.figure(figsize = (7, 6))
for idx, title in enumerate(titles):
plt.subplot(2, 2, idx + 1)
plt.plot(diag[:, 0], PINN_data[idx], c = 'g', linewidth = 2, label = 'PINN + data')
plt.plot(diag[:, 0], FEM[idx], linestyle = '--', c = 'k', linewidth = 3.5, alpha = 0.5, label = 'FEM')
plt.xlabel('x')
plt.ylabel(title)
plt.legend()
plt.tight_layout()
plt.show()
From this example, we can learn two key lessons.
The first is mechanical: you now know how to incorporate observational data into a PINN. It is one more mean squared error term, structurally identical to every other term in the loss, differing only in that the target comes from a measurement rather than from an equation. No new machinery is needed and the change amounts to a handful of lines.
The second is strategic: when data is available, using it improves a PINN, and doing so is generally to be recommended. The reason is worth stating in the general form, because it is the argument of the whole notebook read in reverse. A physics-only model is constrained everywhere but only through the equation, which is insensitive to some features of the solution; a data-only model is constrained sharply but only at the measurement points; combining them constrains the model sharply where the sensors are and everywhere else through the physics. Neither source of information is redundant, because they are informative about different things.
It is worth adding the caveat that this laboratory could not show. Here the data was drawn from the reference solution and was therefore exactly consistent with the equations. Real measurements are noisy and are not exactly consistent with any model, so a large $\omega_{\text{data}}$ will drag the field towards the noise and a small one will ignore the sensors. The weight then becomes a genuine bias-variance decision familiar from regression, rather than the free choice it appears to be here.
Key takeaways
- A data-driven model fails in identifiable ways when data are scarce: physically impossible outputs, spurious attributions, poor extrapolation and no mechanism to inspect. Engineering problems have scarce data and known equations, and physics-informed learning puts the equations where the data would have been, in the objective. The universal approximation theorem licenses the attempt by saying that the solution is representable; it says nothing about whether a given optimisation will find it.
- The physics-informed loss is $L = L_{\text{PDE}} + \lambda_I L_{\text{IC}} + \lambda_B L_{\text{BC}}$, mean squared residuals of the equation at collocation points and of the side conditions on the boundary, each obtained by differentiating the network with respect to its input coordinates. Observations enter as one more mean squared term, and unknown coefficients become tensors handed to the optimizer alongside the weights; together they turn an under-determined inverse problem into an over-determined least-squares one.
- In PyTorch the whole method is three lines: the coordinate tensor carries
requires_grad = True,torch.autograd.grad(u, t, grad_outputs = torch.ones_like(u), create_graph = True)returns $du_\theta/dt$ at every point at once, andcreate_graph = Truelets the loss built from that derivative be differentiated again bybackward(). Activations must be smooth ($\tanh$), and the output layer is linear so that the range of the solution is reachable. - The loss at initialisation can be predicted from the equation: for $du/dt = \cos 2\pi t$ with $u(0) = 1$ it is about $0.5 + 1 = 1.5$, and the run confirms it. Where that run ends, however, is decided by choices that look like housekeeping, resampled versus fixed collocation points and the initialisation scheme, which is the theme Regularization, MAP Estimation, and Loss Weighting and PINN: Limitations and Extensions develop.
- The beam and the plate show the method on a fourth-order equation and on a coupled two-dimensional system. Unit weighting of residuals that carry a factor $E/(1-\nu^2)$ neglects the small transverse displacement, and a handful of measured values repair exactly that field. Physics-Informed Neural Networks: Labs carries the same idea into inverse problems, where data and prior knowledge identify unknown constants and an unknown edge.
Looking ahead. A PINN solves one problem, for one set of coefficients and boundary data, and must be retrained for the next. Deep Neural Operators and Fourier Neural Operator learn the map from coefficients or forcing to the solution itself, so that a new instance costs one forward pass; Embedding Physical Knowledge into Neural Networks broadens the view from equations in the loss to every way domain knowledge can be built into a network.