Guide to partial differential equations

Note: This guide assumes prior knowledge of multivariable and some vector calculus, as well as basic knowledge of vectors and matrices. See Guide to multivariable and vector calculus if you need a refresher/introduction to these topics.

A partial differential equation or PDE is any equation that contains a function and its partial derivatives. The general form of a PDE is:

$$F\left(x,y,\dots ,f(x,y),\dots \frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\dots \frac{{\partial }^{n}f}{\partial x^n},\frac{{\partial }^{n}f}{\partial y^n}\dots \right)=G(x,y,\dots )$$

To shorten derivatives, we often use the $\nabla$ symbol, where:

$${\nabla }^{2}u=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}$$

Some specific examples of PDEs include the heat equation:

$$\frac{\partial u}{\partial t}={\alpha }^2{\nabla }^{2}u$$

The wave equation:

$$\frac{{\partial }^{2}u}{\partial t^2}=c^2{\nabla }^{2}u$$

And Poisson’s equation:

$${\nabla }^{2}u=f$$

PDEs can have many, many solutions. For instance, take the PDE:

$$\frac{{\partial }^{2}u}{\partial x^2}=0$$

The general solution of this PDE is $u(x,y)=xf(y)+g(y)$. And yes, those can be any functions of $y$, whether they be $e^y$ or $y\sin (y)$ or $y+3y-\sqrt{y}$. Thus, just like ODEs need initial conditions to give unique solutions, PDEs need boundary conditions to give unique solutions. The typical set of boundary conditions for a PDE are values at the edges of the domain of the PDE (e.g. $u(x,t)\to 0$ as $x\to \pm \infty$). When the PDE is dependent on time, then the boundary condition for $t\to 0$ is often called an initial condition - remember, for PDEs, initial conditions are considered a type of boundary condition.

Analytically solving PDEs

Technique (ODE)

To solve a first order equation of the form $\frac{dy}{dt}=F(t,y)$ where $F(t,y)=f(t)g(y)$, we proceed as follows.

1. Multiply both sides by $dt$:

$$\begin{array}{c}dt\cdot \frac{dy}{dt}=f(t)g(y)\cdot dt\\ dy=f(t)g(y)dt\end{array}$$

2. Divide both sides by $g(y)$ to separate the variables:

$$\begin{array}{c}\frac{dy}{g(y)}=f(t)dt\end{array}$$

3. Because both sides contain a differential in their respective variables, integrate both sides and add a single constant of integration:

$$\begin{array}{c}\int \frac{dy}{g(y)}=\int f(t)dt+C\end{array}$$

4. Solve for $C$ using initial conditions. For example, if an initial condition where $y(t_0)=y_0$, substitute $t=t_0$ and $y=y_0$ to solve for $C$.

Examples (ODE)

Example 1

Consider the ordinary differential equation given by:

$$\frac{dy}{dt}=2t{e}^{-y}$$

We wish to solve for the general solution. Because $2t{e}^{-y}=2t\cdot e^{-y}$, the equation can be written as $\frac{dy}{dt}=f(t)g(y)$, where $f(t)=2t$ and $g(y)=e^{-y}$. Therefore, we can use separation of variables to solve it, as follows:

$$\begin{array}{rl}\frac{dy}{dt}=2t{e}^{-y}& \Rightarrow dy=2t{e}^{-y}dt\\ & \Rightarrow e^{y}dy=2tdt\\ & \Rightarrow \int e^{y}dy=\int 2tdt\\ & \Rightarrow e^y=t^2+C\\ & \Rightarrow y=\ln (t^2+C)\end{array}$$

Example 2

Consider the initial-value problem given by:

$$\frac{dy}{dt}=\frac{2y-1}{t+3},y(0)=-1$$

Using separation of variables, we proceed as follows:

$$\begin{array}{rl}\frac{dy}{dt}=\frac{2y-1}{t+3}& \Rightarrow dy=\frac{2y-1}{t+3}dt\\ & \Rightarrow \frac{dy}{2y-1}=\frac{dt}{t+3}\\ & \Rightarrow \int \frac{dy}{2y-1}=\int \frac{dt}{t+3}\\ & \Rightarrow \frac{1}{2}\ln |2y-1|=\ln |t+3|+C\end{array}$$

