Numerical solution to ${d\over dt} \nabla^2 p = {dp\over dx}$

Question

I need help numerically solving $$\frac{\partial \nabla^2 p} {\partial t} = \frac{\partial p} {\partial x}$$

I know that to solve

\begin{equation} \frac{\partial p} {\partial t} = \frac{\partial p} {\partial x} \end{equation}

the finite difference formula is \begin{equation} \frac{p_{i,j}^{n+1} - p_{i,j}^n} {\Delta t} = \frac{p_{i+1,j}^n - p_{i-1,j}^n} {2 \Delta x} \end{equation} This simplifies to \begin{equation} p_{i,j}^{n+1} = r \left( p_{i+1,j}^n - p_{i-1,j}^n \right) + p_{i,j}^n \end{equation} where $r = \frac{\Delta t} {2 \Delta x}$, which allows me to explicitly calculate $p$ at the next time step from $p$ at the current time step.

Following the same procedure for $\frac{\partial \nabla^2 p} {\partial t} = \frac{\partial p} {\partial x}$, I get

\begin{multline} p_{i-1,j}^{n+1} + p_{i+1,j}^{n+1} + p_{i,j-1}^{n+1} + p_{i,j+1}^{n+1} - 4p_{i,j}^{n+1} = \\ (1-r) p_{i-1,j}^n + (1+r) p_{i+1,j}^n + p_{i,j-1}^n + p_{i,j+1}^n - 4p_{i,j}^n \end{multline} where $r = \frac{\Delta \Delta t} {2}$ and $\Delta x = \Delta y = \Delta$. But now I must solve $p$ at five points simultaneously. I cannot simply march through the grid, updating each point as I go. How can I solve this equation cheaply and avoid inverting a matrix?

Answer 1

Here is a fast way to solve it. Take the 2D Fourier transform of the equation to get

$$\frac{\partial \hat{p}(\vec{k},t)}{\partial t} = -\frac{i k_x}{k^2}\hat{p}(\vec{k},t)$$

where $\hat{p}(\vec{k},t) = \mathcal{F}(p)$ is the Fourier transform of $p$. The solution of this equation is simply

$$\hat{p}(\vec{k},t) = u_0(\vec{k})e^{-\frac{i k_xt}{k^2}}$$

where I assume you have $u_0(\vec{k}) \equiv \mathcal{F}(p(x,y,0))$ as an initial condition at $t=0$. The solution in real space, at any time $t\geq 0$, can be found from

$$p(x,y,t) = \mathcal{F}^{-1}\left( u_0(\vec{k})e^{-\frac{i k_xt}{k^2}}\right)$$

where $\mathcal{F}^{-1}$ is the inverse Fourier transform. Numerically calculating the transforms is really fast and also easy to implement if you use, for example, the FFTW library. Note that if you only want the solution at one particular time-step then this method only requires two FFTs - one to get $u_0$ and one to get back to $p$ - no matter what time you are interested in (which is for example done in milliseconds for a $50\times 50$ grid).

Answer 2

Presumably you are using some boundary conditions such that $u$ uniquely determines $p$ (otherwise the problem is not well-posed). Since you didn't specify that, I'll assume homogeneous Dirichlet boundary conditions. Solve (in a discretized version) with these boundary conditions $\nabla^2 p = e_{ij}$ for each grid point $(i,j)$, where $e_{ij}$ is $1$ at $(i,j)$ and $0$ everywhere else, obtaining the Green's functions $g(x,y,i,j)$. Then given $u$, the solution of $\nabla^2 p = u$ is $p(x,y) = \sum_{i,j} u(i,j) g(x,y,i,j)$. Now given $p^n$, let $u^n = \nabla^2 p^n$, $u^{n+1} = u^n + (\Delta t) \partial_x p^n$, and $p^{n+1} = \sum_{i,j} u^{n+1}(i,j) g(\cdot, \cdot, i, j)$.