# Solving Poisson-like PDE with FFT

Problem

I have an $$n\times n$$ grid, and each point on the grid is assigned two values: a score, and an (inverse) speed factor. There is a "turtle" moving along the grid, and it's goal is to maximize it's score. Each time it moves to a point $$x$$ it gets the score $$s_x$$ at that point, but it takes $$\epsilon_x$$ turns before it can move again. The turtle can only move in the four cardinal directions (north, east, south, west) and in the "center" direction (remaining where it currently is).

I want to assign a third value $$p_x$$, the "potential" at $$x$$, which helps the turtle find an optimal spot to move to. If the grid were extended to be continuous, the potential function should be continuous with a first derivative (or gradient).

Preliminaries

A natural choice for the potential function satisfies

$$\epsilon_xp_x = s_x + \frac{\sum\limits_{y\in \text{Neighbors}(x)}p_y}{4}\qquad (1)$$

(Notice that this almost exactly defines a qtable, except for the differing $$\epsilon_x$$ values.)

For now, assume that $$\epsilon_x$$ is constant. The equation happens to then be the same as a Discrete Poisson equation. This means that physically $$p_x$$ is the potential field in a system where each point has a $$-s_x$$ charge on it. There are many ways to solve the Discrete Poisson equation, but one of the fastest is to use spectral methods. See here for using the Discrete Cosine Transform (DCT) to solve it. Our grid can be rather large ($$50\times 50$$), so to fit within a one-second time constraint we need an $$O(n^2\log n)$$ algorithm.

However, $$\epsilon_x$$ is not constant. $$\epsilon_x$$ physically represents the electric permittivity at the point $$x$$, so if it were a one dimensional grid we could set every $$\epsilon_x$$ to $$1$$ and insert points to mimic having a higher permittivity. However, we are on a two dimensional grid. Not only that, but $$\epsilon_x$$ is not necessarily an integer—it could be any positive real number. My goal is to use spectral analysis to solve equation $$(1)$$ with an algorithm running in $$O(n^2\log n)$$ time.

Progress

First, rearrange the original equation:

$$\frac{p_{i+1,j}+p_{i-1,j}+p_{i,j+1}+p_{i,j01}}{4}-p_{i,j} + \gamma_{i,j}p_{i,j} = -s_{i,j}\qquad (2)$$ where $$\gamma_{i,j} = 1-\epsilon_{i,j}$$.

Notice that this is the finite difference equation for

$$\nabla^2 p + \gamma p = -s\qquad (3)$$

Now, $$\nabla^2(fp) = f\nabla^2p + 2\nabla f\cdot\nabla p + p\nabla^2f$$ so if I can find a function $$f$$ that satisfies $$2\nabla f\cdot\nabla p = 0, \frac{\nabla^2 f}{f} = \gamma$$ I would be mostly done. I just solve the Poisson equation $$\nabla^2(fp) = -fs$$ and divide by $$f$$. However, I don't know how to do this in $$O(n^2\log n)$$ time.

Alternatively, we can go straight to spectral analysis. Let $$p = \sum \hat{p}_{j,k}e^{i(jx+ky)}$$ $$s = \sum \hat{s}_{j,k}e^{i(jx+by)}$$

Plugging this into $$(3)$$ we get

$$\sum -\hat{p}_{j, k}(j^2+k^2)e^{i(jx+by)} + \gamma\sum \hat{p}_{j, k}e^{i(jx+ky)} = -\sum \hat{s}_{j, k}e^{i(jx+ky)}$$

Removing the summation yields $$\hat{p}_{j,k}(j^2+k^2-\gamma) = \hat{s}_{j, k}\implies \hat{p}_{j,k} = \frac{\hat{s}_{j,k}}{j^2+k^2-\gamma}$$

But $$\hat{p}_{j,k}, \hat{s}_{j,k}$$ are constants while $$\gamma$$ is not, so there is no solution.

Possible Duplicate?

• Spectral methods are best suited for linear problems with constant coefficients. If there is a spatially dependent coefficient then you could use spectral methods by expanding this coefficient in a Fourier series and then use convolutions but that's a pain. However, Eq. (2) is a five-diagonal linear system, and I thought for those the scaling of operation count is linear with the size, e.g., see discussion in 12000.org/my_notes/penta_diagonal_solver_in_matlab/index.htm Jul 14, 2021 at 5:06
• I think you lack boundary conditions.
– Bort
Jul 14, 2021 at 15:14
• I'm using Neumann boundary conditions. Jul 14, 2021 at 15:34
• @MaximUmansky that looks promising, but eq. (2) has diagonals that aren't all in one band. This answer has more info on solving this specific problem. I'll look into it. Jul 14, 2021 at 15:47
• I looked into solving outrigger systems, and for my particular equation it would run in $O(n^3)$ time. It might be fast enough, but I would prefer an $O(n^2\log n)$ algorithm if I can find one. Jul 14, 2021 at 16:07

You can use an iterative method to find $$p$$. Do the following steps:

1. Solve the equation $$\nabla^2 p^1 = -s$$.
2. Solve the new equation $$\nabla^2 p^{n+1} = -(s+\gamma p^n)$$.
3. Repeat step 2 until convergence.

Assuming you use the DCT method you linked to solve the Poisson equations, then steps 1 and 2 only take $$O(n^2\log n)$$ time. The convergence rate of step 3 is exponential, so very few iterations are needed.

The trickiest part is finding the convergence rate. Looking at the DCT method in that link, step 2 involves finding

$$\hat{p^{n+1}}_{x,y} = -\frac{\hat{f^n}_{x,y}}{2(C_x + C_y - 2)}$$

where $$\hat{f^n}_{x, y}$$ represents a coefficient in the DCT of $$s+\gamma p^n$$. Let $$r(x,y) = -2(C_x + C_y - 2)$$, and $$\hat{p^n} = \hat{p} + \delta^n$$ where $$\hat{p}$$ is the true solution. As the DCT is additive, and multiplication in real space corresponds to a convolution in Fourier space, we get

\begin{align*} \hat{p^{n+1}} &= \frac{\hat{s} + \hat{\gamma}*\hat{p^n}}{r}\\ &=\frac{\hat{s}+ \hat{\gamma}*\hat{p} + \hat{\gamma}*\delta^n}{r} \end{align*}

The true solution satisfies

$$\hat{p} = \frac{\hat{s}+ \hat{\gamma}*\hat{p}}{r}$$

so

\begin{align*} \delta^{n+1} &= \frac{\hat{s}+ \hat{\gamma}*\hat{p} + \hat{\gamma}*\delta^n}{r} - \frac{\hat{s}+ \hat{\gamma}*\hat{p}}{r}\\ &= \frac{\hat{\gamma}*\delta^n}{r} \end{align*}

Let $$||x||$$ represent the 1-norm of $$x$$, and $$m=\min(|r(x)|)$$. From Young's theorem,

$$\lVert r\delta^n\rVert = \lVert\hat{\gamma}*\delta^n\rVert \le \lVert\hat{\gamma}\rVert\lVert\delta^n\rVert\le \frac{\lVert\hat{\gamma}\rVert}{m}\lVert r\delta^n\rVert$$

Therefore the method converges exponentially. Unfortunately this implies that unless $$\hat{\gamma}$$ is not close to $$1$$, it won't converge.

• I noticed you are using fixed point iteration, which doesn't have the quickest convergence rate (or may not converge at all if $\lVert \hat{\gamma}\rVert > m$). You can speed this up using a quasi-Newton method. See Woodbury's matrix identity. Jul 16, 2021 at 2:56