This is the general solution. We now use the initial condition to solve for $C$:

$$\begin{array}{rl}\frac{1}{2}\ln |2\cdot (-1)-1|=\ln |0+3|+C& \Rightarrow \frac{1}{2}\ln |-3|=\ln |3|+C\\ & \Rightarrow \frac{1}{2}\ln (3)=\ln (3)+C\\ & \Rightarrow -\frac{1}{2}\ln (3)=C\end{array}$$

Therefore, we find that:

$$\begin{array}{c}\frac{1}{2}\ln |2y-1|=\ln |t+3|-\frac{1}{2}\ln (3)\end{array}$$

We can rearrange to now solve for $y(t)$:

$$\begin{array}{rl}\frac{1}{2}\ln |2y-1|=\ln |t+3|+C& \Rightarrow \frac{1}{2}\ln |2y-1|=\ln |t+3|+-\frac{1}{2}\ln (3)\\ & \Rightarrow \ln |2y-1|=2\ln |t+3|-\ln (3)\\ & \Rightarrow \ln |2y-1|=\ln |(t+3{)}^2|-\ln (3)\\ & \Rightarrow \ln |2y-1|=\ln |\frac{(t+3{)}^2}{3}|\\ & \Rightarrow |2y-1|=|\frac{(t+3{)}^2}{3}|\\ & \Rightarrow 2y-1=\pm \frac{(t+3{)}^2}{3}\end{array}$$

The result gives two possible (implicit) solutions. However, we can use the initial condition to solve for the sign:

$$\begin{array}{rl}2y-1=\pm \frac{(t+3{)}^2}{3}& \Rightarrow 2\cdot (-1)-1=\pm \frac{(0+3{)}^2}{3}\\ & \Rightarrow -3=\pm \frac{3^2}{3}\\ & \Rightarrow -3=\pm 3\end{array}$$

This gives us a unique solution for $y(t)$:

$$\begin{array}{rl}2y-1=\pm \frac{(t+3{)}^2}{3}& \Rightarrow 2y-1=-\frac{(t+3{)}^2}{3}\\ & \Rightarrow 2y=\frac{(t+3{)}^2}{3}+1\\ & \Rightarrow y=\frac{(t+3{)}^2}{6}+\frac{1}{2}\end{array}$$

Technique (PDE)

Separation of variables in PDEs is used similarly as in ODEs. However, there are additional steps to the process. Consider the example initial-boundary value problem (IBVP):

$$\begin{array}{rl}{u}_t& =k{u}_{xx},\\ u(0,t)& =u(L,t)=0,\\ u(x,0)& =f(x)\end{array}$$

We can use separation of variables to determine separated equations and apply initial and boundary conditions to the resulting equations to find the solution.

1. Assume that the function for which we wish to solve is separable, that is, that $u(x,t)$ can be written in the following form:

$$\begin{array}{c}u(x,t)=\varphi (x)T(t)\end{array}$$

2. Solve for the derived terms using the separable equation:

$$\begin{array}{rl}u(x,t)=\varphi (x)T(t)& \Rightarrow u_t(x,t)=\varphi (x)T_t(t)\\ & \Rightarrow u_x(x,t)={\varphi }_x(x)T(t)\\ & \Rightarrow u_{xx}(x,t)={\varphi }_{xx}(x)T(t)\end{array}$$

3. Substitute the derived terms:

$$\begin{array}{c}\varphi (x)T_t(t)=k{\varphi }_{xx}(x)T(t)\end{array}$$

4. Divide both sides by $T(t)$:

$$\begin{array}{c}\varphi (x)\frac{T_t(t)}{T(t)}=k{\varphi }_{xx}(x)\end{array}$$

5. Divide both sides by $\varphi (x)$:

$$\begin{array}{c}\frac{T_t(t)}{T(t)}=k\frac{{\varphi }_{xx}(x)}{\varphi (x)}\end{array}$$

