Let $\Omega \subset \mathbb{R}^n$ be an open, bounded domain with boundary $\partial\Omega$. We seek $u : \bar{\Omega} \to \mathbb{R}$ satisfying:

where $\mathcal{L}$ is a differential operator, $f$ is a known forcing function, and $g$ is the boundary data. A canonical example is the Poisson equation, where $\mathcal{L} = -\nabla^2$.

Since we cannot generally find $u$ exactly, we build an approximation $\tilde{u} \approx u$. Define the residual:

Note that $R = 0$ everywhere if and only if $\tilde{u}$ is the exact solution. For any other $\tilde{u}$, we have $R \not\equiv 0$.

Before choosing a method, we must specify where $\tilde{u}$ and the test functions live. This is not just formalism — it determines whether the method is well-posed and what regularity the solution needs.

We define two spaces: the trial space $V$ (where we seek $\tilde{u}$, which must satisfy the boundary conditions) and the test space $W$ (where weight functions $w$ live).

The general weighted residual statement is:

Answer: By the Dirac delta sifting property:

This enforces the residual to be exactly zero at each collocation point $x_i$. Geometrically, we are pinning the residual to zero at a finite set of chosen points.

We approximate $u$ using a finite linear combination of basis functions $\{\phi_j\}_{j=1}^N$:

where $c_j \in \mathbb{R}$ are the unknown coefficients we want to determine, and $\phi_j \in V$ are prescribed basis functions satisfying the boundary conditions.

Substituting into (1.3), the residual becomes:

Partition $\Omega$ into $N$ non-overlapping subdomains $\Omega_1, \Omega_2, \ldots, \Omega_N$ such that $\bigcup_k \Omega_k = \Omega$.

The test functions are characteristic (indicator) functions:

The weighted residual condition (1.4) for each subdomain $\Omega_k$ becomes:

Substituting (2.2):

This is a linear system $\mathbf{A}\mathbf{c} = \mathbf{b}$ with:

Choose $N$ distinct collocation points $\{x_1, x_2, \ldots, x_N\} \subset \Omega$.

The test functions are Dirac deltas:

From the sifting property, condition (1.4) gives:

Substituting (2.2):

Again a linear system $\mathbf{A}\mathbf{c} = \mathbf{b}$ with:

Here the test space equals the trial space: $W = V$. We take $w = \phi_i$ for each $i = 1, \ldots, N$. Condition (1.4) becomes:

Substituting (2.2):

Again a linear system $\mathbf{A}\mathbf{c} = \mathbf{b}$ with:

For self-adjoint operators (e.g., $\mathcal{L} = -\nabla^2$), one typically integrates by parts to reduce regularity requirements. This gives the weak form, which only requires $\tilde{u} \in H^1_0(\Omega)$ instead of $C^2$:

Answer:

(a) In the subregion method, $A_{kj}$ is an integral of $\mathcal{L}[\phi_j]$ over a subdomain $\Omega_k$. In the collocation method, $A_{kj}$ is a pointwise evaluation $\mathcal{L}[\phi_j](x_k)$ — no integration. In the Galerkin method, $A_{ij}$ is a weighted integral of $\mathcal{L}[\phi_j]$ against the test function $\phi_i$ over the entire domain $\Omega$.

(b) The collocation method requires the highest regularity. Since $\mathcal{L}[\tilde{u}]$ must be evaluated pointwise, all derivatives in $\mathcal{L}$ must exist classically. For $\mathcal{L} = -\nabla^2$, this demands $\tilde{u} \in C^2(\Omega)$. The Galerkin weak form, by contrast, only needs $\tilde{u} \in H^1_0(\Omega)$.

In Section 2, our approximation was:

This is linear in the parameters $c_j$. The basis functions $\phi_j$ are fixed and chosen a priori (polynomials, piecewise linear hat functions, trigonometric functions, etc.).

Two consequences:

(a) The weighted residual equation always produces a linear system $\mathbf{A}\mathbf{c} = \mathbf{b}$.

(b) The quality of the approximation depends entirely on how well the fixed basis $\{\phi_j\}$ can represent $u$ — this is the approximation power of the space $V^h = \text{span}\{\phi_1, \ldots, \phi_N\}$.

Now replace the linear combination with a neural network. Define:

where $\theta \in \mathbb{R}^p$ collects all weights and biases, and $\mathcal{N}$ is the network architecture (e.g., a feedforward multilayer perceptron).

For a single hidden layer with $m$ neurons and activation $\sigma$:

where $\theta = \{\alpha_k,\, w_{ki},\, b_k,\, \beta\}$.

The set of all functions representable by a given architecture forms:

