Solving Poisson-like PDE with FFT

Question

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?

Answer 1

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.