6. Because the sides are dependent on different variables, they are both equal a constant; if you increase $t$, the right side will still be the same. Therefore, let the equation equal a constant $-\lambda$:

$$\begin{array}{c}\frac{T_t(t)}{T(t)}=k\frac{{\varphi }_{xx}(x)}{\varphi (x)}=-\lambda \end{array}$$

7. Use each side to form a homogeneous equation:

$$\begin{array}{rl}\frac{{\varphi }_{xx}(x)}{\varphi (x)}=-\lambda & \Rightarrow {\varphi }_{xx}(x)=-\lambda \varphi (x)\\ & \Rightarrow {\varphi }_{xx}(x)+\lambda \varphi (x)=0\\ \frac{T_t(t)}{T(t)}=-\lambda & \Rightarrow T_t(t)=-\lambda T(t)\\ & \Rightarrow T_t+\lambda T(t)=0\end{array}$$

8. Use the initial conditions to solve for initial conditions for the new equations:

$$\begin{array}{rl}u(x,t)=\varphi (x)T(t)& \Rightarrow u(0,t)=\varphi (0)T(t)=0\end{array}$$

Now note that if $T(t)=0$, then $u(x,t)=0$, so $T(t)\ne 0$. Therefore, $\varphi (0)=0$. Likewise:

$$\begin{array}{rl}u(x,t)=\varphi (x)T(t)& \Rightarrow u(L,t)=\varphi (L)T(t)=0\end{array}$$

Using the same logic, $\varphi (L)=0$.

Example for the Helmholtz Equation

The Helmholtz Equation is an eigenvalue problem that occurred while using separation of variables for a PDE with a Laplacian operator ${\nabla }^2$ in a multi-dimensional, bounded region. However, for this guide, we will be using only 2 dimensions in Cartesian coordinates $(x,y)$ because it is easier to learn and display. It takes the form of:

$$\begin{array}{c}{\nabla }^2\varphi +\lambda \varphi =0\end{array}$$

It can have standard boundary conditions such as a Dirichlet boundary condition, $\varphi (x,y)=x_0$, or a Neumann boundary condition, ${\varphi }_x(x,y)=x_0$, but because it is used in functions with multiple variables, in 2-dimension, there can be the boundary condition:

$$\begin{array}{c}a\varphi +b\nabla \varphi \cdot \hat{n}=0\end{array}$$

where $a,b$ are functions of every element on the boundary, $\nabla$ is the gradient, and $\hat{n}$. Examining the condition, if $b=0$, then is becomes a Dirichlet boundary condition. If $a=0$, then the condition can be rewritten into $\nabla \varphi \cdot \hat{n}=0$. This translates into a Neumann boundary condition at every variable $\varphi$ depends on. Therefore, this boundary condition is a linear combination of Dirichlet and Neumann boundary conditions. According to the Stern-Liouville Theorem, - All eigenvalues are real - There are infinitely many eigenvalues with a smallest eigenvalue, but not a biggest eigenvalue - Each eigenvalue may have multiple eigen functions - The eigen functions for a complete set, i.e. piecewise, smooth functions can converge to linear combinations of eigenfunctions.

$$\begin{array}{c}f(x,y)\sim \sum _{\lambda }{a}_{\lambda }{\varphi }_{\lambda }(x,y)\end{array}$$

Eigenfunctions for different eigenvalues are orthogonal to each other. The eigenvalues are given by the Rayleigh quotient formula:

$$\begin{array}{c}\lambda =\frac{-\oint \varphi \nabla \varphi \cdot ds+\int {\int }_R|\nabla {\varphi }^2|dxdy}{\int {\int }_R{\varphi }^{2}dxdy}\end{array}$$

where $R$ is the boundaries of the region. Using the properties of this theorem and multi-variable calculus, you can solve for $\lambda$.

Example for the wave equation

Now, consider the IBVP of the wave equation, which is given by:

$$\frac{{\partial }^{2}u}{\partial t^2}=c^2{\nabla }^{2}u(x,y,t),0<x<L,0<y<H$$