This is not a linear subspace — it is a nonlinear manifold in function space.

However, the Universal Approximation Theorem (Cybenko 1989, Hornik 1991) guarantees:

If the activation function $\sigma$ is smooth (e.g., $\sigma = \tanh$, sigmoid, softplus), then:

because $u_\theta$ is a composition of smooth maps. All derivatives exist and can be computed via automatic differentiation.

Answer:

(a) No. Since $u_\theta$ is nonlinear in $\theta$, the residual $R(x;\,u_\theta) = \mathcal{L}[u_\theta](x) - f(x)$ is a nonlinear function of $\theta$. We cannot factor it as $\mathbf{A}\boldsymbol{\theta} - \mathbf{b}$. The linearity that gave us equation (2.2) relied on $\tilde{u}$ being a linear combination of fixed basis functions — that structure is gone.

(b) We replace the linear solve with optimization. We define a scalar loss function $\mathcal{J}(\theta)$ that measures how large the residual is, and minimize it via gradient-based methods (gradient descent, Adam, L-BFGS, etc.) using gradients computed by automatic differentiation.

PINNs combine the neural network approximation of Section 3 with the collocation philosophy of Section 2.3. The key move: instead of solving $R(x_k;\,u_\theta) = 0$ as a nonlinear system, we minimize the squared residual over a set of collocation points.

Consider the boundary value problem from Section 1:

We represent the solution by a neural network $u_\theta(x) = \mathcal{N}(x;\,\theta)$ and define the PDE residual at any point $x \in \Omega$:

and the boundary residual at any point $x \in \partial\Omega$:

Sample two sets of collocation points:

(i) Interior points: $\{x_r^{(i)}\}_{i=1}^{N_r} \subset \Omega$ (for the PDE residual)

(ii) Boundary points: $\{x_b^{(i)}\}_{i=1}^{N_b} \subset \partial\Omega$ (for the boundary condition)

Define the component losses:

The total PINN loss is:

where $\lambda_r, \lambda_b > 0$ are weighting hyperparameters that balance PDE fidelity against boundary satisfaction.

The PINN training problem is:

This is a nonconvex optimization problem in $\mathbb{R}^p$, where $p$ can be $10^4$–$10^7$ depending on the network size. We solve it iteratively via gradient descent:

where $\eta > 0$ is the learning rate. In practice, variants such as Adam or L-BFGS are used.

The gradient $\nabla_\theta \mathcal{J}$ requires differentiating through the PDE operator $\mathcal{L}$. By the chain rule, for a single collocation point $x_r$:

The term $\nabla_\theta[\mathcal{L}[u_\theta](x_r)]$ involves second-order derivatives: first, spatial derivatives of $u_\theta$ (to evaluate $\mathcal{L}$), and then parameter derivatives of the result (to compute the gradient). This is handled by automatic differentiation — specifically, by differentiating through the computational graph of $u_\theta$.

One of the strengths of PINNs is the ability to fuse physics (the PDE) with observational data. Suppose we have $N_d$ measurements $\{(x_d^{(i)},\, u_d^{(i)})\}_{i=1}^{N_d}$ of the true solution. We add a data-fidelity term:

The full PINN loss becomes:

Recall from Section 3.3 that the neural network hypothesis class $\mathcal{F}_\theta$ is a nonlinear manifold, not a linear subspace. This has important consequences:

(a) No Céa's lemma: In the Galerkin/FEM setting, Céa's lemma guarantees that the discrete solution is the best approximation in $V^h$ (up to a constant). No such result exists for PINNs — minimizing $\mathcal{J}(\theta)$ does not guarantee we find the closest function in $\mathcal{F}_\theta$ to the true solution.

(b) No guaranteed convergence: The loss landscape is nonconvex. Gradient descent can get trapped in local minima or saddle points. The quality of the solution depends on initialization, architecture, optimizer, and learning rate schedule.

(c) Regularity is controlled by architecture: With smooth activations ($\tanh$, $\sin$), $u_\theta \in C^\infty$. But the effective approximation quality depends on width, depth, and the spectral bias of the network (tendency to learn low-frequency components first).

Answer:

(a) With $\lambda_r = 0$, the PDE is completely absent from the loss. What remains is boundary matching ($\mathcal{J}_b$) and data fitting ($\mathcal{J}_d$). This is pure data fitting / interpolation — the network learns to pass through the data and satisfy boundary values, but has no knowledge of the governing physics. It is no longer "physics-informed."

