import matplotlib.pyplot as plt
import numpy as np
A Soft Introduction to Partial Differential Equations (in Python)¶
What are Partial Differential Equations?¶
Partial Differential Equations (PDEs) are equations with more than one independent variables, $x, y, ....$. They consist of a dependent variable, $u$, which is some unknown function that depends on these independent variables, $u(x, y, ...)$, as well as, partial derivatives of the unknown function with respect to these independent variables, $\frac{\partial u}{\partial x} = u_x$. A PDE can be thought of as a function that relates the independent variables $x, y, ....$, dependent variable $u$, and the partial derivatives of $u$.
The most general form of a first-order PDE with two independent variables $(x, y)$ is: \begin{equation} \tag{1} F(x, y, u(x, y), u_{x}(x, y), u_{y}(x, y)) = F(x, y, u, u_{x}, u_{y}) = 0 \end{equation}
The order of a PDE is the highest derivative that appears. For example, a second-order PDE may be written as: \begin{equation} \tag{2} F(x, y, u, u_{x}, u_{y}, u_{xx}, u_{xy}, u_{yy}) = 0 \end{equation}
A solution is the function, $u(x, y, ...)$, that satisfies the PDE exactly (or to some tolerance) for each of the $x, y, ...$ variables.
Initial and Boundary Conditions¶
PDEs are not unique, that is, there exists multiple solutions $u(x, y, ...)$ for any multiple $x, y, ...$ variables. To obtain a single solution, we must impose some auxiliary condition. These conditions are common in many science and engineering applications, and are referred to as initial, (ICs), and boundary, (BCs), conditions.
An IC describes the beginning state(s), $\phi(\mathbf{x})$, of the system at a time $t_0$. $\phi(\mathbf{x}) = \phi(x, y, ...)$ is a function that depends only on the spatial variables $(x, y, ...)$. For a heat flow problem, this represents the initial temperature of the system. It can be written as: \begin{equation} \tag{3} u(x, t_0) = \phi(\mathbf{x}) \end{equation} There may be one or more IC describing the states of the system at $t_0$, such as initial position and initial velocity.
PDEs in science and engineering are only valid within some domain, $D$. For a one-dimensional problem, $D$ can be defined as the interval $0 < x < l$ (a domain with only two endpoints). These endpoints are the boundaries of our domain and in order to solve our PDE, some BC must be specified for each. Three most common BCs commonly used are:
- Dirichlet condition, (D), whereby some value for $u$ is specified,
- Neumann condition, (N), whereby the normal derivative, $\frac{\partial u}{\partial n}$, is specified,
- Robin condition, (R), which is a combination of the Dirichlet and Neumann conditions, $u + \frac{\partial u}{\partial n}$ is specified.
For a one-dimensional problem, the boundaries can be: \begin{equation} \tag{D} u(0, t) = g(t) \quad u(l, t) = h(t) \end{equation}
\begin{equation} \tag{N} \frac{\partial u}{\partial n}(0, t) = g(t) \quad \frac{\partial u}{\partial n}(l, t) = h(t) \end{equation} where $g$ and $h$ are some known functions defined on the boundaries.
Numerical Solutions of PDEs¶
Analytical (exact) solutions of PDEs are often not possible, requiring the use of numerical methods to approximate the solutions of PDEs. Several numerical methods exist, such as the Finite Difference Method, the Finite Element Method, the Finite Volume Method and the Spectral Method. Each method has its own advantages and disadvantages, and choice is often problem dependent. We introduce the Finite Difference Method for numerical solutions, and refer eager readers to Partial Differential Equations: An Introduction for further details.
Finite differences is the most common method used for approximating solutions to PDEs, and involves replacing the derivatives with difference approximations. Consider a function $u(x)$ and a mesh size of $\Delta x$. A solution for $u$ on this mesh can then be defined as $u_j = u(j \cdot \Delta x)$.
nodes = np.arange(start=0, stop=9)
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 2))
fig.suptitle(r"$u(x) = u_j$")
ax.set_xlabel("j")
ax.set_yticklabels([]) # Turn off
ax.plot(nodes, np.zeros_like(nodes), "o-")
By taking the Taylor expansion of $u$ in terms of $x$: \begin{equation} \tag{PX} u(x + \Delta x) = u(x) + u_{x}(x) \cdot \Delta x + u_{xx}(x) \cdot \frac{(\Delta x)^2}{2} + O(\Delta x)^3 \end{equation}
\begin{equation} \tag{NX} u(x - \Delta x) = u(x) - u_{x}(x) \cdot \Delta x + u_{xx}(x) \cdot \frac{(\Delta x)^2}{2} + O(\Delta x)^3 \end{equation}
By discarding second and higher order terms, we can approximate the first derivative, $u_{x}$, by three difference methods - Forward difference (FD) using equation (PX), Backward difference (BD) using equation (NX), and Central difference (CD) using both equations (PX) and (NX):
\begin{equation} \tag{FD} u_{x}(x) = \frac{u(x + \Delta x) - u(x)}{\Delta x} + O(\Delta x) \end{equation}
\begin{equation} \tag{BD} u_{x}(x) = \frac{u(x) - u(x - \Delta x)}{\Delta x} + O(\Delta x) \end{equation}
\begin{equation} \tag{CD} u_{x}(x) = \frac{u(x + \Delta x) - u(x - \Delta x)}{2 \cdot \Delta x} + O(\Delta x)^2 \end{equation}
For functions of two or more variable, $u(x, y, ...)$, a mesh size for each variable must be chosen: $(\Delta x, \Delta y, ...)$. Consider a function, $u(x, y, t)$, with mesh sizes, $\Delta x, \Delta y, \Delta t$, then a point on the mesh can be defined as: $$ u(i \cdot \Delta x, j \cdot \Delta y, n \cdot \Delta t) = u_{i, j}^{n} $$ where $n$ is a superscript for $t$.
We can then approximate the derivatives of $u$ for any of the variables $(x, y, t)$ as above. For example, the forward difference approximations are given as: $$ \frac{\partial u}{\partial t}(i \cdot \Delta x, j \cdot \Delta y, n \cdot \Delta t) = \frac{u_{i, j}^{n + 1} - u_{i, j}^{n}}{\Delta t} \quad or \quad \frac{\partial u}{\partial x}(i \cdot \Delta x, j \cdot \Delta y, n \cdot \Delta t) = \frac{u_{i + 1, j}^{n} - u_{i, j}^{n}}{\Delta x} $$
fig, ax = plt.subplots(
nrows=1, ncols=1, figsize=(8, 6), subplot_kw={"projection": "3d"}
)
# Define function
x = np.linspace(0, 10, 40)
y = np.linspace(0, 10, 40)
X, Y = np.meshgrid(x, y)
U = 4.0 * np.sin(np.sqrt(X**2 + Y**2))
surf = ax.plot_surface(
X,
Y,
U,
rstride=1,
cstride=1,
cmap="plasma",
edgecolor="k",
linewidth=0.005,
antialiased=False,
)
fig.colorbar(surf, shrink=0.5, aspect=10)
ax.set_title(r"$u(x, y, t) = u_{i, j}^n$")
References¶
- Strauss, W.A., 2007. Partial differential equations: An introduction. John Wiley & Sons.
- Partial Differential Equations in Python
- Taylor Series