with the initial and boundary conditions:

$$\begin{array}{rl}u(0,y,t)& =0,\\ u(L,y,t)& =0,\\ u(x,0,t)& =0,\\ u(x,H,t)& =0,\\ u(x,y,0)& =\alpha (x,y),\\ \frac{\partial u}{\partial t}u(x,y,0)& =\beta (x,y)\end{array}$$

Because the PDE is linear and homogeneous, $u(x,y,t)$ is separable. Therefore,

$$\begin{array}{c}u(x,y,t)=h(t)\varphi (x,y)\end{array}$$

Even though $\varphi (x,y)$ is separable because the region is a simple shape, let’s leave them together for later use. We can then compute the derivatives as follows:

$$\begin{array}{rl}u(x,y,t)=h(t)\varphi (x,y)& \Rightarrow \frac{\partial u}{\partial t}=\frac{dh}{dt}\varphi (x,y)\\ & \Rightarrow \frac{\partial u}{\partial x}=h(t)\frac{\partial \varphi }{\partial x}\\ & \Rightarrow \frac{{\partial }^{2}u}{\partial x^2}=h(t)\frac{{\partial }^2\varphi }{\partial x^2}\\ & \Rightarrow \frac{\partial u}{\partial y}=h(t)\frac{\partial \varphi }{\partial y}\\ & \Rightarrow \frac{{\partial }^{2}u}{\partial y^2}=h(t)\frac{{\partial }^2\varphi }{\partial y^2}\end{array}$$

Now, substituting the derivatives gives us:

$$\begin{array}{rl}\frac{{\partial }^{2}u}{\partial t^2}=c^2{\nabla }^{2}u(x,y,t)& \Rightarrow \frac{{\partial }^{2}u}{\partial t^2}=c^2(\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2})\\ & \Rightarrow \frac{dh}{dt}\varphi (x,y)=c^2(h(t)\frac{{\partial }^2\varphi }{\partial x^2}+h(t)\frac{{\partial }^2\varphi }{\partial y^2})\\ & \Rightarrow \frac{dh}{dt}\varphi (x,y)=c^{2}h(t)(\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2})\\ & \Rightarrow h^{\prime }(t)\varphi (x,y)=c^{2}h(t)({\varphi }_{xx}(x,y)+{\varphi }_{yy}(x,y))\\ & \Rightarrow \frac{h^{\prime }(t)}{c^{2}h(t)}\varphi (x,y)={\varphi }_{xx}(x,y)+{\varphi }_{yy}(x,y)\\ & \Rightarrow \frac{h^{\prime }(t)}{c^{2}h(t)}=\frac{{\varphi }_{xx}(x,y)}{\varphi (x,y)}+\frac{{\varphi }_{yy}(x,y)}{\varphi (x,y)}\\ & \Rightarrow \frac{h^{\prime }(t)}{c^{2}h(t)}=\frac{{\varphi }_{xx}(x,y)}{\varphi (x,y)}+\frac{{\varphi }_{yy}(x,y)}{\varphi (x,y)}=-\lambda \end{array}$$

This gives us two differential equations. First, we have the ODE for $h(t)$, given by:

$$\begin{array}{rl}\frac{h^{\prime }(t)}{c^{2}h(t)}=-\lambda & \Rightarrow h^{\prime }(t)=-\lambda c^{2}h(t)\end{array}$$

Second, we have the PDE for $\varphi (x,y)$:

$$\begin{array}{rl}\frac{{\varphi }_{xx}(x,y)}{\varphi (x,y)}+\frac{{\varphi }_{yy}(x,y)}{\varphi (x,y)}=-\lambda & \Rightarrow {\varphi }_{xx}(x,y)+{\varphi }_{yy}(x,y)=-\lambda \varphi (x,y)\end{array}$$

This is actually the Helmholtz equation in disguise, since it can be written as ${\nabla }^2\varphi +\lambda \varphi =0$. Now we can apply the boundary conditions:

$$\begin{array}{rl}u(x,y,t)=h(t)\varphi (x,y)& \Rightarrow u(0,y,t)=h(t)\varphi (0,y)=0\\ & \Rightarrow \varphi (0,y)=0\\ & \Rightarrow u(L,y,t)=h(t)\varphi (L,y)=0\\ & \Rightarrow \varphi (L,y)=0\\ & \Rightarrow u(x,0,t)=h(t)\varphi (x,0)=0\\ & \Rightarrow \varphi (x,0)=0\\ & \Rightarrow u(x,H,t)=h(t)\varphi (x,H)=0\\ & \Rightarrow \varphi (x,H)=0\end{array}$$

Therefore, the B.C.’s for $\varphi (x,y)$ reduce to:

$$\begin{array}{c}\varphi (0,y)=0,\varphi (L,y)=0,\varphi (x,0)=0,\varphi (x,H)=0\end{array}$$

Now, we can perform separation of variables on $\varphi$. Let $\varphi (x,y)=f(x)g(y)$. Then, the partial derivatives are given by:

$$\begin{array}{rl}\varphi (x,y)=f(x)g(y)& \Rightarrow {\varphi }_x(x,y)=f^{\prime }(x)g(y)\\ & \Rightarrow {\varphi }_{xx}(x,y)=f^{\prime \prime }(x)g(y)\\ & \Rightarrow {\varphi }_y(x,y)=f(x)g^{\prime }(y)\\ & \Rightarrow {\varphi }_{yy}(x,y)=f(x)g^{\prime \prime }(y)\end{array}$$

Substituting ${\varphi }_{xx}$ and ${\varphi }_{yy}$ into the previous homogeneous equations, we get:

$$\begin{array}{rl}{\varphi }_{xx}(x,y)+{\varphi }_{yy}(x,y)=-\lambda \varphi (x,y)& \Rightarrow f^{\prime \prime }(x)g(y)+f(x)g^{\prime \prime }(y)=-\lambda f(x)g(y)\\ & \Rightarrow \frac{f^{\prime \prime }(x)}{f(x)}g(y)+g^{\prime \prime }(y)=-\lambda g(y)\\ & \Rightarrow \frac{f^{\prime \prime }(x)}{f(x)}+\frac{g^{\prime \prime }(y)}{g(y)}=-\lambda \\ & \Rightarrow \frac{f^{\prime \prime }(x)}{f(x)}=-\frac{g^{\prime \prime }(y)}{g(y)}-\lambda \\ & \Rightarrow \frac{f^{\prime \prime }(x)}{f(x)}=-\frac{g^{\prime \prime }(y)}{g(y)}-\lambda =-\mu \end{array}$$

Where $\lambda$ is the separation constant. This gives us two differential equations for $f$ and $g$:

$$\begin{array}{rl}\frac{f^{\prime \prime }(x)}{f(x)}=-\mu & \Rightarrow f^{\prime \prime }(x)=-\mu f(x)\\ -\frac{g^{\prime \prime }(y)}{g(y)}-\lambda =-\mu & \Rightarrow g^{\prime \prime }(y)-\lambda g(y)=\mu g(y)\\ & \Rightarrow g^{\prime \prime }(y)=(\mu +\lambda )g(y)\end{array}$$

Therefore, the homogeneous equations in standard form are given by:

$$\begin{array}{c}{f}^{\prime \prime }(x)=-\mu f(x),g^{\prime \prime }(y)=(\mu +\lambda )g(y)\end{array}$$

Substituting in the B.C.’s gives us:

$$\begin{array}{rl}\varphi (x,y)=f(x)g(y)& \Rightarrow \varphi (0,y)=f(0)g(y)=0\\ & \Rightarrow f(0)=0\\ & \Rightarrow \varphi (L,y)=f(L)g(y)=0\\ & \Rightarrow f(L)=0\\ & \Rightarrow \varphi (x,0)=f(x)g(0)=0\\ & \Rightarrow g(0)=0\\ & \Rightarrow \varphi (x,H)=f(x)g(H)=0\\ & \Rightarrow g(H)=0\end{array}$$

The B.C.’s here reduce to:

$$\begin{array}{c}f(0)=0,f(L)=0,g(0)=0,g(H)=0\end{array}$$