(b) With $\lambda_d = 0$, we have $\mathcal{J} = \lambda_r\,\mathcal{J}_r + \lambda_b\,\mathcal{J}_b$, which is the same structure as classical collocation but with two key differences: (i) the trial space is a nonlinear neural network manifold $\mathcal{F}_\theta$ instead of a linear span, and (ii) enforcement is soft (minimize $|R|^2$) rather than hard ($R = 0$ exactly). This is the "vanilla" PINN — a neural-network collocation method with least-squares residual minimization.

In Sections 1–4, we always solved one instance of a PDE: given a specific $f$ and $g$, find the corresponding $u$. If $f$ changes, we must retrain the entire PINN from scratch.

In many applications — design optimization, uncertainty quantification, real-time control — we need to solve the same PDE for many different inputs $f$. The natural mathematical object is the solution operator:

The operator $\mathcal{G}$ maps an input function $f$ (the forcing) to an output function $u$ (the solution). It acts between function spaces, not between finite-dimensional vectors.

Let $\mathcal{A}$ and $\mathcal{U}$ be Banach spaces of functions (e.g., $\mathcal{A} = L^2(\Omega)$ for the input forcing, $\mathcal{U} = H^1_0(\Omega)$ for the solution). The true solution operator $\mathcal{G}: \mathcal{A} \to \mathcal{U}$ is typically nonlinear and defined implicitly through the PDE.

We seek a parametric approximation:

Given a training dataset of input–output pairs $\{(f^{(i)},\, u^{(i)})\}_{i=1}^{N_{\text{train}}}$ (generated by a classical solver), we train by minimizing:

In practice, the norm $\|\cdot\|_{\mathcal{U}}$ is approximated by evaluating on a discrete grid:

The theoretical foundation for neural operators comes from a generalization of the universal approximation theorem to operators.

The Deep Operator Network (DeepONet, Lu et al. 2021) is a direct realization of Theorem 5.1. It decomposes the operator into two sub-networks:

Branch network $\mathbf{b}_\theta$: encodes the input function $f$. The input function is evaluated at $m$ fixed sensor locations $\{\xi_1, \ldots, \xi_m\} \subset \Omega$:

Trunk network $\mathbf{t}_\theta$: encodes the output query location $y \in \Omega$.

Both networks output vectors of the same dimension $q$. The operator evaluation is:

or in compact notation:

The Fourier Neural Operator (Li et al. 2021) takes a different approach: instead of separating branch and trunk, it processes the input function through layers that mix information in Fourier space.

Let $v^{(\ell)}(x)$ be the hidden representation at layer $\ell$. Each FNO layer applies:

where $W^{(\ell)}$ is a pointwise linear transformation and $\mathcal{K}^{(\ell)}$ is an integral kernel operator implemented via the Fourier transform:

Here $\mathcal{F}$ and $\mathcal{F}^{-1}$ are the Fourier and inverse Fourier transforms (implemented via FFT), and $R^{(\ell)} \in \mathbb{C}^{d_v \times d_v \times k_{\max}}$ is a learnable tensor of Fourier coefficients, truncated at $k_{\max}$ modes.

A natural question: can we combine the operator learning framework with the physics-informed loss from PINNs? Yes — this gives physics-informed neural operators (PINO, Li et al. 2024).

The PINO loss combines data-driven and physics-driven terms:

Answer:

(a) In the classical approximation (2.1), $\tilde{u}(x) = \sum_{j=1}^N c_j\,\phi_j(x)$: the basis functions $\phi_j$ are fixed (chosen a priori — e.g., polynomials, hat functions), and the coefficients $c_j$ are fixed scalars determined by solving the linear system. Both are independent of the input $f$; if $f$ changes, we must solve the system again to get new $c_j$.

In DeepONet (5.7), $\mathcal{G}_\theta(f)(y) = \sum_{k=1}^q b_k(f)\,t_k(y)$: the trunk outputs $t_k(y)$ play the role of basis functions, but they are learned (not chosen a priori). The branch outputs $b_k(f)$ play the role of coefficients, but they are functions of the input $f$ — so the same trained network instantly produces new coefficients for any new forcing. In summary:

(b) FNO requires a uniform grid because it relies on the Fast Fourier Transform (FFT) to implement the spectral convolution (5.10). The FFT algorithm assumes uniformly spaced samples — this is a fundamental requirement of the discrete Fourier transform, which exploits the periodicity and uniform spacing structure to achieve $O(M \log M)$ complexity. On a non-uniform grid, the FFT cannot be applied directly (one would need the NUFFT, which is more expensive and approximate).

DeepONet's trunk network, by contrast, takes the query coordinate $y$ as a continuous input to a standard neural network — it can evaluate $t_k(y)$ at any point $y \in \Omega$, regardless of how the input $f$ was discretized. This gives the trunk full resolution invariance on the output side.