Differential equations: ODEs, PDEs, and the forward solve.
In mathematics, we often hear about differential equations - equations that relate a function to its own rate of change. Within differential equations, you’ll find a distinction between Ordinary Differential Equations (ODEs) and Partial Differential Equations (PDEs), but you may find yourself asking: What is the difference between the two?
Well, at its simplest (and it may sound silly to say), the difference between and ODE and a PDE is that an ODE is the number of independent variables.
Very simply, an ODE relates a function to its own rate of change via a single independent variable and a full derivative. Perhaps one of the simplest examples of this is the function: $u\left(x\right)=\frac{d u\left(x\right)}{dx}$, where $u\left(x\right) = e^x$. Another example is an oscillating mass on a spring, given by: $m\frac{d^2x}{dt^2} + kx = 0$.
A PDE relates a function to its own rate of change via multiple independent variables using partial derivates. For this, think of the heat equation, which depends on location and time: $\frac{\partial u \left(x,t\right)}{\partial t} = \alpha \frac{\partial^2 u\left(x,t\right)}{\partial x ^2}$. In this case, $u\left(x,t\right)$ is a function of both x and t, therefore, we must take the partial derivative of the function u, meaning we hold the other independent variables as constant.
ODEs essentially model systems that evolve in a single dimension while PDEs model systems that evolve in multiple dimensions simultaneously. PDEs therefore are often (but not always) more difficult to solve, and require numerical methods such as finite differences, finite elements, finite volumes, spectral methods, etc. or they may only have solutions for specific boundary or initial conditions.
Forward solving PDEs
Let's consider a general form of a PDE: $$ \mathcal{L}_a u\left(\bf{x}\right)=f\left(\bf{x}\right),$$ subject to $b\left(\bf{x}\right) = u\left(\bf{x}\right)$. $\mathcal{L}$ is a linear differential operator, with $a$ being some coefficients that relate to the physical setup of the system and can vary across the system, $u$ being the response that we want to find and $f$ being the forcing function imposed on the system. The solutions to these equations are then defined by the interplay being $u$, $a$, and $f$.
We want to find the solution, which is how $u$ responds given the operator $\mathcal{L}$, the coefficients $a \left( x \right)$, and the forcing function $f\left(x\right)$. Put mathematically, we are looking for: $$ u\left(x\right) = \mathcal{L}^{-1}_a f\left(x\right).$$ This is known as the forward problem.
The Poisson Equation
Let's work through a concrete example of this (the well-known Poisson equation) with the following defined operator: $$ -\frac{\partial}{\partial x}\left(a(x)\frac{\partial u(x)}{\partial x}\right) = f(x), \qquad u(0)=u(1)=0. $$
Here, $\mathcal{L}=-\frac{\partial }{\partial x}\left( a\left(x\right) \frac{\partial }{\partial x}\right)$.
Next, we consider the simplest case of constant coefficients, e.g. the coefficients $a$ are the same across the domain: $a\in A = 1$. First, we expand the derivative in the operator via the product rule: $$\mathcal{L}=-\frac{\partial }{\partial x}\left( a\left(x\right) \frac{\partial }{\partial x}\right) = - a\left(x\right)\frac{\partial ^2u}{\partial x^2} - \frac{\partial a\left(x\right)}{\partial x}\frac{\partial u}{\partial x}.$$
In this scenario with $a=1$ everywhere on the domain, the right hand side of the expanded operator falls away to 0 (because $\frac{\partial }{\partial x}\left(1\right)=0$), leaving us with:
$$ -\frac{\partial ^2u}{\partial x^2} = f\left(x\right) ; u\left(0\right) = u\left(1\right) = 0.$$
Finite differences
Now we have the second derivative of $u$ set equal to our forcing function, $f$. To solve this, we use finite differences, where we uniformly discretise the domain into N points $x_i$, with forcings $f_i$ and responses $u_i$. We then approximate the second derivative as: $$ \frac{\partial ^2 u\left(x_i\right)}{\partial x^2}\approx \frac{-u_{i-1} + 2u_i - u_{i+1}}{h^2}=f_i, $$ where $h=\frac{1}{N+1}$ is the spacing of the uniform grid in $x$.
This is the finite differences solution. We can write this as a linear system equations: $\mathbf{A}u=f$, and this can be solved very efficiently using the forward elimination and backwards substitution method.
Green's functions
Instead of doing a full numerical solve each time, we might want to explore other methods. One such method is the Green's function. At its core, Green's functions are useful because they allow us to solve a differential equation as an integral equation. We are essentially trying to find the inverse of the differential operator: $\mathcal{G}~\mathcal{L}^{-1}$, but how do we get there?
Given our earlier equation: $\mathcal{L}u\left(x\right)=f\left(x\right)$, a Green's function $\mathcal{G}\left(x, y\right)$ is defined as: $$ \mathcal{L}\mathcal{G}\left(x, y\right)= \delta\left(x-y\right),$$ where $\delta\left(x-y\right)$ is a Dirac delta function localized at $y$. Thus, $\mathcal{G}\left(x,y\right)$ is the solution of our PDE produced by an infinitetesimal impulse by a point source at $y$.
With a little bit of re-ordering of our definition of a Green's function above we can see that: $$ \mathcal{G}\left(x,y\right) = \mathcal{L}^{-1} \delta\left(x-y\right). $$
Now, because our equation is linear (as defined earlier), any arbitrary source $f\left(x\right)$ can be decomposed into a continuous sum of the contribution of our forcing at discrete point sources: $$ f\left(x\right) = \int f\left(y\right) \delta\left(x-y\right)dy.$$
We then have from our original definition: $$\mathcal{L} u\left(x\right) = f\left(x\right) = \int f\left(y\right) \delta\left(x-y\right)dy, $$
because $\mathcal{L}$ is linear and the integral is over $y$ with fixed $x$, we can quickly write: $$\mathcal{L}u\left(x\right) = f\left(x\right) = \int f\left(y\right) \mathcal{L}^{-1} \delta\left(x-y\right)dy = \int f\left(y\right) \mathcal{G}\left(x,y\right)dy. $$
Why is this so handy? Well, once you numerically compute the operator $\mathcal{G}$ once, you have it! he Green’s function encodes the intrinsic response properties of the system (geometry, coefficients $ a(x) $, boundary conditions). The forcing $ f $ is just the external input. Because the system is linear, the mapping from input to output is fixed — and that fixed mapping is precisely what $ \mathcal{G} $ represents. This is why you only construct $ \mathcal{G} $ once.
Linear algebra analogy
Consider the finite-difference version from earlier.
You discretize the domain and get a matrix equation: $$\mathbf{A} \mathbf{u} = \mathbf{f}$$ where $ \mathbf{A} $ is the discrete version of $ \mathcal{L} $.
Solving this once gives you $ \mathbf{u} = \mathbf{A}^{-1} \mathbf{f} $. The inverse matrix $ \mathbf{A}^{-1} $ is universal: you compute it (or factor it) once, and then for any new right-hand side $ \mathbf{f} $, you just do a matrix-vector multiply: $ \mathbf{u} = \mathbf{A}^{-1} \mathbf{f} $.
The Green’s function $ \mathcal{G}(x,y) $ is exactly the continuous analogue of the inverse matrix $ \mathbf{A}^{-1} $. The entry $ \mathcal{G}(x_i, y_j) $ tells you the response at grid point $ i $ due to a unit force at grid point $ j $. The integral $ u(x) = \int \mathcal{G}(x,y) f(y) , dy $ is the continuous version of the matrix-vector product. Just as one inverse matrix works for every possible $ \mathbf{f} $, one Green’s function works for every possible $ f(x) $.
Alternative solvers
Just as a comment, the finite diferences (FD) method isn't the only means of numerically solving PDEs. Two other helpful methods are the Finite Elements (FE) and Finite Volumes (FV) methods, which both have different assumptions.
Finite Element Method
Rather than approximating pointwise-derivatives (as is done in the FD method), the finite element (FE) method asks a slightly different question: instead of enforcing $\mathcal{L}_a u = f$ exactly at every grid point, can we enforce it on average against a family of test functions?
To see this, take our equation and multiply both sides by an arbitrary smooth function $v(x)$ (this is the test function we were talking about), also with Dirilicht boundary conditions ($v(0) = v(1) = 0$), and integrate over the domain:
$$\int_0^1 -\frac{\partial }{\partial x}\left(a(x)\frac{\partial u}{\partial x}\right) v(x) dx = \int_0^1 f(x) v(x) dx.$$
Integrating the left-hand side by parts shifts one derivative off $u$ and onto $v$, giving the weak form${^[1]}$:
$$\int_0^1 a(x) u'(x) v'(x) dx = \int_0^1 f(x) v(x) dx.$$
This is a weaker requirement than the original PDE — $u$ only needs to be differentiable once rather than twice — which turns out to be a feature, not a bug. We then partition $[0,1]$ into small subintervals (elements) and approximate $u$ as a sum of piecewise polynomial basis functions $\phi_j(x)$, each supported over just a few adjacent elements:
$$u(x) \approx \sum_j u_j \phi_j(x).$$
Substituting this into the weak form and choosing $v = \phi_i$ for each basis function in turn produces a linear system $\mathbf{K}\mathbf{u} = \mathbf{f}$, where $K_{ij} = \int_0^1 a(x) \phi_i'(x) \phi_j'(x) dx$ is the stiffness matrix$^{[2]}$. For our 1D Poisson problem with constant $a$, FEM and FD on the same uniform grid produce identical systems — so the advantage here is not immediately obvious. It becomes clear once $a(x)$ varies across the domain: the integral formulation handles spatially heterogeneous coefficients naturally, and the element-by-element assembly extends cleanly to irregular geometries in higher dimensions where a uniform grid is no longer practical.
Finite Volume Method
The finite volume (FV) method takes yet another perspective, this time rooted directly in the conservation structure of the equation. Rather than approximating the PDE at points or in a weighted-integral sense, we integrate it exactly over small control volumes centred on each grid point.
Integrating $-\frac{d}{dx}\left(a(x)\frac{du}{dx}\right) = f(x)$ over the control volume $[x_{i-1/2}, x_{i+1/2}]$ gives us:
$$-\left\lbrack a(x)\frac{du}{dx}\right\rbrack_{x_{i-1/2}}^{x_{i+1/2}} = \int_{x_{i-1/2}}^{x_{i+1/2}} f(x), dx.$$
The left-hand side is simply the net flux leaving the control volume — the difference between what flows out of the right face and what flows in from the left. Approximating each face flux and the right-hand side integral gives:
$$-a_{i+1/2}\frac{u_{i+1} - u_i}{h} + a_{i-1/2}\frac{u_i - u_{i-1}}{h} = f_i h,$$
which again yields a (tridiagonal) linear system. The crucial property here is that the flux leaving one cell is exactly the flux entering its neighbour — conservation holds at the discrete level, by construction, regardless of grid resolution. This makes FV method the natural choice when $a(x)$ has sharp spatial jumps, such as the layered permeabilities that occur in Darcy flows. In those cases, pointwise approximations of derivatives can behave poorly at discontinuities, whereas flux-based accounting remains well-defined.
Concluding comments
FD, FE, and FV all provide forward solves of PDEs using different assumptions. In the simple case we have been working in (a linear operator, no time dependence, constant coefficients) this is essentially the whole story. Although these methods differ in how they handle geometry, discontinuities, and conservation, the underlying structure is the same.
Things get considerably more involved when any of these assumptions are relaxed.
Non-linear cases. Our equation is linear because $\mathcal{L}_a$ acts linearly on $u$ — doubling $f$ doubles $u$. If instead the operator depends on $u$ itself, for example $-\frac{d}{dx}\left(a(u)\frac{du}{dx}\right) = f$, linearity is lost. The assembled matrix $\mathbf{A}$ now depends on the very solution we are trying to find, so we cannot simply solve $\mathbf{A}\mathbf{u} = \mathbf{f}$ in one step. Instead we iterate: guess a solution $u^{(k)}$, assemble $\mathbf{A}(u^{(k)})$, solve for a better $u^{(k+1)}$, and repeat until convergence. Newton's method and Picard iteration are the standard approaches. Here, convergence is no longer guaranteed, and the cost per solve multiplies with the number of iterations required.
Time-dependent cases. Adding a time derivative turns our elliptic equation into a parabolic one. The heat equation is the canonical example:
$$\frac{\partial u}{\partial t} - \frac{\partial}{\partial x}\left(a(x)\frac{\partial u}{\partial x}\right) = f(x, t).$$
The spatial discretisation (FD, FE, or FV) still produces a system $\mathbf{A}\mathbf{u} = \mathbf{f}$, but now $\mathbf{u}$ and $\mathbf{f}$ evolve in time and we must additionally choose a time integration scheme. Explicit schemes (e.g. forward Euler) are simple but impose a stability constraint on the time steps (e.g., the CFL condition) which can force very small steps (which in turn means many steps) when the spatial grid is fine. Implicit schemes (e.g. backward Euler, Crank-Nicolson) are stable for much larger time steps but require solving a linear system at every step. The choice between them becomes a trade-off between cost per step and the maximum allowable step size for whatever problem you're considering$^{[3]}$.
Variable coefficients. The coefficients of the operator $a(x)$ can vary across the domain we're working on. This is actually what happens in most cases. In this situation, the spatial structure of $\mathbf{A}$ changes and its entries depend on local values of $a$, so the matrix must be reassembled for each new coefficient field. For a single solve this is manageable. The issues come when we need solutions for many different $a(x)$. In those settings, the cost of repeated assembly and factorisation becomes prohibitive. This is precisely the regime where the Green's function and, more recently, machine learning methods become attractive: rather than solving afresh for every $a$, we want to learn the map $a \mapsto u$ directly. This will be the subject of a separate post.
$^{[1]}$: Whereas the strong form of an equation must hold at all points in the domain, the weak form of an equation must only hold on average.
$^{[2]}$: Stiffness refers to how strongly one point influences its neighbours.
$^{[3]}$: Sounding a bit more art-than-a-science now, eh?