Solving for $f(x)$, we have:

$$\begin{array}{rl}{f}^{\prime \prime }(x)=-\mu f(x)& \Rightarrow r^2{e}^{rx}=-\mu e^{rx}\\ & \Rightarrow r^2=-\mu \\ & \Rightarrow r=\pm \sqrt{-\mu }\end{array}$$

After solving through the cases of $\mu >0,\mu =0,$ and $\mu <0$, you will find $\mu =(\frac{n\pi }{L}{)}^2$ and $f(x)=B\sin (\frac{n\pi y}{L})$ where $n=1,2,\dots$ and $B$ is an arbitrary constant. Solving for $g(y)$ we get:

$$\begin{array}{rl}{g}^{\prime \prime }(y)=(\lambda +\mu )f(x)& \Rightarrow r^2{e}^{ry}=(\mu +\lambda )e^{ry}\\ & \Rightarrow r^2=\mu +\lambda \\ & \Rightarrow r=\pm \sqrt{\mu +\lambda }\end{array}$$

After solving the cases of $\lambda >\mu$, $\lambda =\mu$, and $\lambda <\mu$, you will find $\lambda =\mu +(\frac{m\pi }{H}{)}^2=(\frac{n\pi }{L}{)}^2+(\frac{m\pi }{L}{)}^2$ and $g(y)=D\sin (\frac{m\pi x}{H})$ where $D$ is arbitrary, $n=1,2,\dots$, $m=1,2,\dots$ Therefore, $\varphi (x,t)=f(x)g(y)=B\sin (\frac{n\pi y}{L})\sin (\frac{m\pi x}{H})$. Solving for $h(t)$, we have:

$$\begin{array}{rl}{h}^{\prime }(t)=-\lambda c^{2}h(t)& \Rightarrow \frac{h^{\prime }(t)}{h(t)}=-\lambda c^2\\ & \Rightarrow \int \frac{h^{\prime }(t)}{h(t)}dt=-\int \lambda c^{2}dt\\ & \Rightarrow \ln |h(t)|=-\lambda c^{2}t+C\\ & \Rightarrow h(t)=C{e}^{-\lambda c^{2}t}\end{array}$$

Afterwards, you can use superposition to solve for $u(x,y,t)$, giving a general solution of the form:

$$u(x,y,t)=\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }{B}_{mn}\sin (\frac{n\pi y}{L})\sin (\frac{m\pi x}{H})e^{-\lambda c^{2}t}$$

Numerically solving PDEs

The most common methods of solving PDEs numerically are the finite difference, finite element, and boundary element methods. We will start with the finite difference method (FDM). The basic idea is to transform a PDE into a linear algebra problem by encoding derivatives as matrices. Take one-dimensional wave equation as example. It is given by:

$$\frac{{\partial }^{2}u}{\partial t^2}=c^2\frac{{\partial }^{2}u}{\partial x^2}$$

We want to discretize the second partial derivative. Recall that the derivative is the limit of a difference quotient, which reads, in the discrete case, as:

$$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{u(x-h,t)+2u(x,t)+u(x+h)}{h^2}$$

Or, alternatively written:

$$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{u_{i-1,j}+2u_{i,j}+u_{i+1,j}}{h^2}$$

This works in all cases except at the left and right boundaries of the domain, as there is not another point to the left of the boundary to take the central difference. Instead, we use the single-sided difference - for the left, we have:

$$\frac{2u_{i,j}-5u_{i+1,j}+4u_{i+2,j}-u_{i+3,j}}{h^2}$$

And for the right, we have:

$$\frac{-u_{i-3,j}+4u_{i-2,j}-5u_{i-1,j}+2u_{i,j}}{h^2}$$

Note: see https://web.media.mit.edu/~crtaylor/calculator.html for a calculator for these values (enter 0, 1, 2, 3 for the left-handed 2nd derivative and -3, -2, -1, 0 for the right-handed 2nd derivative).

Now, partial derivatives are linear operators, just like matrices. So we can spatially discretize the equation by turning the 2nd spatial partial derivative into a matrix $A$ acting on the solution vector $U$:

$$\frac{{\partial }^{2}u}{\partial x^2}\approx AU$$

Where:

$$AU=\left[\begin{array}{ccccccc}2& -5& 4& -1& & & \\ 1& -2& 1& & & & \\ & 1& -2& 1& & & \\ & & \ddots & \ddots & \ddots & & \\ & & & 1& -2& 1& \\ & & -1& 4& -5& 2& \end{array}\right]\left[\begin{array}{c}{u}_{1,1}\\ u_{1,2}\\ ⋮\\ u_{1,n-1}\\ u_{1,n}\\ u_{2,1}\\ ⋮\\ u_{m,n}\end{array}\right]=\left[\begin{array}{c}2u_{1,1}-5u_{1,2}+4u_{1,3}-u_{1,4}\\ u_{1,3}-2u_{1,2}+u_{1,1}\\ u_{1,4}-2u_{1,3}+u_{1,2}\\ ⋮\\ u_{i,1}-5u_{i,2}+4u_{i,3}-u_{i,4}\\ u_{i,3}-2u_{i,2}+u_{i,1}\\ ⋮\\ u_{i,n}-2u_{i,n-1}+u_{i,n-1}\\ -u_{i,n-3}+4u_{i,n-2}-5u_{i,n-1}+2u_{i,n}\end{array}\right]$$

By discretizing the spatial derivative, the partial time derivative simply becomes an ordinary time derivative. Thus we can rewrite as:

$$U^{\prime \prime }=c^{2}AU$$

Which is an ordinary differential equation that can be solved with conventional ODE solvers. The only thing left is to compute an initial condition $U_0(x)$ and $U_0^{\prime }(x)$ (if you know the initial condition you can differentiate it to get the initial time derivative).

This method works well even in the higher-dimensional case. For instance, consider the generalized wave equation:

$$\frac{{\partial }^{2}u}{\partial t^2}=c^2{\nabla }^{2}u$$

By encoding the Laplacian ${\nabla }^{2}u$ as a matrix $A$ acting on $U$, the same approach can be used as with the other approaches previously mentioned. In fact, this approach works for all linear partial differential equations that have one time and one or several space derivatives.

For other cases of linear PDEs that are not time-based, such as Poisson’s or Laplace’s equation, they can be generally written in the form:

$$AU=B$$

And solved using standard linear algebra techniques for solving linear systems, of which there are many different solvers in popular libraries like NumPy, SciPy, Eigen, and others. Finally, for nonlinear PDEs, they can be generally written in the form:

$$AU-f(U)=0$$

Which is a root-finding problem that can be solved with Newton’s method:

$$U_{i+1}=U_i-J^{-1}U(x_0)$$

Solving vector differential equations

To solve vector differential equations, it is necessary to break up the vector differential equation into its components, and then solve for each of the components. For example, the first two of Maxwell’s equations result in 2 trivial equations. The second two are more complex, as the curl of a vector field produces another vector field. Thus, we must instead rewrite the curl component-by-component. The complete equations are:

$$\begin{array}{rl}\frac{\partial E_x}{\partial x}+\frac{\partial E_y}{\partial y}+\frac{\partial E_z}{\partial z}& =\frac{\rho }{{ϵ}_0}\\ \frac{\partial B_x}{\partial x}+\frac{\partial B_y}{\partial y}+\frac{\partial B_z}{\partial z}& =0\\ \frac{\partial E_z}{\partial y}-\frac{\partial E_y}{\partial z}& =-\frac{\partial B_x}{\partial t}\\ \frac{\partial E_x}{\partial z}-\frac{\partial E_z}{\partial x}& =-\frac{\partial B_y}{\partial t}\\ \frac{\partial E_y}{\partial x}-\frac{\partial E_x}{\partial y}& =-\frac{\partial B_z}{\partial t}\\ \frac{\partial B_z}{\partial y}-\frac{\partial B_y}{\partial z}& ={\mu }_0\left(J_x+{ϵ}_0\frac{\partial E_x}{\partial t}\right)\\ \frac{\partial B_x}{\partial z}-\frac{\partial B_z}{\partial x}& ={\mu }_0\left(J_y+{ϵ}_0\frac{\partial E_y}{\partial t}\right)\\ \frac{\partial B_y}{\partial x}-\frac{\partial B_x}{\partial y}& ={\mu }_0\left(J_z+{ϵ}_0\frac{\partial E_z}{\partial t}\right)\end{array}$$

So Maxwell’s equations result in 8 coupled PDEs. However, only 6 of them are independent - any system that satisfies the last 6 PDEs also must satisfy the first 2, so long as the boundary conditions satisfy the conservation of charge.

Here, note that for a solution, the current density must be set, which consists of one function each for $J_x,J_y,J_z$.

The Maxwell equations are essential in our understanding of electromagnetism, and our work heavily depends on being able to solve them (both analytically and numerically), allowing us to analyze phenomena as diverse as electromagnetic waves to electron beams. See Solving for the fields in a hollow waveguide for an example.

List of common PDEs

This list contains a number of very common PDEs in physics, mathematics, and the sciences for reference purposes.

Wave equation:

$$\frac{{\partial }^{2}u}{\partial t^2}=c^2{\nabla }^{2}u$$

Heat/diffusion equation:

$$\frac{\partial T}{\partial t}={\alpha }^2{\nabla }^{2}T$$

Laplace’s equation:

$${\nabla }^{2}u=0$$

Poisson’s equation:

$${\nabla }^{2}u=f(u)$$

Burger’s equation:

$$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=\nu \frac{{\partial }^{2}u}{\partial x^2}$$

Continuity equation:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot u=0$$

Navier-Stokes:

$$\nabla \cdot u=0$$ $$\rho \left(\frac{\partial u}{\partial t}+(u\cdot \nabla )u\right)=-\nabla p+\mu {\nabla }^{2}u+\rho g$$

Maxwell’s equations (in vector form):

$$\begin{array}{rl}\nabla \cdot \mathbf{E}& =\frac{\rho }{{ϵ}_0}\\ \nabla \cdot \mathbf{B}& =0\\ \nabla \times \mathbf{E}& =-\frac{\partial \mathbf{B}}{\partial t}\\ \nabla \times \mathbf{B}& ={\mu }_0\left(J+{ϵ}_0\frac{\partial \mathbf{E}}{\partial t}\right)\end{array}$$

Einstein field equations:

$$G_{\mu \nu }=\frac{8\pi G}{c^4}{T}_{\mu \nu }$$

Schrödinger equation:

$$i\hslash \frac{\partial \psi }{\partial t}=-\frac{{\hslash }^2}{2m}{\nabla }^2\psi +V(\mathbf{x})\psi$$

Euler-Lagrange equations:

$$\frac{\partial L}{\partial q}-\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}}\right)=0$$

Hamilton’s equations:

$$\frac{\partial H}{\partial p}=\dot{q},\frac{\partial H}{\partial q}=-\dot{p}$$

Black-Scholes equation:

$$\frac{\partial V}{\partial t}+\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^{2}V}{\partial S^2}+rS\frac{\partial V}{\partial S}-rV(t,S)=0$$

Sources and useful reading

• https://www.diva-portal.org/smash/get/diva2:1438836/FULLTEXT01.pdf
• https://youtu.be/hxGA1Je1P-s
• https://github.com/lukepolson/youtube_channel/blob/main/Python%20GPU/schrodinger.ipynb
• https://gist.github.com/c0rychu/6fc08cf5cdec2cc78a1f9334d103a869#file-schrodingereq_1d_tutorial-ipynb
• https://mathematica.stackexchange.com/questions/221207/finite-difference-method-for-1d-wave-equation
• https://nextjournal.com/sosiris-de/pde-2018
• https://empossible.net/wp-content/uploads/2019/08/Lecture-4a-Finite-Difference-Method.pdf
• https://aquaulb.github.io/book_solving_pde_mooc/solving_pde_mooc/notebooks/01_Introduction/01_00_Preface.html